PHOTON USE
   AUTOPLOT
   file
   PHI 5
 
   cl;d 1 h1;d 1 ha;col3 1;blb4 2;redr
   msg    temperature profile; press  to continue
   pause
   cl;d 1 radx;d 1 ra;col3 1;blb4 2;redr
   msg    radiation-flux profile; press e to end
   pause;end
   ENDUSE
 
    GROUP 1. Run title and other preliminaries
 
TEXT(Radiative Transfer + Heat Source  
TITLE
 
  DISPLAY
   The problem considered is that of pure radiative heat
   transfer in a 1d plane-parallel slab containing a uniformly
   distributed heat source. The medium may absorb, emit and
   scatter radiation and the boundaries of the slab are diffuse
   emitters and reflecters kept at the same uniform temperature.
   The medium is gray and because the wall temperatures are
   equal, symmetry can be exploited.
  ENDDIS
 
   Since the energy transfer is by pure radiation the
   energy equation is given by:
 
        -d/dx (Qrad) + Qvol = 0
 
   where Qrad is the radiative heat flux and Qvol is the uniform
   volumetric heat generation rate in the medium. The equation
   for the composite radiative heat flux is given by:
 
        d/dx ( 1/(a+s) d/dx (Rx) ) + a (E - Rx) = 0
 
   where a is the absorption coefficient, s is the scattering
   coefficient, Rx is the composite radiation flux defined as
   the average of the +ve and -ve radiation fluxes, and E is
   the black-body emissive power. It may be noted that the
   radiative heat flux is given by:
 
        d/dx (Qrad) = 2a (E - Rx)
 
   The black-body emissive power E=sig*T**4  where sig is the
   Stefan-Boltzmann constant and T is the temperature of the
   medium. The problem is the determination of the temperature
   and composite radiative-flux distributions, as given by the
   following analytical solutions:
 
     Rx = Ew + 0.5*Qvol*L* [ 2/emw-1+0.5*(a+s)*L*{1-(x/L)**2} ]
 
     E  = Ew + Qvol*L*[ 1./(2*a*L) + 1./emw - 0.5
 
                      + 0.25*(a+s)*L*{1 - (x/L)**2} ]
 
    where Ew is the emissive power at the wall, L is the slab
    width from symmetry plane to wall, and emw is the emissivity
    of the wall. For the case considered below, Qvol is taken to
    be Qvol = Ew/L/(1./emw - 0.5).
 
    The locally-defined parameters are as follows:
       GSIGMA         Stefan-Boltzmann constant    {W/m**2/K**4}
       SCAT           Scattering coefficient       {1/m}
       ABSORB         Absorption coefficient       {1/m}
       EMIW           emissivity of the wall
       TWAL           hot wall temperature         { K }
       QVOL           volumetric internal heat source {W/m**3}
 
CHAR(CH1);REAL(GSIGMA,SCAT,ABSORB,EMIW,TWAL,QVOL)
GSIGMA=5.6697E-8;SCAT=0.5;ABSORB=0.5;EMIW=1.0;TWAL=1000.0
    GROUPS 3,4,5. X,Y,Z-direction grid specification
REAL(LENGTH);LENGTH=1.0;GRDPWR(X,50,LENGTH,1.0)
    GROUP 7. Variables stored, solved & named
    
    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(FLUX,ABSORB,SCAT,:CH1:);STORE(EMPO)

  CP1=1.0
  MESG( Enter required energy variable ? (TEM1 or H1)
  IF(:CH1:.EQ.' ') THEN
   READVDU(CH1,CHAR,H1)
  ENDIF  
  IF(:CH1:.EQ.TEM1) THEN
  + MESG( TEM1 solution selected
  ELSE
  + MESG( H1 solution selected
  + TMP1=LINH;TMP1B=1.0/CP1
  ENDIF
RADIAT(FLUX,ABSORB,SCAT,:CH1:);STORE(EMPO)
    GROUP 8. Terms (in differential equations) & devices
   ** Deactive built-in sources, convection and conduction
TERMS(:CH1:,N,N,N,N,P,P)
    GROUP 11. Initialization of variable or porosity fields
FIINIT(RADX)=GSIGMA*TWAL**4;FIINIT(:CH1:)=TWAL
  ** Define Qvol = Ew/L/ (1./emw - 0.5)
QVOL=FIINIT(RADX)/XULAST/(1./EMIW-0.5)
  ** analytical solution
REAL(EW,EG,QRADA,ALF,BET,GAM,ACON,TA,RAN,GX,GXRAT)
STORE(HA,RA);INTEGER(JJM1)
EW=GSIGMA*TWAL**4;QRADA=2.*EW;ALF=2./EMIW-1.0;GAM=0.5*QVOL*XULAST
BET=0.5*(ABSORB+SCAT)*XULAST;ACON=1.0/(2.*ABSORB*XULAST)+0.5*ALF
DO JJ=1,NX
+PATCH(IN:JJ:,INIVAL,JJ,JJ,1,NY,1,NZ,1,1)
+GX=0.5*XFRAC(JJ)
IF(JJ.NE.1) THEN
+JJM1=JJ-1;GX=XFRAC(JJM1)+0.5*(XFRAC(JJ)-XFRAC(JJM1))
ENDIF
+GX=GX*XULAST;GXRAT=(GX/XULAST)**2
+EG=EW+QVOL*XULAST*(ACON+0.5*BET*(1.-GXRAT))
+RAN=EW+GAM*(ALF+BET*(1.0-GXRAT));TA=(EG/GSIGMA)**0.25
+INIT(IN:JJ:,RA,ZERO,RAN);INIT(IN:JJ:,HA,ZERO,TA)
ENDDO
    GROUP 13. Boundary conditions and special sources
   ** Net radiation flux from wall
PATCH(WALLR,EAST,NX,NX,1,NY,1,NZ,1,1)
COVAL(WALLR,RADX,EMIW/(2.0-EMIW),GSIGMA*TWAL**4)
   ** 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=50
    GROUP 16. Termination of iterations
RESREF(:CH1:)=1.E-12*QVOL*LENGTH;RESREF(RADX)=0.5*RESREF(:CH1:)
    GROUP 17. Under-relaxation devices
RELAX(:CH1:,FALSDT,100./QVOL)
    GROUP 22. Spot-value print-out
IXMON=NX/2;NPLT=5;NXPRIN=NX/10
    GROUP 23. Field print-out and plot control
OUTPUT(:CH1:,Y,N,N,Y,Y,Y)
PATCH(XWISE,PROFIL,1,NX,1,NY,1,NZ,1,1)
PLOT(XWISE,:CH1:,0.0,0.0);PLOT(XWISE,RADX:,0.0,0.0)
    GROUP 24. Dumps for restarts
TSTSWP=-1