#immersol)
#declare
#cls
  DISPLAY
  IMMERSOL  2D radiative heat exchange in a duct flow with
            fixed wall temperatures.

   Laminar or turbulent; H1 - T3 solved. Two solid plates (might
   be of different materials) in the middle. The emissivities of
   the duct walls, and of each solid can be set to different values.

   Chemical reaction in the gas is simulated by the SCRS model.
   
  ENDDIS
#pause   
  **************************************************************
BOOLEAN(LTURB); LTURB= T
  **************************************************************
  PHOTON USE
   p ; ; ; ; ;
 
   msg Computational Domain:
   gr k 1
   use patgeo
   msg Press Any Key to Continue...
   pause
   cl
   set vec av off
   msg Velocity Vectors:
   vec k 1 sh
   msg Press Any Key to Continue...
   pause
   cl
   msg Contours of Pressure:
   con p1 k 1 fi;0.005
   pause
   cl
   msg Contours of TMP1:
   con tmp1 k 1 fi;0.005
   pause
   cl
   msg Contours of T3:
   con t3 k 1 fi;0.005
   msg Press E  to exit PHOTON ...
  ENDUSE
  **************************************************************
TEXT(2D radiative heat exchange, + SCRS  :211
TITLE
INTEGER(MAT1,MAT2)
LENG1= 0.2;  LENG2= 0.6;  GAP= 0.1;  WPLT= 0.01
    * Stoichiometic ratio and heat of reaction:
STOIC= 17.24;  FSTOI= 1./(1.+STOIC);  HFU= 4.9E7
    * Specific heats:
        CPAIR= 1.5E3;  CPPR= 1.5E3;  CPFU= 1.5E3;  GRHO= 3.606
CPAIR= 1005;  CPPR= 1005;  CPFU= 3800
        CPAIR= 1005;  CPPR= 1005;  CPFU= 1005
    * Molecular weights:
WAIR= 29.0;  WFU= 16.0;  WPR= 28.0
    * Flow parameters:
TAIR = 300.0;       TFUEL= 300.0;  TWALL= 1000.
HAIRIN= CPAIR*TAIR;  FINF = FSTOI;  HINF = CPAIR*TFUEL+FINF*HFU
IF(LTURB) THEN
+ UIN  = 10.0
+ EGWF = T;  WALLCO= GRND2
ELSE
+ REYNO = 200.;  UIN= 1.0;  ENUL= UIN*GAP/REYNO
+ WALLCO= 1.0
ENDIF
    * Define emissivities of the domain walls:
REAL(EMIW1, EMIW2, EMIP1, EMIP2)    
EMIW1= 0.9;  EMIW2= 0.9
    * Define material and emissivity of its surface for the
    * 1st plate:
MAT1= 111;  EMIP1= 0.8
    * Define material and emissivity of its surface for the
    * 2nd plate:
MAT2 =112;  EMIP2= 0.1
    * Define optical thickness for the gap between plates:
OPTHI= 0.1;    KROSS= OPTHI/GAP;  SCATT= 0.0;  EMISS= KROSS-SCATT
#geom3
    GROUP 7. Variables stored, solved & named
NAME(C1)= MIXF;  SOLVE(P1,U1,V1,H1,FUEL,MIXF)
STORE(PRPS,OXID,PROD,TMP1,H0_1,SPH1,KOND,DEN1)
#radflux
IF(LTURB) THEN
+TURMOD(KEMODL); STORE(GEN1,ENUT)
ENDIF
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,Y,Y,N,Y,N)
    GROUP 9. Properties of the medium (or media)
SETPRPS(1,0)
RHO1=GRND6; RHO1A=WFU; RHO1B=WAIR; RHO1C=WPR
TMP1= SCRSNONEQ
TMP2A= FSTOI;  TMP2B= HFU;  PRESS0= 1.0E5
CP1=GRND10; CP1A= CPFU; CP1B= CPAIR; CP1C= CPPR
GRHO = PRESS0/((8314/WAIR)*TAIR)
    GROUP 11. Initialization of variable or porosity fields
INIADD= F;  FIINIT(U1)= UIN
FIINIT(MIXF)= FINF;FIINIT(FUEL)= FINF
FIINIT(H1)= HAIRIN+FIINIT(FUEL)*HFU
IF(LTURB) THEN
+ FIINIT(KE)= 0.1125*FIINIT(U1)*FIINIT(U1)
+ FIINIT(EP)= 0.1643*FIINIT(KE)**1.5/(0.27*YVLAST)
ENDIF
FIINIT(EMIS)= EMISS; FIINIT(SCAT)= SCATT
    *** Solid plates:
