TALK=T;RUN( 1, 1) ** LOAD(x113) from the x Input Library PHOTON USE p; phi 15 1 set con fi dep 3 con fuel x 1 fi;0.01 gr ou x 1 te 1 fuel mass-fraction contours 0.13265E+04 0.21215E+04 CR msg Pressto continue pause cl con prod x 1 fi;0.01 gr ou x 1 te 1 product contours 0.13265E+04 0.21215E+04 CR msg Press to continue pause cl con mixf x 1 fi;0.01 gr ou x 1 te 1 mixture fraction contours 0.13265E+04 0.21215E+04 CR msg Press to continue pause cl con tmp1 x 1 fi;0.01 gr ou x 1 te 1 temperature contours 0.13265E+04 0.21215E+04 CR msg Press to continue pause cl vec x 1 sh gr ou x 1 te 1 velocity vectors 0.13265E+04 0.21215E+04 CR msg Press to continue pause end enduse TEXT(2-D CH4/AIR combustion+thermal NOX: C113 TITLE mesg(PC486/50 time last reported as 1.min DISPLAY The problem considered is two-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 9 jzdbck 18 jsdbck 10 gn2ox 0.767 gn2pr 0.725 gco2pr 0.151 gh2opr 0.124 crkbug f dbgnox f noxdend Variable declarations BOOLEAN(NOXCAL) -- Change NOXCAL to T to perform a NOx post-processing run CHAR(NOXC) MESG( MESG(Main combustion calculation (M) or NOX Post-processing (P) MESG( (Default = M) READVDU(NOXC,CHAR,M) IF(:NOXC:.EQ.P)THEN TEXT(NOX Post-processing calculation : C113 +NOXCAL=T ELSE +NOXCAL=F ENDIF 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=4.E5 NS=7;NR=3 REAL(FLOWFU,FLOWOX,EQUIV,RINFU,AREAFU,AREAOX,PI) PI=3.141592 ** EQUIV is the equivalence ratio EQUIV=1.1;WINFU = 0.6 TFUEL = 289.0;TOX = 810.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 3. X-direction grid specification CARTES=F;XULAST=0.1 GROUP 4. Y-direction grid specification REAL(ROUTER);ROUTER=0.01;RINFU=0.5*ROUTER INTEGER(NYFU,NYOX);NYFU=12;NYOX=12 NREGY=2 IREGY=1;GRDPWR(Y,NYFU,RINFU,-1.5) IREGY=2;GRDPWR(Y,NYOX,ROUTER-RINFU,1.5) GROUP 5. Z-direction grid specification NZ=35;ZWLAST=60.*ROUTER GRDPWR(Z,NZ,ZWLAST,1.0) GROUP 7. Variables stored, solved & named IF(NOXCAL) THEN + STORE(P1,V1,W1,KE,EP,MIXF,H1) ELSE + SOLVE(P1,V1,W1,MIXF,H1) + SOLUTN(P1,Y,Y,Y,P,P,P) + TURMOD(KEMODL) ENDIF STORE(TMP1,RHO1,OXID,PROD,FUEL,ENUT) IF(NOXCAL) THEN + SOLVE(XNO,XN) ELSE + STORE(XNO,XN) ENDIF STORE(XH,XO,XOH,XO2,CRKT,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 CP1=GRND10 CP1A=CPFU;CP1B=CPOX;CP1C=CPPR GROUP 10. Inter-phase-transfer processes and properties GROUP 11. Initialization of variable or porosity fields ** fuel flow rate for the solution domain AREAFU=0.5*XULAST*RINFU*RINFU FLOWFU=RHOFU*WINFU*AREAFU ** stoichometric air supply FLOWOX=STOIC*FLOWFU ** actual air supply FLOWOX=EQUIV*FLOWOX AREAOX=0.5*XULAST*ROUTER*ROUTER-AREAFU WINOX=FLOWOX/(AREAOX*RHOOX) FIINIT(W1) = WINFU; FIINIT(H1) = HFUEL REAL(KEINF,EPINF,KEINOX,EPINOX) KEINF=(0.05*WINFU)**2 EPINF=0.1643*KEINF**1.5/(0.1*RINFU) KEINOX=(0.05*WINOX)**2 EPINOX=0.1643*KEINOX**1.5/(0.1*(ROUTER-RINFU)) FIINIT(KE)=KEINF;FIINIT(EP)=EPINF 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/2+1,1,1,1,1) + COVAL(FUELNOX,XNO,RHOFU*WINFU,0.0) + COVAL(FUELNOX,XN,RHOFU*WINFU,0.0) + PATCH(AIRNOX,LOW,1,NX,NY/2+1,NY,1,1,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/2,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 ) + COVAL(INLETF,KE,ONLYMS,KEINF) + COVAL(INLETF,EP,ONLYMS,EPINF) * air inlet + PATCH(INLETA,LOW,1,NX,NY/2+1,NY,1,1,1,1) + COVAL(INLETA,P1 ,FIXFLU, RHOOX*WINOX) + COVAL(INLETA,W1 ,ONLYMS, WINOX ) + COVAL(INLETA,H1 ,ONLYMS, HOX ) + COVAL(INLETA,MIXF,ONLYMS, 0.0 ) + COVAL(INLETA,KE,ONLYMS,KEINOX) + COVAL(INLETA,EP,ONLYMS,EPINOX) * outlet boundary + PATCH(OUTLET,HIGH,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) * wall boundary + WALL(WFNN,NORTH,1,NX,NY,NY,1,NZ,1,1) ENDIF GROUP 15. Termination of sweeps FSWEEP=1;LSWEEP=250 REAL(MDOT);MDOT=FLOWFU+FLOWOX RESREF(P1)=MDOT*1.E-12 RESREF(W1)=(FLOWFU*WINFU+FLOWOX*WINOX)*1.E-12 RESREF(V1)=RESREF(W1) RESREF(H1)=(FLOWFU*HFUEL+FLOWOX*HOX)*1.E-12 RESREF(KE)=RESREF(P1)*KEINF;RESREF(EP)=RESREF(P1)*EPINF RESREF(MIXF)=FLOWFU*1.E-12 RESREF(XNO)=RESREF(P1);RESREF(XN)=RESREF(P1) LITER(P1 )=25 GROUP 17. Under-relaxation devices REAL(DTF);DTF=ZWLAST/WINOX/NZ RELAX(P1 ,LINRLX,1.0) RELAX(V1,FALSDT,0.1*DTF);RELAX(W1,FALSDT,DTF) RELAX(KE,FALSDT,DTF);RELAX(EP,FALSDT,DTF) RELAX(H1,FALSDT,10.*DTF);RELAX(MIXF,FALSDT,10.*DTF) RELAX(RHO1,LINRLX,0.1);RELAX(XNO,FALSDT,DTF*0.25) GROUP 18. Limits on variables or increments to them VARMIN(TMP1)=0.1*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)=50.0;VARMIN(RHO1)=0.001 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 OUTPUT(TMP1,Y,N,Y,Y,Y,Y) GROUP 22. Spot-value print-out IXMON=1;IYMON=9;IZMON=18 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(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=150 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 TSTSWP=-1;DENPCO=T DISTIL=T IF(NOXCAL) THEN + EX(P1 )=5.829E+01;EX(V1 )=3.029E-02;EX(W1 )=1.095E+01 + EX(KE )=7.682E-01;EX(EP )=2.113E+02;EX(H1 )=6.619E+06 + EX(NOSR)=6.377E+03;EX(CRKT)=1.967E+03;EX(FUEL)=7.205E-02 + EX(PROD)=7.303E-01;EX(OXID)=1.977E-01;EX(RHO1)=7.681E-01 + EX(TMP1)=2.059E+03;EX(MIXF)=1.121E-01;EX(XO2 )=2.199E-03 + EX(XO )=2.951E-05;EX(XH )=2.821E-06;EX(XN )=2.021E-09 + EX(XNO )=1.487E-04;EX(XOH )=7.033E-04 ELSE + EX(P1 )=5.829E+01;EX(V1 )=3.029E-02;EX(W1 )=1.095E+01 + EX(KE )=7.682E-01;EX(EP )=2.113E+02;EX(H1 )=6.619E+06 + EX(FUEL)=7.205E-02;EX(PROD)=7.303E-01;EX(OXID)=1.977E-01 + EX(RHO1)=7.279E-01;EX(TMP1)=2.059E+03;EX(MIXF)=1.121E-01 + EX(ENUT)=4.299E-04;EX(NOSR)=1.000E-10 ENDIF LIBREF = 113 STOP