TEXT(CHEMKIN-1DY Plug Flow H2-O2 Reactor 201
  DISPLAY
  The example considered is the steady combustion of an H2-O2-N2
  mixture in a 1d plug-flow reactor (PSR) at atmospheric pressure.
  The case is a steady-state emulation of the transient zero-d PSR
  test case described in the CHEMKIN SANDIA Report 'PSR: A Fortran
  program for modelling well-stirred reactors', P.Glarborg et al,
  SAND86-8209, (1992). Species diffusion and heat conduction are
  absent. As with all CHEMKIN calculations, cgs units are used.
 
  The first cell is intended to represent the reactor volume of
  67.4cm^3, and the nominal residence time is 0.03msec. The inlet
  mass flow rate is 364g/s and the inlet mole fractions are:
  H2=0.275; O2=0.1375; and N2=0.5875. The inlet temperature is 298K
  and there is no heat loss from the reactor. The H2-O2 reaction
  mechanism, which is defined in the file HO11.CKM, is simpler than
  that used in the SANDIA report. The reaction scheme used is taken
  from; "Numerical Methods in Laminar Flame Propagation", Vol. 6,
  Notes on Numerical Fluid Mechanics, Ed. Peters, N. & Warnatz, J.
  (1992).
  ENDDIS
  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
MESG(
TEXT(CHEMKIN-1DY Plug Flow H2-O2 Reactor 201
REAL(DTF,DTC,FLORAT,REAVOL,YH2IN,YO2IN,YLEN)
BOOLEAN(TWOPNT)
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 4. Y-direction grid specification
NY=2;YLEN=2.0;GRDPWR(Y,NY,YLEN,1.0)
    GROUP 7. Variables stored, solved & named
  ** conduction terms removed from energy equation
SOLVE(P1,V1);STORE(DEN1);SOLVE(TEM1);TERMS(TEM1,N,P,N,P,Y,N)
  ** store: elemental mass fractions of H, N & O; the heat release
            rate per unit volume; and the mean effective specific
            heat of the mixture
STORE(ELH,ELN,ELO,QDOT,SPH1)
INTEGER(KK,JJ);KK=8;ARRAY(CIN,REAL,KK)
CHEMKIN(SPECIES,H2,H,O2,O,OH,HO2,H2O,N2)
  ** CHEMKIN(SPECIES,... sets CHSOB=16 (=C1) & LSG61=T
       CHSOB = PHOENICS Variable No of 1st CHEMKIN species
       LSG61 = T  activates calls to GXCHKI from GREX3
  ** Diffusion terms removed from mass-fraction equations
DO II=1,KK-1
+ JJ=II+CHSOB-1
+ TERMS(C:JJ:,N,P,N,P,Y,N)
ENDDO
    GROUP 9. Properties of the medium (or media)
CP1=GRND9; RHO1=CHEMIST
  ** Sets the reference pressure CHSOC (in Atmospheres) and
     CSG4='ho11' which identifies the CHEMKIN link file:
     ho11ckln; and the transport-properties link file: ho11mcln.
CHEMKIN(PROP,ho11,1.0)
    GROUP 11. Initialization of variable or porosity fields
DO II=1,KK
+ JJ=II+CHSOB-1
+ FIINIT(:JJ:)=0.0
ENDDO
YH2IN=0.025701;YO2IN=0.205607; FIINIT(H2)=YH2IN; FIINIT(O2)=YO2IN
REAL(SUMY,MDOT);SUMY=0.0
DO II=1,KK-1
+ JJ=II+CHSOB-1
+ SUMY=SUMY+FIINIT(:JJ:)
ENDDO
FIINIT(N2)=1.-SUMY; FIINIT(TEM1)=1700.
DO II=1,KK
+ JJ=II+CHSOB-1
+ CIN(:II:)=FIINIT(:JJ:)
ENDDO
  ** 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,SOUTH,1,NX,1,1,1,NZ,1,LSTEP)
VALUE(NOCPCK1,P1,MDOT); VALUE(NOCPCK1,V1,GRND9)
VALUE(NOCPCK1,TEM1,GRND9);TMP1A=298.
DO II=1,KK
+ JJ=II+CHSOB-1
+ VALUE(NOCPCK1,:JJ:,CIN(:II:))
ENDDO
OUTLET(OUT1,NORTH,1,NX,NY,NY,NZ,NZ,1,LSTEP); VALUE(OUT1,P1,0.0)
 
IF(TWOPNT) THEN
  ** Activate the implicit TWOPNT PBP algorithm for the energy
     and species equations, i.e. CHSOA=GRND9
+ CHEMKIN(REACT,TWOPNT,TEM1)
ELSE
+ CHEMKIN(REACT,PHOENICS,TEM1)
ENDIF
    GROUP 15. Termination of sweeps
IF(TWOPNT) THEN
+ LSWEEP=40;NPLT=5;TSTSWP=-1
ELSE
+ LSWEEP=1000;NPLT=50;TSTSWP=-50
ENDIF
    GROUP 16. Termination of iterations
DO II=1,KK
+ JJ=II+CHSOB-1
+ ENDIT(:JJ:)=1.E-6
ENDDO
ENDIT(TEM1)=1.E-10
    GROUP 17. Under-relaxation devices
DTF=3.E-5; RELAX(V1,FALSDT,DTF); VARMIN(TEM1)=200.;DENPCO=T
     NB: SELREF=F in CHEMKIN interface if TWOPNT solver activated.
IF(TWOPNT) THEN
+ DTC=10.*DTF
+ RESREF(TEM1)=-1.E-3; RELAX(TEM1,FALSDT,DTF)
ELSE
  ** set FIXVAL=1.E15 to prevent solution cut-out of TEM1
+ FIXVAL=1.E15;DTC=2.*DTF
+ RELAX(TEM1,FALSDT,0.007*DTF)
ENDIF
DO II=1,KK
+ JJ=II+CHSOB-1
+ RELAX(:JJ:,FALSDT,DTC)
ENDDO
    GROUP 19. Data communicated by satellite to GROUND
  ** The link files ho11ckln & ho11mcln were created from the
     mechanism link file HO11.CKM, which 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
  **
  ** The link files ho11ckln, ho11mcln & ho11.ckm reside in the
     directory: d_earth\d_opt\d_chem.
    GROUP 20. Preliminary print-out
ECHO=T
    GROUP 22. Spot-value print-out
ITABL=3; OUTPUT(DEN1,Y,Y,Y,Y,Y,Y)
    GROUP 24. Dumps for restarts
INIFLD=T;IYMON=1;NYPRIN=1