PHOTON USE
  AUTOPLOT 
  file
  phi 5
 
  cl
  da 1
  temp
  scale
  colf 1
  msg temp of gas - red
  msg Press  to continue
  pause
  cl
  da 1
  w1
  col3 1
  msg gas velocity - blue
  msg Press  to continue
  pause
  cl
  da 1
  fuel
  da 1
  mixf
  scale
  col3 1
  colf 2
  msg volume fraction of methane - blue
  msg mixture fraction - red
  msg Press  to end
  pause
  end
  enduse
 
TITLE
mesg(PC486/50 Time Last Reported As 1.min
  DISPLAY
   The problem considered is one dimensional and concerns
   the calculation of thermal nox emissions. A two-step
   approach is used whereby the first step involves the main
   exothermic reaction of methane in air, and the second
   involves the solution of the Zeldovich reaction scheme for
   the formation of the major nitric oxide: nitrogen oxide, NO.
   The first PHOENICS calculation considers the methane
   combustion only, so that the logical variable NOXCAL=F. The
   second and final PHOENICS calculation sets NOXCAL=T which
   deactivates the solution of all dependent variables
   excepting those involved in the nox chemistry.
  ENDDIS
 
   In the first PHOENICS calculation a single-step SCRS
   scheme is taken for methane combustion, which is presumed
   to be diffusion controlled. Thus, the thermochemistry of
   the hydrocarbon oxidation is described by the solution of
   equations for the mixture fraction and the mixture enthalpy,
   wherein the definition of the latter includes the chemical
   potential of the reactants.
 
   In the second PHOENICS calculation, the species the product
   species  O2, O, OH, H, H2, CO & CO2 are calculated in such
   proportions as are appropriate to equilibrium stoichometric
   combustion at the prevailing pressure and enthalpy. This
   calculation is performed from a minimization of the Gibbs
   free energy via the CREK program of Pratt & Wormeck (1976).
   The calculation then computes the mass fractions of NO & N
   via the Zeldovich reaction mechanism by solving species
   conservation equations which have as convection and diffusion
   fluxes the ones determined from the first PHOENICS calculation.
    *
  noxdbegin
    ns       7
    nr       3
    jxdbck   1
    jydbck   1
    jzdbck   1
    jsdbck   30
 
    gn2ox    0.767
    gn2pr    0.725
    gco2pr   0.151
    gh2opr   0.124
 
    crkbug   t
    dbgnox   t
  noxdend
    *
    Variable declarations
BOOLEAN(NOXCAL);NOXCAL=F
REAL(CMOL,HMOL,NMOL,OMOL,O2MOL,N2MOL,CH4MOL,H2OMOL,COMOL,CO2MOL)
REAL(NOMOL,OHMOL,OXMOL,FUMOL,PRMOL)
REAL(CO2PRF,H2OPRF,N2OXF,N2PRF,O2OXF,CH4MFU)
    Molecular masses
CMOL = 12.01115;HMOL=1.00797;NMOL=14.0067
OMOL = 15.9994;CH4MOL=CMOL+4.0*HMOL; O2MOL = 2.0*OMOL
N2MOL = 2.0*NMOL;H2OMOL=2.0*HMOL+OMOL
COMOL = CMOL+OMOL;CO2MOL= 2.0*OMOL+CMOL
NOMOL = NMOL+OMOL; OHMOL = OMOL+HMOL
    Mass fractions
O2OXF = 0.233;N2OXF = 0.767;CO2PRF= 0.151
H2OPRF= 0.124;N2PRF = 0.725;CH4MFU= 1.0
OXMOL = O2OXF*O2MOL+N2OXF*N2MOL; FUMOL = CH4MOL
PRMOL = CO2PRF*CO2MOL + H2OPRF*H2OMOL + N2PRF*N2MOL
    Additional constants & initialization
REAL(GASCON,TFUEL,TADIAB,TOX,CPOX,CPFU,CPPR,STOIC,FSTOI)
REAL(SETMIN,HFU,HOX,HFUEL,WINOX,WINFU,RHOOX,RHOFU)
INTEGER(NS,NR)
GASCON=8314.0;PRESS0=1.E5
NS=7;NR=3
REAL(FLOWFU,FLOWOX,EQUIV); EQUIV=1.1
WINFU = 10.0;TFUEL = 289.0;TOX   = 800.0
HFU = 4.9E7; STOIC=17.24; FSTOI=1./(1.+STOIC)
CPOX = 1500.0 ; CPPR = 1500.0 ; CPFU = 1500.0
HOX  = CPOX*TOX;HFUEL= CPFU*TFUEL + HFU
RHOOX = PRESS0*OXMOL/(GASCON*TOX)
RHOFU = PRESS0*FUMOL/(GASCON*TFUEL)
TADIAB=(HFU+HFUEL+STOIC*HOX)/(CPPR*(1.+STOIC))
SETMIN=1.E-10
    GROUP 4. Y-direction grid specification
NY=1;YVLAST=1.0
    GROUP 5. Z-direction grid specification
NZ=34;ZWLAST=34.;ZFRAC(1)=-34 ;ZFRAC(2)=1./34.
    GROUP 7. Variables stored, solved & named
IF(NOXCAL) THEN
+ STORE(P1,W1,MIXF,H1)
ELSE
+ SOLVE(P1,W1,MIXF,H1)
+ SOLUTN(P1,Y,Y,Y,P,P,P)
ENDIF
STORE(TMP1,RHO1,OXID,PROD,FUEL)
IF(NOXCAL) THEN
+ SOLVE(XNO,XN)
ELSE
+ STORE(XNO,XN)
ENDIF
STORE(XH,XO,XOH,XO2,CRKT,PRDT,NOSR)
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,Y,Y,N,Y,N)
    GROUP 9. Properties of the medium (or media)
