#$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)
REAL(QRAD,QVOL,KRAD,OTHICK,LENGTH,QRADA,TWNA,QWNA,TGCA,TGWA)
CHAR(CH1);INTEGER(JJM1)
MESG( Enter optical thickness  0 < Kr*L << 10.0 (default 5)
READVDU(OTHICK,REAL,5.0)
LENGTH=1.0;KRAD=OTHICK/LENGTH
SCATT=0.0;EMISS=KRAD-SCATT
EMWN=1.0;TWN=1000.0
EWN=SIGMA*TWN**4;QRAD=EWN;QVOL=EWN/LENGTH
  ** 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)
READVDU(CH1,CHAR,H1)
IF(:CH1:.EQ.TEM1) THEN
+ MESG( TEM1 solution selected
ELSE
+ MESG( H1 solution selected
+ TMP1=LINH;TMP1B=1.0/CP1
ENDIF
RADIAT(RADI,EMISS,SCATT,:CH1:);STORE(EMPO)
    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
FIINIT(:CH1:)=TWN;FIINIT(SRAD)=EWN
  ** analytical solution
STORE(HA);ACON=3.0*KRAD*QVOL/8.0
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)
COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN)
   ** 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
SELREF=F;RESREF(:CH1:)=1.E-4*QRAD;RESREF(SRAD)=RESREF(:CH1:)
    GROUP 17. Under-relaxation devices
RELAX(:CH1:,LINRLX,0.3);RELAX(SRAD,LINRLX,0.6)
    GROUP 18. Limits on variables or increments to them
VARMIN(:CH1:)=0.5*TWN
    GROUP 19 to 21
    GROUP 22. Spot-value print-out
IYMON=1;NPLT=20;NYPRIN=1;ITABL=3
    GROUP 23. Field print-out and plot control
OUTPUT(:CH1:,Y,N,N,Y,Y,Y);PATCH(YWISE,PROFIL,1,1,1,NY,1,1,1,1)
PLOT(YWISE,:CH1:,0.0,0.0);PLOT(YWISE,SRAD,0.0,0.0)
    GROUP 24. Dumps for restarts
TSTSWP=-1