TALK=T;RUN( 1, 1)
** LOAD(x204) from the x Input Library
DISPLAY
1D Laminar Premixed Burner-Stabilised H2-Air Flame
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The case considered is a steady, 1D, laminar, burner-stabilised,
pre-mixed, stoichiometric H2-air flame at atmospheric pressure.
The reactants enter the domain at 298 K with a velocity of 120
cm/s. The flame is stabilised on a porous-plug burner in the
sense that the heat loss necessary to stabilise the flame is
achieved through conductive heat transfer to the inlet boundary.
The present case involves 8 species and 11 reactions, and Ficks
law is used to represent species diffusion. The computation may
be performed using the PHOENICS point-by-point (PBP) solver, or
alternatively the more efficient CHEMKIN TWOPNT solver which
solves the species and TEM1 equations simultaneously in a PBP
manner.
ENDDIS
The case is based on the test case for laminar-flame propagation
given in the following proceedings: "Influence of Transport Models
and Boundary Conditions on Flame Structure", J.Warnatz, pp.87-111,
in "Numerical Methods in Laminar Flame-Propagation", Ed. N.Peters
& J.Warnatz, Vol.6, Notes on Numerical Fluid Mechanics (1982).
The parameters for the case are:
Pressure = 1 Atm; Tin = 298 K ; mass fractions
Y(H2)in = 0.028522; Y(O2)in = 0.226364;
Y(N2)in = 0.745114
and the mass inflow rate is 0.1026 g/s. The calculation is
performed on a domain of 0.6cm length, and the mesh of 40 cells is
concentrated towards the burner. Further mesh refinement is
required in the vicinity of the flame for a numerically accurate
solution.
The species are: H2 H O2 O OH HO2 H2O N2 and the following
11 reactions are considered:
Kj Bj Ej
H2 + OH = H2O + H 2.2E13 0.0 21.5E3
O2 + H = O + OH 2.2E14 0.0 70.3E3
H2 + O = H + OH 1.8E10 1.0 37.2E3
OH + OH = H2O + O 6.3E12 0.0 4.6E3
H + H + M = H2 + M 6.4E17 -1.0 0.0
OH + H + M = H2O + M 5.0E16 0.0 0.0
O2 + H + M = HO2 + M 1.4E15 0.0 7.9E3
HO2 + H = OH + OH 2.5E14 0.0 7.9E3
HO2 + H = H2 + O2 2.5E13 0.0 2.9E3
OH + HO2 = H2O + O2 1.5E13 0.0 0.0
HO2 + O = O2 + OH 6.3E13 0.0 2.9E3
where Kj=Aj*(T**Bj)**exp(-Ej/{R*T}) and the units are g-mole,
cm**3 & Joules.
If several distinct flow rates are used, and the corresponding
heat fluxes to the porous plug are calculated, the adiabatic
flame speed can be deduced by extraoplation to zero heat flux.
For direct computation of the adiabatic flame speed, an unsteady
calculation may be carried out until a steadily propogating
flame develops with time. The adiabatic flame speed is about
184 cm/s for the present case.
AUTOPLOT USE
file
PHI 5
clear
d 1 MH2;d 1 MO2;d 1 MH2O;d 1 TEM1; mult y 0.0001 4
plot 1 3;col3 4;colf 4;blb3 3;scale x 0. 0.2
msg H2, O2, H2O mole-fraction & Temp(/1000) profiles
msg Press to continue
pause
clear
d 1 MH;d 1 MO;d 1 MOH;d 1 MHO2
mult y 10. 2;mult y 10. 3;mult y 1000. 4
col3 1 4;col6 2;col9 3;colf 4;scale x 0. 0.2
msg H, O(*10), OH (*10) & HO2(*1000) mole fraction profiles
msg Press to continue
pause
clear
msg Heat-release rate profile in erg/s/cm**3
d 1 QDOT;scale x 0. 0.2;plot 1;blb1 1
enduse
GROUP 1. Run title and other preliminaries
TEXT( CHEMKIN - 1DY Premixed H2-Air Flame
TITLE
BOOLEAN(TWOPNT,FICK);CHAR(CHSO)
REAL(YLEN,TIN,SUMY,MDOT,YH2IN,YO2IN)
MESG( Enter required method of solution for chemistry
MESG( Default: CHEMKIN Point-by-Point solver
MESG( The alternative is:
MESG( Enter PBP for PHOENICS Point-by-Point solver
READVDU(CHSO,CHAR,CHS)
CASE :CHSO: OF
WHEN PBP,3
+ TWOPNT=F
+ MESG(PHOENICS Point-by-point Solver
WHEN CHS,3
+ TWOPNT=T
+ MESG(CHEMKIN Point-by-Point Solver
ENDCASE
GROUP 4. Y-direction grid specification
YLEN=0.6;GRDPWR(Y,40,YLEN,1.6)
GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,TEM1);SOLUTN(TEM1,P,P,P,P,P,N)
TERMS(TEM1,N,P,P,P,Y,N)
There are 8 species in the chemical scheme, and N2 is
stored while the remainder are solved. CHEMKIN(SPECIES,...)
also sets CHSOB=C1=16 and LSG61=T. CHSOB defines the PHOENICS
variable that corresponds with the first CHEMKIN species.
CHEMKIN(SPECIES,H2,H,O2,O,OH,HO2,H2O,N2)
INTEGER(KK,JJ);KK=8
DO II=1, KK-1
+ JJ=CHSOB+II-1
+ SOLUTN(:JJ:,P,P,Y,P,P,N)
ENDDO
IF(TWOPNT) THEN
** Sets CHSOA=GRND9 & RESREF's to -1.E-6
+ CHEMKIN(REACT,TWOPNT,TEM1)
ELSE
** Sets CHEMK1 GRND9 Patches & covals for TEM1 & species
+ CHEMKIN(REACT,PHOENICS,TEM1)
ENDIF
STORE(VISL,DEN1,QDOT,ENTH,KOND,SPH1)
Store variables for the elemental and mole composition
STORE(ELH,PRPS,ELO,ELN,MH2,MO2,MH2O,MOH,MH,MO,MHO2)
GROUP 8
DIFCUT=0.
GROUP 9. Properties of the medium (or media)
Activate the CHEMKIN physical property calculations
REAL(PROPGR);PROPGR=GRND9
ENUL=PROPGR;CP1=PROPGR; RHO1=PROPGR
** NB: Only Ficks Law should be used at present due to an
error in the Curtis-Hirschfelder species-diffusion
model (mrm 08/02/95)
FICK=T
IF(.NOT.FICK) THEN
* Select the Curtis-Hirschfelder formulation for
species diffusion
+ ENULA=GRND9
* Select thermal-diffusion (Soret Effect) terms in the transport
equations (if using the Curtiss-Hirschfelder formulation )
+ ENULB=GRND9
+ DO II=1,KK
+ JJ=CHSOB+II-1
+ PRNDTL(:JJ:)=-GRND9
+ ENDDO
MESG(WARNING! Curtis-Hirschfelder implementation in error
ENDIF
PRNDTL(TEM1)=-PROPGR
** Sets the stem of the names for the CKLINK and TPLINK files
(CSG4='ho11') and the reference pressure CHSOC in Atmospheres.
CHEMKIN(PROP,ho11,1.0)
CSG10='q1'
FIINIT(PRPS)=71
MATFLG=T;NMAT=1
71 GRND9 GRND9 GRND9 GRND9 0.
0.0
0.0
0.0
0.0
GROUP 11. Initialization of variable or porosity fields
ARRAY(CIN,REAL,KK)
FIINIT(V1)=120.
DO II=1,KK
+ JJ=CHSOB+II-1
+ FIINIT(:JJ:)=1.E-4;CIN(:II:)=0.0
ENDDO
** Set the inlet mass-fractions of H2 & O2 and then calculate
N2 from overall continuity.
SUMY=0.0;YH2IN=0.028522;YO2IN=0.226364
CIN(1)=YH2IN;CIN(2)=0.0;CIN(3)=YO2IN
FIINIT(H)=1.E-3; FIINIT(O)=0.02; FIINIT(OH)=0.015
FIINIT(H2O)=18*(CIN(1)/2-CIN(2)-CIN(3)/32-CIN(4)/16-CIN(5)/17)
DO II=1,KK
+ JJ=CHSOB+II-1
+ SUMY=SUMY+FIINIT(:JJ:)
ENDDO
FIINIT(N2)=1.-SUMY
** Initialise temperature to final flame temperature
FIINIT(TEM1)=2250.
Initiallise the density to a value consistent with the
other gas properties
INIADD=F
PATCH(ICHEMK1,INIVAL,1,NX,1,NY,1,NZ,1,1)
INIT(ICHEMK1,DEN1,0.0,GRND1)
GROUP 13. Boundary conditions and special sources
** Mass inflow boundary condition
Set the inflowing mass-flux to the inflow density*inflow
speed, and use NOCPCK... PATCHes to set up the inlet
conditions where the inflow speed is derived from the mass
flux and density, and the inflowing thermal enthalpy is
derived from the gas compostion and temperature .
INLET(NOCPCK1,SOUTH,1,NX,1,1,1,NZ,1,LSTEP)
VALUE(NOCPCK1,P1,GRND9); VALUE(NOCPCK1,V1,120.)
DO II=1,KK
+ JJ=CHSOB+II-1
+ COVAL(NOCPCK1,:JJ:,ONLYMS,CIN(:II:))
ENDDO
** Define inlet temperature via SPEDAT for NOCPCK1
( earlier versions of CHEMKIN used TMP1A=298, which is still
supported ).
TIN=298.0
SPEDAT(SET,NOCPCK1,TINLET,R,TIN)
COVAL(NOCPCK1,TEM1,ONLYMS,GRND9)
** Allow heat loss to porous plug for flame stabilisation
PATCH(FLUXIN,SWALL,1,NX,1,1,1,NZ,1,LSTEP)
COVAL(FLUXIN,TEM1,1.0,TIN)
** Outflow boundary
OUTLET(OUT1,NORTH,1,NX,NY,NY,1,NZ,1,LSTEP)
VALUE(OUT1,P1,0.0)
GROUP 15. Termination of sweeps
GROUP 16. Termination of iterations
DO II=1,KK
+ JJ=CHSOB+II-1
+ ENDIT(:JJ:)=1.E-6
ENDDO
ENDIT(TEM1)=1.E-6
GROUP 17. Under-relaxation devices
IF(TWOPNT) THEN
+ LSWEEP=250;NPLT=25
+ CHEMKIN(RELAX,6.E-2)
+ RELAX(TEM1,FALSDT,1.E-5); RELAX(V1,FALSDT,6.E-5)
ELSE
+ LSWEEP=3500;NPLT=100; RESFAC=1.E-4
+ CHEMKIN(RELAX,7.E-3)
+ RELAX(TEM1,FALSDT,1.E-6); RELAX(V1,FALSDT,5.E-5)
ENDIF
GROUP 18. Limits on variables or increments to them
** Prevent TEM1 falling below reference temperature of 0 C.
VARMIN(TEM1)=273.
GROUP 19. Data communicated by satellite to GROUND
GROUP 21. Print-out of variables
Only a 1D case so look at all cells; print every 50 sweeps
NYPRIN=1;IYPRF=1;IYPRL=NY;NPRINT=LSWEEP
OUTPUT(DEN1,Y,Y,Y,N,Y,N); OUTPUT(PRPS,N,P,P,P,P,P)
GROUP 22. Spot-value print-out
TSTSWP=-1;IYMON=NY-1;ITABL=3
GROUP 23. Field print-out and plot control
Generate plots of H2, O2 & TEM1 and H, O, OH, H2O and HO2
PATCH(PROFSY,PROFIL,NX/2+1,NX/2+1,1,NY,NZ/2+1,NZ/2+1,1,LSTEP)
COVAL(PROFSY,H2,0.0,0.0);COVAL(PROFSY,O2,0.0,0.0)
PATCH(PROFSYB,PROFIL,NX/2+1,NX/2+1,1,NY,NZ/2+1,NZ/2+1,1,LSTEP)
COVAL(PROFSYB,TEM1,0.0,0.0);COVAL(PROFSYB,N2,0.0,0.0)
PATCH(PROFSYA,PROFIL,NX/2+1,NX/2+1,1,NY,NZ/2+1,NZ/2+1,1,LSTEP)
COVAL(PROFSYA,H,0.0,0.0);COVAL(PROFSYA,O,0.0,0.0)
COVAL(PROFSYA,OH,0.0,0.0);COVAL(PROFSYA,HO2,0.0,0.0)
COVAL(PROFSYA,H2O,0.0,0.0)
GROUP 24. Dumps for restarts
LIBREF = 204
DISTIL=T
IF(TWOPNT) THEN
+ EX(P1 )=1.330E+01;EX(V1 )=6.926E+02;EX(O2 )=3.609E-02
+ EX(H2O )=2.065E-01;EX(N2 )=7.451E-01;EX(MH )=1.020E-02
+ EX(MH2O)=2.747E-01;EX(MO2 )=2.486E-02;EX(MH2 )=5.084E-02
+ EX(ELN )=7.451E-01;EX(ELO )=2.264E-01;EX(ELH )=2.852E-02
+ EX(SPH1)=1.232E+07;EX(KOND)=1.691E+04;EX(ENTH)=1.910E+05
+ EX(QDOT)=1.003E+10;EX(VISL)=4.671E+00;EX(TEM1)=1.975E+03
+ EX(H2 )=4.629E-03;EX(H )=4.477E-04;EX(O )=1.525E-03
+ EX(OH )=5.686E-03;EX(HO2 )=3.197E-06;EX(MHO2)=2.134E-06
+ EX(MO )=2.219E-03;EX(MOH )=7.992E-03;EX(DEN1)=1.822E-04
ELSE
+ EX(P1)=1.330E+01;EX(V1)=6.926E+02;EX(O2)=3.609E-02
+ EX(H2O)=2.065E-01;EX(N2)=7.451E-01;EX(MH)=1.020E-02
+ EX(MH2O)=2.747E-01;EX(MO2)=2.486E-02;EX(MH2)=5.084E-02
+ EX(ELN )=7.451E-01;EX(ELO)=2.264E-01;EX(ELH)=2.852E-02
+ EX(SPH1)=1.232E+07;EX(KOND)=1.691E+04;EX(ENTH)=1.912E+05
+ EX(QDOT)=1.003E+10;EX(VISL)=4.670E+00;EX(TEM1)=1.975E+03
+ EX(H2 )=4.629E-03;EX(H )=4.477E-04;EX(O )=1.525E-03
+ EX(OH )=5.686E-03;EX(HO2 )=3.197E-06;EX(N2 )=7.451E-01
+ EX(MHO2)=2.134E-06;EX(MO )=2.219E-03;EX(MOH )=7.993E-03
+ EX(DEN1)=1.822E-04
ENDIF
STOP