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