```
#\$r002
PHOTON USE
AUTOPLOT
file
PHI 5

cl;d 1 h1;d 1 ha;col3 1;blb4 2;redr
msg    temperature profile; press  to continue
ENDUSE
TEXT(1D RADIATION+HEAT SOURCE IN A SLAB:122
#cls
TITLE
DISPLAY
The problem considered is that of radiative heat transfer
in a stationary emitting and absorbing gray medium containing
a uniformly-distributed volumetric heat source. The medium is
bounded by two infinite, plane parallel walls at y=0 and y=L,
which are kept at the same uniform temperatures. Thus,
symmetry is exploited in the calculations. The energy transfer
is by pure radiation, so that the energy equation is simply:

-d/dy (Qr) + Qv = 0

where Qr is the radiative heat flux and Qv is the uniform
volumetric heat generation rate in the medium. Thermal radiation
is modelled by solving the equation for the radiosity R, as
follows:
d/dy ( 1/(a+s) dR/dy ) + 4*a*(E - R) = 0

where a is the absorption coefficient, s is the scattering
coefficient, and E is the black-body emissive power.
#pause
The problem is to solve for the temperature distribution given
the wall temperature Tw, the optical thickness Kr*L; and the wall
emissivity emw. Kr is the Rosseland mean absorption coefficient
which is given by: Kr=(a+s).

This problem has been solved by Deissler [ASME J.Heat Transfer
P241-246, (1964)] who used the Diffusion approximation with
jump boundary conditions to obtain the following solutions:

Qw  = 0.5*L*Qv

(Egc - Ew)/Qw = (3./16)*Kr*L+1./emw-0.5+3./(4.*Kr*L)

(Egw - Ew)/Qw = (1./emw-0.5+3.0/(4.*Kr*L))

Eg  = Egc - (3./8)*Kr*Qv*y**2

where Qw is the wall heat flux, and Egc and Egw are the gas
emissive powers at the centre line and wall.
#pause
ENDDIS
REAL(EMWN,TWN,EWN,TA,GY,EG,EGC,EGW,ACON)
CHAR(CH1);INTEGER(JJM1)
MESG( Enter optical thickness  0 < Kr*L << 10.0 (default 5)
EMWN=1.0;TWN=1000.0
** analytical solution
QWNA=QVOL*0.5*LENGTH
EGC=EWN+QWNA*(3.*OTHICK/16.+1./EMWN-0.5+3./(4.*OTHICK))
EGW=EWN+QWNA*(1./EMWN-0.5+3./(4.*OTHICK))
TGCA=(EGC/SIGMA)**0.25;TGWA=(EGW/SIGMA)**0.25
ACON=QWNA/(EGC-EWN)
GROUP 3,4,5. X,Y,Z-direction grid specification
GRDPWR(Y,50,0.5*LENGTH,1.0)
GROUP 6. Body-fitted coordinates or grid distortion
GROUP 7. Variables stored, solved & named
CP1=1.0
MESG( Enter required energy variable ? (TEM1 or H1)
IF(:CH1:.EQ.TEM1) THEN
+ MESG( TEM1 solution selected
ELSE
+ MESG( H1 solution selected
+ TMP1=LINH;TMP1B=1.0/CP1
ENDIF
GROUP 8. Terms (in differential equations) & devices
** Deactivate conduction & any built-in sources
TERMS(:CH1:,N,N,N,N,P,P)
GROUP 9. Properties of the medium (or media)
GROUP 10. Inter-phase-transfer processes and properties
GROUP 11. Initialization of variable or porosity fields
** analytical solution
DO JJ=1,NY
+PATCH(IN:JJ:,INIVAL,1,NX,JJ,JJ,1,NZ,1,1)
+GY=0.5*YFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1;GY=YFRAC(JJM1)+0.5*(YFRAC(JJ)-YFRAC(JJM1))
ENDIF
+GY=GY*YVLAST;EG=EGC-ACON*GY*GY
+TA=(EG/SIGMA)**0.25;INIT(IN:JJ:,HA,ZERO,TA)
ENDDO
GROUP 13. Boundary conditions and special sources
** Net radiation flux from wall
PATCH(WALLRB,NORTH,1,NX,NY,NY,1,NZ,1,1)
** uniformly-distributed volumetric heat source
PATCH(QHEAT,VOLUME,1,NX,1,NY,1,NZ,1,1)
COVAL(QHEAT,:CH1:,FIXFLU,QVOL)
GROUP 15. Termination of sweeps
LSWEEP=500;LITC=5
GROUP 16. Termination of iterations