#$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 RADIATIVE EQUILIBRIUM IN SLAB : 121
#cls
TITLE
  DISPLAY
   The problem considered is radiative heat transfer in a stationary
   emitting and absorbing gray medium bounded by two infinite, plane
   parallel walls. The boundary surfaces at y=0 and y=L are kept at
   fixed temperatures Ts and Tn. The energy transfer is by pure
   radiation, so that the energy equation is simply:
 
                  -d/dy (Qr) = 0
 
   where Qr is the radiative heat transfer per unit area. Thermal
   radiation is modelled by solving an 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.
 
   The problem is to solve for the temperature distribution, and
   hence Qr, given Ts, Tn, the optical thickness Kr*L; and the wall
   emissivities emws and emwn. Kr is the Rosseland mean absorption
   coefficient which is given by: Kr=(a+s).
#pause
   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:
 
     Qr = (Ews - Ewn)/(0.75*Kr*L + 1./emws + 1./emwn - 1.0)
 
     (Ews - Egs)/Qrad = 1./emws - 0.5
 
     (Egn - Ewn)/Qrad = 1./emwn - 0.5
 
     E  = Egs - 0.75*Kr*Qr*y
 
   where Ews and Ewn are the emissive powers at the wall, and Egs
   and Egn are the emissive powers in the gas at the wall. The
   results of the radiosity model show excellent agreement with
   Deissler's results for the whole range of optical thickness
   considered by Deissler, i.e. 0 < Kr*L < 10.0.
  ENDDIS
#pause
REAL(EMWS,EMWN,TWS,TWN,EWS,EWN,ACON,GY)
REAL(QRAD,KRAD,OTHICK,LENGTH,QRADA,TGNA,TGSA,EGS,EG,TA)
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
EMWS=1.0;EMWN=1.0;TWS=1500.0;TWN=1000.0
EWS=SIGMA*TWS**4;EWN=SIGMA*TWN**4;QRAD=EWS-EWN
  ** analytical solution
QRADA=1.0/(0.75*OTHICK+1./EMWS+1./EMWN-1.0)
QRADA=QRADA*(EWS-EWN)
TGSA =((EWS-QRADA*(1./EMWS-0.5))/SIGMA)**0.25
TGNA =((EWN+QRADA*(1./EMWN-0.5))/SIGMA)**0.25
    GROUP 3,4,5. X,Y,Z-direction grid specification
GRDPWR(Y,50,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
   ** Deactive 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:)=TWS;FIINIT(SRAD)=EWS
  ** analytical solution
STORE(HA);ACON=0.75*KRAD*QRADA;EGS=SIGMA*TGSA**4
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=EGS-ACON*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(WALLRA,SOUTH,1,NX,1,1,1,NZ,1,1)
COVAL(WALLRA,SRAD,2.*EMWS/(2.0-EMWS),EWS)
PATCH(WALLRB,NORTH,1,NX,NY,NY,1,NZ,1,1)
COVAL(WALLRB,SRAD,2.*EMWN/(2.0-EMWN),EWN)
    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*TWS
    GROUPS 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