RHO1 = 3GASES ; RHO1A=FUMOL ; RHO1B=OXMOL ; RHO1C=PRMOL
TMP1 = SCRSEQ ; TMP1A=CPFU  ; TMP1B=CPOX  ; TMP1C=CPPR
tmp2a= FSTOI ; tmp2b=HFU
    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
FLOWFU=RHOFU*WINFU
  ** Stoichmetric air supply
FLOWOX=STOIC*FLOWFU
  ** Actual air supply
FLOWOX=EQUIV*FLOWOX;WINOX=FLOWOX/(RHOOX*ZWLAST)
FIINIT(W1) = WINFU; FIINIT(H1) = HFUEL
    GROUP 12. Unused
    GROUP 13. Boundary conditions and special sources
IF(NOXCAL) THEN
+ PATCH(NOXSR,FREEVL,1,NX,1,NY,1,NZ,1,1)
+ COVAL(NOXSR,XNO,FIXFLU,GRND)
+ COVAL(NOXSR,XN,GRND,GRND)
+ PATCH(FUELNOX,LOW,1,NX,1,NY,1,1,1,1)
+ COVAL(FUELNOX,XNO,RHOFU*WINFU,0.0)
+ COVAL(FUELNOX,XN,RHOFU*WINFU,0.0)
+ PATCH(AIRNOX,NORTH,NX,NX,NY,NY,1,NZ,1,1)
+ COVAL(AIRNOX,XNO,RHOOX*WINOX,0.0)
+ COVAL(AIRNOX,XN,RHOOX*WINOX,0.0)
ELSE
    * fuel inlet
+ PATCH(INLETF,LOW,1,NX,1,NY,1,1,1,1)
+ COVAL(INLETF,P1  ,FIXFLU, RHOFU*WINFU)
+ COVAL(INLETF,W1  ,ONLYMS, WINFU  )
+ COVAL(INLETF,H1  ,ONLYMS, HFUEL  )
+ COVAL(INLETF,MIXF,ONLYMS, 1.0    )
    * air inlet
+ PATCH(INLETA,NORTH,NX,NX,NY,NY,1,NZ,1,1)
+ COVAL(INLETA,P1  ,FIXFLU, RHOOX*WINOX)
+ COVAL(INLETA,W1  ,ONLYMS, WINOX )
+ COVAL(INLETA,H1  ,ONLYMS, HOX   )
+ COVAL(INLETA,MIXF,ONLYMS, 0.0   )
    * outlet boundary
+ PATCH(OUTLET,CELL,1,NX,1,NY,NZ,NZ,1,1)
+ COVAL(OUTLET,P1  , FIXP , 0.0)
+ COVAL(OUTLET,H1  ,ONLYMS,SAME)
+ COVAL(OUTLET,MIXF,ONLYMS,SAME)
+ COVAL(OUTLET,FUEL,ONLYMS,SAME)
ENDIF
    GROUP 15. Termination of sweeps
FSWEEP=1;LSWEEP=50
REAL(MDOT);MDOT=RHOFU*WINFU
RESREF(P1)=MDOT*1.E-12
RESREF(W1)=WINFU*RESREF(P1)
RESREF(H1)=HFUEL*RESREF(P1); RESREF(MIXF)=RESREF(P1)
RESREF(XNO)=RESREF(P1); RESREF(XN)=RESREF(P1)
LITER(P1)=25
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
REAL(RELAXN); RELAXN=ZWLAST/WINFU/NZ
RELAX(P1,LINRLX,1.0); RELAX(W1,FALSDT,RELAXN)
RELAX(H1,FALSDT,RELAXN*1.E5)
RELAX(MIXF,FALSDT,RELAXN*1.E5)
RELAX(RHO1,LINRLX,0.1)
RELAX(XNO,FALSDT,RELAXN*0.25)
    GROUP 18. Limits on variables or increments to them
