TALK=T;RUN( 1, 1) r221

     **** PIL fragment for IMMERSOL method
REAL(SIGMA); SIGMA = 5.6697E-8
STORE(EMIS,SCAT)
EMISS=0.0;SCATT=0.0;FIINIT(EMIS)=EMISS;FIINIT(SCAT)=SCATT
SOLVE(T3);RELAX(T3,LINRLX,0.5);STORE(WGAP)
DISWAL
TERMS(T3,N,N,Y,N,Y,N)
char(temps,direct,declare);temps=$r198;direct=$r199;declare=$r200
char(geom1,geom2,radflux);geom1=$r195;geom2=$r196;radflux=$r197
char(geom3);geom3=$r194
MESG(IMMERSOL settings have been loaded
 ** LOAD(x  1) from the x Input Library
 ** LOAD(x  1) from the x Input Library
 
  ** declarations needed for many immersol examples **
REAL(QRAD,WSL1,WSL2,WGAP,THOT,TCLD,EMISH,EMISC,EMISB)
REAL(OPTHI,KROSS,TGHT,TGCL,KGAS,NN1,NN2,AA1,AA2)
REAL(HHOT,HCLD,CP111,CP112)
REAL(SIZX,SIZZ,SIZY,UIN)
REAL(QCNV,AEMISH,AEMISC,NNH,NNC,TGH,TGC)
REAL(LENG1,LENG2,GAP,WPLT)
REAL(EMISM1,EMISM2,EMISW1,EMISW2,REYNO)
REAL(SPH0,TAIR,TFUEL,TWALL,HAIRIN,CPAIR,CPPR,DTF)
REAL(CPFU,WAIR,WFU,WPR,GRHO,STOIC,FSTOI,HFU,TKEI1,EMISI1,FINF,HINF)
 ** LOAD(x200) from the x Input Library
 ** LOAD(x200) from the x Input Library
   IMMERSOL Emitting and absorbing gray media
            (1D, TEM1)
TEXT(IMMERSOL 1D Radiative equilibrium   :202
TITLE
  **************************************************************
   PHOTON USE
   AUTOPLOT
   FILE
   PHI 5
 
   CL;DA 1 TEM1;DA 1 T3;COL3 1;BLB4 2;REDR
   MSG    Temperature (curve) and T3 (crosses) profiles
   ENDUSE
  **************************************************************
#cls
  DISPLAY
   The problem is as for case r202; but 
   1. the solid walls at the boundaries are both thicker and divided
      into five cells each;
      
   2. The heat conduction due to true (ie not T3) temperature 
      gradients, which was switched off in r202, is swiched back on
      again;
      
   3. In-Form "stored" statements are used for eliciting print-out
      of coefficients, which is useful for checking correctness. 
      
   4. Also printed is the variable #3-1, which is the value of 
      t3-tem1. This should ideally be zero in the solid regions; 
      but differences of a few degrees will be noted. Increasing 
      the value of RSG41 to above 2.0 should diminish these 
      differences.
  ENDDIS
  
  ** sequence used in several immersol examples
INTEGER(IDIR,I1,I2,J1,J2,K1,K3)
MESG(  The problem may be stated and solved in any direction
MESG(  ENTER  1 for X; 2 (or RETURN) for Y, 3 for Z.
READVDU(IDIR,INT,2)
MESG(  ENTER Temperature of hot solid (default 800K):
READVDU(THOT,REAL,800.)
MESG(  ENTER Temperature of cold solid (default 400K):
READVDU(TCLD,REAL,400.)
MESG(  ENTER Surface emissivity of hot solid (default 0.3):
READVDU(EMISH,REAL,0.3)
MESG(  ENTER Surface emissivity of cold solid (default 0.1):
READVDU(EMISC,REAL,0.1)
 ** LOAD(x198) from the x Input Library
 ** LOAD(x198) from the x Input Library
MESG(  ENTER optical thickness Kr*Wg (default 1.):
READVDU(OPTHI,REAL,1.0)
  settings of wall thicknesses and the gap between them
WSL1 = 0.01;  WSL2= 0.01;  WGAP= 1.0;
  conductivity of medium swet to zero
KGAS= 0.0
  QRAD first used to store EH-EC (see above)
QRAD= SIGMA*(THOT**4 - TCLD**4)
  KROSS deduced from above-set optical thickness; scattering
  coefficient set to zero; then medium emissivity deduced
KROSS= OPTHI/WGAP;  SCATT= 0.0;  EMISS= KROSS-SCATT
MESG( Expected radiative heat flux (W/m**2):
  QRAD now computed, based on above settings.
  Note that KGAS has been set to zero, so that NN1 and NN2 = 0
NN1 = KGAS*KROSS/4./SIGMA/THOT**3
NN2 = KGAS*KROSS/4./SIGMA/TCLD**3
QRAD= QRAD + KGAS/WGAP*(THOT-TCLD)*0.75*OPTHI
AA1 = (1./EMISH-0.5)/(1.+0.75*NN1)
AA2 = (1./EMISC-0.5)/(1.+0.75*NN2)
QRAD= QRAD/(0.75*OPTHI + AA1 + AA2)
QRAD
MESG( Expected gas temperatures near walls are:
TGCL = ((SIGMA*TCLD**4 + QRAD*AA2)/SIGMA)**0.25
TGHT = ((SIGMA*THOT**4 - QRAD*AA1)/SIGMA)**0.25
TGCL
TGHT
#pause
SOLVE(TEM1)
STORE(PRPS)
  ** geom1 **

IF(IDIR.EQ.1) THEN
 NREGX= 3
 IREGX= 1; GRDPWR(X,5,0.5*wgap,1.0)
 IREGX= 2; GRDPWR(X,10,WGAP,1.0)
 IREGX= 3; GRDPWR(X,5,0.5*wgap,1.0)
 PATCH(SOL1,INIVAL, #1,#1,1,NY,1,NZ,1,LSTEP)
 PATCH(SOL2,INIVAL,#3,#3,1,NY,1,NZ,1,LSTEP)
 PATCH(HOT,CELL,1, 1,1,NY,1,NZ,1,LSTEP)
 PATCH(COLD,CELL,nx,nx,1,NY,1,NZ,1,LSTEP)
ENDIF
IF(IDIR.EQ.2) THEN
 NREGY= 3; IREGY= 1; GRDPWR(Y,5,0.5*wgap,1.0)
 IREGY= 2; GRDPWR(Y,10,WGAP,1.0); 
 IREGY= 3; GRDPWR(Y,5,0.5*wgap,1.0)
 PATCH(SOL1,INIVAL,1,NX, #1,#1,1,NZ,1,LSTEP)
 PATCH(SOL2,INIVAL,1,NX,#3,#3,1,NZ,1,LSTEP)
 PATCH( HOT,CELL,1,NX, 1,1,1,NZ,1,LSTEP)
 PATCH(COLD,CELL,1,NX,ny,ny,1,NZ,1,LSTEP)
ENDIF
IF(IDIR.EQ.3) THEN
 NREGZ= 3; IREGZ= 1; GRDPWR(Z,5,0.5*wgap,1.0)
 IREGZ= 2; GRDPWR(Z,10,WGAP,1.0)
 IREGZ= 3; GRDPWR(Z,5,0.5*wgap,1.0)
 PATCH(SOL1,INIVAL,1,NX,1,NY, #1, #1,1,LSTEP)
 PATCH(SOL2,INIVAL,1,NX,1,NY,#3,#3,1,LSTEP)
 PATCH( HOT,CELL,1,NX,1,NY,1,1,1,LSTEP)
 PATCH(COLD,CELL,1,NX,1,NY,nz,nz,1,LSTEP)
ENDIF
 ** LOAD(x195) from the x Input Library
 ** LOAD(x195) from the x Input Library
  ** radflux **
if(nx.gt.1) then
 store(qrx)
endif
if(ny.gt.1) then
 store(qry)
endif
if(nz.gt.1) then
 store(qrz)
endif
 ** LOAD(x197) from the x Input Library
 ** LOAD(x197) from the x Input Library
TERMS(TEM1,N,N,y,N,Y,N);  PRNDTL(TEM1)=1.e10
INIADD= F;  FIINIT(TEM1)= (THOT+TCLD)/2.;  FIINIT(PRPS)= 0.
FIINIT(T3)=FIINIT(TEM1)
 COVAL(HOT,TEM1,FIXVAL,THOT); COVAL(COLD,TEM1,FIXVAL,TCLD)
 COVAL(HOT,TEM1,1.e5,THOT); COVAL(COLD,TEM1,1.e5,TCLD)
 INIT(SOL1,PRPS,0.0,111.); INIT(SOL2,PRPS,0.0,112.)
LSWEEP= 500; SELREF=F; RESREF(T3)= 1.E-3
IXMON=NX/2+1;IYMON=NY/2+1;IZMON=NZ/2+1;
NXPRIN=1; NYPRIN=1; NZPRIN=1
OUTPUT(LTLS,N,N,N,N,N,N); OUTPUT(WDIS,N,N,N,N,N,N)
 
  load macro for storing and setting emissivity and scattering
  coefficient, using EMISS, SCATT, EMISH and EMISC
 
  use of 3D store for emissivity,
  over-writing what may be supplied by SPEDAT
STORE(EMIS,SCAT)
FIINIT(EMIS)=EMISS;FIINIT(SCAT)=SCATT
IF(IDIR.EQ.1) THEN
 PATCH( HOTINI,INIVAL, #1, #1,1,NY,1,NZ,1,LSTEP)
 PATCH(COLDINI,INIVAL,#3,#3,1,NY,1,NZ,1,LSTEP)
  line-printer output
 PATCH(PROFILES,PROFIL,1,NX,1,1,1,1,1,1)
 COVAL(PROFILES,QRX,0.0,0.0)
ENDIF
IF(IDIR.EQ.2) THEN
 PATCH( HOTINI,INIVAL, 1, NX,#1,#1,1,NZ,1,LSTEP)
 PATCH(COLDINI,INIVAL,1,NX,#3,#3,1,NZ,1,LSTEP)
  line-printer output
 PATCH(PROFILES,PROFIL,1,1,1,NY,1,1,1,1)
 COVAL(PROFILES,QRY,0.0,0.0)
ENDIF
IF(IDIR.EQ.3) THEN
 PATCH( HOTINI,INIVAL, 1, NX,1,NY,#1,#1,1,LSTEP)
 PATCH(COLDINI,INIVAL,1,NX,1,NY,#3,#3,1,LSTEP)
  line-printer output
 PATCH(PROFILES,PROFIL,1,1,1,1,1,NZ,1,1)
 COVAL(PROFILES,QRZ,0.0,0.0)
ENDIF
 
  setting absorptivity/emissivity of solids
COVAL(HOTINI,EMIS,0.0,EMISH)
COVAL(COLDINI,EMIS,0.0,EMISC)
 
COVAL(PROFILES,T3,0.0,0.0)
IF(LIBREF.NE.201.AND.LIBREF.NE.204) THEN
  COVAL(PROFILES,TEM1,0.0,0.0)
ENDIF
 ** LOAD(x193) from the x Input Library
 ** LOAD(x193) from the x Input Library
 
VARMAX(T3)=THOT;VARMIN(T3)=TCLD
VARMAX(TEM1)=THOT;VARMIN(TEM1)=TCLD
  
  
  inform7begin
(stored var #3-1 is t3-tem1)
(stored var #rs1 is resi(tem1))
(stored var #rs3 is resi(t3))
(stored var #ap1 is apco(tem1))
(stored var #ap3 is apco(t3))
if(idir.eq.1) then
(stored var #ae1 is aeco(tem1))
(stored var #ae3 is aeco(t3))
endif 
if(idir.eq.2) then
(stored var #an1 is anco(tem1))
(stored var #an3 is anco(t3))
endif 
if(idir.eq.3) then
(stored var #ah1 is ahco(tem1))
(stored var #ah3 is ahco(t3))
endif 
  inform7end
  
  
  
STOP