TEXT(CHEMKIN -1DZ Plug Flow H2-O2 Reactor
DISPLAY
The example considered is the steady combustion of an H2-O2-N2
mixture in a 1d plug-flow reactor (PSR) at atmospheric pressure.
This case is the same a case 201, except that the calculation is
performed in the z-direction rather than the y-direction. In
addition, the option exists for using the parabolic, rather than
the elliptic solution procedure of PHOENICS. The same results
will be obtained irrespective of the method of solution.
ENDDISPLAY
The results obtained in the 1st cell agree adequately the
following values taken from the SANDIA PSR report:
exit mol fractions
O = 1.54E-2 O2 = 3.48E-2 H = 6.75E-2 H2 = 5.80E-2
OH = 1.34E-2 HO2 = 9.27E-6 H2O = 0.192 H2O2 = 2.23E-5
N2 = 0.619
temperature = 1700 K (fixed)
The run is set up so that the calculation may be performed either
with the CHEMKIN PBP TWOPNT solver (the default), or alternatively
with the PHOENICS solver.
GROUP 1. Run title and other preliminaries
TEXT(CHEMKIN -1DZ Plug Flow H2-O2 Reactor
TITLE
REAL(DTV,DTC,DTT,FLORAT,REAVOL,YH2IN,YO2IN,ZLEN)
BOOLEAN(TWOPNT);CHAR(ANS2)
MESG( Do you want to use CHEMKIN TWOPNT solver? (y/n)
READVDU(ANS,CHAR,Y)
IF(:ANS:.EQ.Y) THEN
+ TWOPNT=T
+ MESG(CHEMKIN PBP SOLVER activated
ELSE
+ TWOPNT=F
+ MESG(PHOENICS SOLVER activated
ENDIF
GROUP 2. Transience; time-step specification
STEADY=T
GROUP 3. X-direction grid specification
REAVOL=67.4;GRDPWR(X,1,REAVOL,1.0)
GROUP 5. Z-direction grid specification
NZ=2;ZLEN=2.0;GRDPWR(Z,NZ,ZLEN,1.0)
GROUP 7. Variables stored, solved & named
MESG( Do you want a PARABOLIC solution? (y/n)
READVDU(ANS2,CHAR,Y)
IF(:ANS2:.EQ.Y) THEN
+ PARAB=T
+ MESG(PARABOLIC solution activated
+ IF(TWOPNT)THEN
+ TEXT(Plug Flow Reactor, Parabolic PBP solver
+ ELSE
+ TEXT(Plug Flow Reactor, Para. PHOENICS solver
+ ENDIF
ELSE
+ PARAB=F
+ MESG(ELLIPTIC solution activated
+ IF(TWOPNT)THEN
+ TEXT(Plug Flow Reactor, Elliptic PBP solver
+ ELSE
+ TEXT(Plug Flow Reactor, Elli. PHOENICS solver
+ ENDIF
ENDIF
SOLVE(P1,W1);STORE(DEN1); OUTPUT(DEN1,Y,Y,Y,Y,Y,Y)
** conduction terms removed from energy equation
SOLVE(TEM1);TERMS(TEM1,N,P,N,P,Y,N)
** Diffusion terms removed from mass-fraction equations
INTEGER(KK,JJ);KK=8
CHEMKIN(SPECIES,H2,H,O2,O,OH,HO2,H2O,N2)
DO II=1, KK-1
+ JJ=CHSOB+II-1
+ TERMS(:JJ:,N,P,N,P,Y,N)
ENDDO
** store: elemental mass fractions of H and N; the heat release
rate per unit volume; the specific enthalpy; the mean
effective specific heat of the mixture; the production
rate of O2; and the rate of reaction 8, i.e HO2 + H
= OH + OH.; the later two are available only when
TWOPNT=T, i.e. when the CHEMKIN TWOPNT solver is
activated.
STORE(ELH,ELN,QDOT,ENTH,SPH1,O2+,8&)
GROUP 9. Properties of the medium (or media)
CP1=GRND9; RHO1=CHEMIST
** Sets the reference pressure (CHSOC in Atmospheres) and CSG4 so
as to identify the CHEMKIN link file: ho11ckln; and the
transport-properties link file: ho11mcln.
** The link files ho11ckln & ho11mcln were created from the
mechanism link file HO11.CKM, which has the form:
CHEMKIN(PROP,ho11,1.0)
GROUP 11. Initialization of variable or porosity fields
DO II=1,KK
+ JJ=CHSOB+II-1
+ FIINIT(:JJ:)=0.
ENDDO
YH2IN=0.025701;YO2IN=0.205607; FIINIT(H2)=YH2IN; FIINIT(O2)=YO2IN
REAL(SUMY,MDOT);SUMY=0.0
DO II=1,KK
+ JJ=CHSOB+II-1
+ SUMY=SUMY+FIINIT(:JJ:)
ENDDO
FIINIT(N2)=1.-SUMY; FIINIT(TEM1)=1700.
** initialise the density field
PATCH(ICHEMK1,INIVAL,1,NX,1,NY,1,NZ,1,1)
INIT(ICHEMK1,DEN1,0.0,GRND1)
GROUP 13. Boundary conditions and special sources
**
FLORAT=364.0;MDOT=FLORAT/REAVOL
INLET(NOCPCK1,LOW,1,NX,1,NY,1,1,1,LSTEP)
VALUE(NOCPCK1,P1,MDOT); VALUE(NOCPCK1,W1,GRND9)
VALUE(NOCPCK1,TEM1,GRND9);TMP1A=298.
DO II=1,KK
+ JJ=CHSOB+II-1
+ VALUE(NOCPCK1,:JJ:,FIINIT(:JJ:))
ENDDO
*** For the PARABolic solution no outlet is required
IF (.NOT.PARAB) THEN
+ OUTLET(OUT1,HIGH,1,NX,1,NY,NZ,NZ,1,LSTEP); VALUE(OUT1,P1,0.0)
ENDIF
IF(TWOPNT) THEN
** If CHSOA=GRND9, ie. the implicit TWOPNT PBP algorithm is used
for the energy and species equations
+ CHEMKIN(REACT,TWOPNT,TEM1)
ELSE
+ CHEMKIN(REACT,PHOENICS,TEM1)
ENDIF
GROUP 15. Termination of sweeps
IF (.NOT.PARAB) THEN
+ IF(TWOPNT) THEN
+ LSWEEP=40;NPLT=5
+ ELSE
+ LSWEEP=800;NPLT=50
+ ENDIF
ENDIF
GROUP 16. Termination of iterations
DO II=1,KK
+ JJ=CHSOB+II-1
+ ENDIT(:JJ:)=1.E-6
ENDDO
ENDIT(TEM1)=1.E-6
IF (PARAB) THEN
+ IF(TWOPNT) THEN
+ LITHYD=40;NPRMON=20;NPLT=5
+ ELSE
+ LITHYD=600;NPRMON=300;NPLT=50
+ ENDIF
ENDIF
GROUP 17. Under-relaxation devices
DTV=3.E-5;DENPCO=T
NB: SELREF=F in CHEMKIN interface if TWOPNT solver activated.
IF(TWOPNT) THEN
+ DTT=3.E-5;DTC=3.E-4
ELSE
** set FIXVAL=1.E15 to prevent solution cut-out of TEM1
+ FIXVAL=1.E15;DTC=2.*DTV;DTT=0.007*DTV
ENDIF
CHEMKIN(RELAX,DTC)
RELAX(TEM1,FALSDT,DTT); RELAX(W1,FALSDT,DTV)
RESREF(TEM1)=-1.E-3; VARMIN(TEM1)=200.
GROUP 19. Data communicated by satellite to GROUND
** The link files ho11ckln, ho11mcln & ho11.ckm reside in the
directory: d_earth\d_opt\d_chem. ho11.ckm has the form:
ELEMENTS O H N END
SPECIES H2 H O2 O OH HO2 H2O N2 END
REACTIONS JOULES/MOLE
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
END
**
GROUP 22. Spot-value print-out
TSTSWP=-1;ITABL=2
GROUP 24. Dumps for restarts
INIFLD=T;IZMON=1;NZPRIN=1
distil=t
if(parab) then
+ if(twopnt) then
EX(P1 )=1.015E+05;EX(W1 )=3.190E+04;EX(H2 )=3.445E-03
EX(H )=1.481E-03;EX(O2 )=3.162E-02;EX(O )=5.401E-03
EX(OH )=7.878E-03;EX(HO2 )=1.054E-06;EX(H2O )=1.815E-01
EX(N2 )=7.687E-01;EX(8& )=1.903E-04;EX(O2+ )=1.449E-01
EX(SPH1)=1.207E+07;EX(ENTH)=2.338E+02;EX(QDOT)=3.379E+10
EX(ELN )=7.687E-01;EX(ELH )=2.570E-02;EX(TEM1)=1.695E+03
EX(DEN1)=1.693E-04
+ else
EX(P1 )=1.015E+05;EX(W1 )=3.190E+04;EX(H2 )=3.445E-03
EX(H )=1.481E-03;EX(O2 )=3.162E-02;EX(O )=5.401E-03
EX(OH )=7.878E-03;EX(HO2 )=1.054E-06;EX(H2O )=1.815E-01
EX(N2 )=7.687E-01;EX(8& )=1.000E-10;EX(O2+ )=1.000E-10
EX(SPH1)=1.208E+07;EX(ENTH)=8.644E+03;EX(QDOT)=3.359E+10
EX(ELN )=7.687E-01;EX(ELH )=2.570E-02;EX(TEM1)=1.695E+03
EX(DEN1)=1.693E-04
+ endif
else
+ if(twopnt) then
EX(P1 )=5.076E+04; EX(W1 )=2.495E+04; EX(H2 )=4.469E-03
EX(H )=2.084E-03; EX(O2 )=4.504E-02; EX(O )=5.945E-03
EX(OH )=5.643E-03; EX(HO2 )=1.829E-06; EX(H2O )=1.681E-01
EX(N2 )=7.687E-01; EX(8& )=6.828E-04; EX(O2+ )=4.698E-01
EX(SPH1)=1.156E+07; EX(ENTH)=2.662E+02; EX(QDOT)=5.431E+10
EX(ELN )=7.687E-01; EX(ELH )=2.570E-02; EX(TEM1)=1.489E+03
EX(DEN1)=1.929E-04
+ else
EX(P1 )=9.147E+04;EX(W1 )=4.002E+04;EX(H2 )=4.803E-03
EX(H )=1.955E-03;EX(O2 )=3.507E-02;EX(O )=9.756E-03
EX(OH )=1.722E-02;EX(HO2 )=1.438E-06;EX(H2O )=1.578E-01
EX(N2 )=7.687E-01;EX(8& )=1.000E-10;EX(O2+ )=1.000E-10
EX(SPH1)=1.310E+07;EX(ENTH)=1.718E+06;EX(QDOT)=4.652E+10
EX(ENTH)=1.718E+06;EX(ELH )=2.570E-02;EX(TEM1)=2.188E+03
EX(DEN1)=1.287E-04;EX(ELN )=7.734E-01
+ endif
endif