PATCH(PLT1,INIVAL,#2,#2,#2,#2,1,NZ,1,LSTEP)
 INIT(PLT1,PRPS,0.0,MAT1)
 INIT(PLT1,EMIS,0.0,EMIP1)
 INIT(PLT1,H1,  0.0,226*TAIR)
PATCH(PLT2,INIVAL,#2,#2,#4,#4,1,NZ,1,LSTEP)
 INIT(PLT2,PRPS,0.0,MAT2)
 INIT(PLT2,EMIS,0.0,EMIP2)
 INIT(PLT2,H1,  0.0,473*TAIR)
    GROUP 13. Boundary conditions and special sources
    *** Inlet:
PATCH(IN1,WEST,$1,$1,#1,#NREGY,1,NZ,1,LSTEP)
 COVAL(IN1,P1,FIXFLU,GRHO*UIN); COVAL(IN1,H1,ONLYMS,HINF)
 COVAL(IN1,  U1,ONLYMS, UIN); COVAL(IN1,  V1,ONLYMS, 0.0)
 COVAL(IN1,MIXF,ONLYMS,FINF); COVAL(IN1,FUEL,ONLYMS,FINF)
    *** Outlet:
PATCH(OUT,EAST,%NREGX,%NREGX,#1,#NREGY,1,NZ,1,LSTEP)
 COVAL(OUT,P1,1000.,0.0)
    *** IMMERSOL-walls at the domain boundaries:
PATCH(IMSWL1,SWALL,#1,#NREGX,$1,$1,1,NZ,1,LSTEP)
 COVAL(IMSWL1,H1,WALLCO,GRND4); COVAL(IMSWL1,T3,GRND4,TWALL)
 COVAL(IMSWL1,U1,WALLCO,  0.0)
PATCH(IMSWL2,NWALL,#1,#NREGX,%NREGY,%NREGY,1,NZ,1,LSTEP)
 COVAL(IMSWL2,H1,WALLCO,GRND4); COVAL(IMSWL2,T3,GRND4,TWALL)
 COVAL(IMSWL2,U1,WALLCO,  0.0)
    *** Set emissivity of wall surfaces:
SPEDAT(SET,EMISSIVITY,OF IMSWL1,R,:EMIW1:)
SPEDAT(SET,EMISSIVITY,OF IMSWL2,R,:EMIW2:)

IF(LTURB) THEN
+COVAL(IN1,KE,ONLYMS,FIINIT(KE)); COVAL(IN1,EP,ONLYMS,FIINIT(EP))
+COVAL(OUT,KE,ONLYMS, SAME); COVAL(OUT,EP,ONLYMS, SAME)
+COVAL(IMSWL1,KE,WALLCO,LOGLAW); COVAL(IMSWL1,EP,WALLCO,LOGLAW)
+COVAL(IMSWL2,KE,WALLCO,LOGLAW); COVAL(IMSWL2,EP,WALLCO,LOGLAW)
    *** Eddy-breakup reaction-rate source:
+ PATCH(CHSOTB,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(CHSOTB,FUEL,GRND9,GRND9)
+ CHSOA= FSTOI;  CHSOB= 1.0
ELSE
    *** Power low temperature dependent reaction-rate source:
+ PATCH(CHSOLM,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
+ COVAL(CHSOLM,FUEL,GRND5,0.0)
+ CHSOA= 1.0;  CHSOB= -1./FSTOI;  CHSOC= 1./FSTOI-1.
+ CHSOE= 1.0;  CHSOB= 1.0
ENDIF
    GROUP 15. Termination of sweeps
LSWEEP= 1000;  TSTSWP= -1
    GROUP 16. Termination of iterations
SELREF=T; RESFAC= 0.0001
    GROUP 17. Under-relaxation devices
        DTF=XULAST/UIN/NX
RELAX(P1,LINRLX,0.1)
RELAX(T3,LINRLX,0.05)
RELAX(TMP1,LINRLX,0.05)
    GROUP 18. Limits on variables or increments to them
VARMIN(P1)=-PRESS0+10
    GROUP 22. Spot-value print-out
IXMON= NX/2+1;  IYMON= NY/2+1;  IZMON= 1
LIBREF=211
VARMIN(T3)=TAIR;  VARMAX(T3)=TWALL*3.0
VARMIN(TMP1)=TAIR;VARMAX(TMP1)=TWALL*3.0