VARMIN(TMP1)=TFUEL; VARMAX(TMP1)=4000
VARMAX(OXID)=1.0; VARMIN(OXID)=0.0
VARMAX(PROD)=1.0; VARMIN(PROD)=0.0
VARMAX(MIXF)=1.0; VARMIN(MIXF)=0.0
VARMAX(FUEL)=1.0; VARMIN(FUEL)=0.0
VARMAX(RHO1)=10.0; VARMIN(RHO1)=0.01
VARMAX(XNO )=1.0; VARMIN(XNO)=0.0
VARMAX(XN)=1.0; VARMIN(XN)=0.0
    GROUP 19. Data communicated by satellite to GROUND
LSG63=NOXCAL
  NAMGRD=NOX;USEGRD=F;USEGRX=T
  RG(1)=CO2PRF; RG(2)=H2OPRF; RG(3)=N2PRF
  RG(4)=O2OXF; RG(5)=N2OXF; RG(6)=CH4MOL
  RG(7)=CO2MOL; RG(8)=HMOL; RG(9)=H2OMOL
  RG(10)=NMOL; RG(11)=NOMOL; RG(12)=N2MOL
  RG(13)=OMOL; RG(14)=OHMOL; RG(15)=O2MOL
  RG(16)=PRESS0; RG(17)=GASCON; RG(18) = TOX
  RG(19)=TFUEL; RG(20)= HOX; RG(21) = HFUEL
  RG(22)=STOIC; RG(23)=TADIAB; RG(24) = SETMIN
  IG(1)=NS;IG(2)=NR
    * lg(1) no longer used
  LG(1)=F
    * temperature selection for nox calc.
  BOOLEAN(CRKTMP,PRDTMP,FLDTMP)
  CRKTMP=T;PRDTMP=F; FLDTMP=F
  LG(2)=CRKTMP;LG(3)=PRDTMP;LG(4)=FLDTMP
    *..............NOX CALCULATION SWITCH(This is obsolete.  sas)
       LG(5)=NOXCAL
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=1;IZMON=12;NPRMON=1
    GROUP 23. Field print-out and plot control
ITABL=3;IPLTL=LSWEEP
IPLTF=FSWEEP;YZPR=T;NPRINT=LSWEEP
NZPRIN=1;NUMCLS=5;NPLT=5
IF(NOXCAL) THEN
+ OUTPUT(XNO,Y,Y,Y,Y,Y,Y); OUTPUT(XN,Y,Y,Y,Y,Y,Y)
ELSE
+ OUTPUT(XH,N,N,N,N,N,N); OUTPUT(XO,N,N,N,N,N,N)
+ OUTPUT(XOH,N,N,N,N,N,N); OUTPUT(XO2,N,N,N,N,N,N)
+ OUTPUT(CRKT,N,N,N,N,N,N); OUTPUT(PRDT,N,N,N,N,N,N)
+ OUTPUT(XNO,N,N,N,N,N,N); OUTPUT(XN,N,N,N,N,N,N)
+ OUTPUT(NOSR,N,N,N,N,N,N)
ENDIF
    GROUP 24. Dumps for restarts
IF(NOXCAL) THEN
+ RESTRT(ALL);LSWEEP=30
ENDIF
FIINIT(XNO)=0.0; FIINIT(XN)=0.0
FIINIT(XH)=0.0; FIINIT(XO2)=0.0
FIINIT(XO)=0.0; FIINIT(XOH)=0.0
FIINIT(CRKT)=0.0; FIINIT(PRDT)=0.0
TSTSWP=-1;DENPCO=T
  BOOLEAN(DBGNOX)
  INTEGER(JXDBCK,JYDBCK,JZDBCK,JSDBCK)
  ** NOX & equilibrium debug switch
  DBGNOX=T;JXDBCK=1;JYDBCK=1;JZDBK=12
  JSDBCK=1
  LG(6)=DBGNOX
  IG(8)=JXDBCK;IG(9)=JYDBCK;IG(10)=JZDBCK
  IG(11)=JSDBCK
  #nk minimum value for pressure& NEW DEF. OF CP1
cp1=grnd10
cp1a=CPFU ;cp1b=CPOX;cp1c=CPPR
VARMIN(P1))=-0.5*PRESS0