1 PHOENICS INTERFACE TO CHEMKIN - CONTENTS (i)
PANELS 1 - 2 Contents 3 - 4 The CHEMKIN System 5 - 7 The PHOENICS CHEMKIN Interface 8 Data Input for CHEMKIN and TRANLIB 9 The Mechanism File 10 An Example of a Mechanism File 11 Input for PHOENICS (General Features) 12 Energy Equation 13 Transport Model 14 - 16 Physical Properties 17 - 18 Reaction Sources 19 - 20 Boundary Conditions - Inflow 21 - 22 Example of Inflowing Boundary Conditions 23 Wall Conditions 24 - 29 Further Problem Specification Facilities 30 Additional Print-out rev 1.1DA
2 PHOENICS INTERFACE TO CHEMKIN - CONTENTS (ii)
PANEL 31 General Advice on Convergence Control 32 Summary and Concluding Remarks 33 Sources of Further Information
3 THE CHEMKIN SYSTEM
The CHEMKIN system is supplied by Sandia National Laboratories.
It consists of:
* a thermodynamics database;
* a library of Fortran subroutines for the retrieval of thermodynamic data, reaction-rate data, and various other thermo-chemical data; and
* a 'stand-alone' interpreter program which reads a user-input file containing all the relevant information about the thermo-chemical system under investigation ( Such as chemical elements, chemical species, the reactions etc.).
4 THE CHEMKIN SYSTEM II
* A molecular-transport properties system which comprises:
a transport database;
a library of Fortran subroutines enabling retrieval of viscosities, thermal conductivities, diffusion coefficients, and thermal diffusion ratios or coefficients; and
a curve-fitting program which provides polynomial fits to transport properties data.
5 The PHOENICS CHEMKIN Interface
The PHOENICS CHEMKIN Interface facilitates:
* calculation of thermodynamic and transport properties of gas mixtures;
* improvements to the standard PHOENICS mass-transport model, including the Curtiss-Hirschfelder model and thermal-diffusion (the Soret effect) for species diffusion;
* calculation of the source terms resulting from chemical reactions in laminar flow for chemical species and enthalpy;
* a fully implicit solver for the solution of the chemical species and enthalpy variables in tightly-coupled (stiff) reacting sytems;
6 The PHOENICS CHEMKIN Interface II
* the calculation of inlet properties, such as the inflowing enthalpy of the gas-mixture, the mass-flux from an inlet gas velocity, or an inlet gas speed from an inlet mass-flux; and
* the use of the CHEMKIN two-point boundary-problem solver, TWOPNT, for systems in which reaction source terms dominate transport terms.
The thermodynamic properties and reactions rates are obtained by calling routines in the CHEMKIN Library.
The transport properties are obtained by calling routines in the TRANLIB library.
The subroutine GXCHKI provides the interface to CHEMKIN and TRANLIB.
Subroutine GXCHKI is called from a number of points in GREX3, from GXPRPS, GXRHO and GXENUL.
7 The PHOENICS CHEMKIN Interface III
For turbulent flows, species and heat diffusion is normally dominated by the eddy viscosity and so molecular diffusion may be represented by a mixture viscosity and appropriate Prandtl/Schmit numbers.
The CHEMKIN reaction rates may be used in turbulent flow, but interactions between turbulence and chemical reaction are not yet included in the CHEMKIN interface. Consequently, its use in turbulent flows is subject to large uncertainties.
8 Data input for CHEMKIN and TRANLIB
The user must supply the following file for CHEMKIN and TRANLIB
CKLINK CHEMKIN link file generated from a mechanism file ###.CKM by running CKINTERP program
TPLINK Transport Properties link file generated from the CKLINK file by running TRANFIT program
Hence the EARTH run requires CKLINK and TPLINK files in addition to the EARDAT file. The PIL variable CSG4 in the Q1 file identifies the name for the LINK files which would be the same as the name given to the mechanism file with ckln and mcln extension. For example if HOL.CKM is the mechanism file then CSG4='HOL'.
A script file CKI.BAT ( for PC ) is supplied with the PHOENICS CHEMKIN installation which will create the link files based on the mechanism file supplied by the user.
9 The Mechanism File
The mechanism file has the extension *.CKM and its contents provide a symbolic description of the chemical reaction.
This would include information on elements, species, thermodynamic data, and the reaction mechanism.
The order in which the CHEMKIN interpreter reads in the description is: element data; species data; thermodynamic data (optional), and reactions.
The rules governing the specification of the required data are provided in the SANDIA Report SAND89-8009b [1993]
10 An Example of a Mechanism File
The mechanism file listed below provides element, species and reaction data for single step octane-oxygen reaction
ELEMENTS H N O C END
SPECIES C8H18 O2 CO2 H2O N2 END
REACTIONS JOULES/MOLE C8H18 + 25O2 => 8CO2 + 9H2O 7.803E9 0 167.467 END
11 Input for PHOENICS (General Features)
The problem definition for PHOENICS is as usual through the Q1 file.
At present the units employed are restricted to those of CHEMKIN ie cgs.
Hence the geometry should be specified in cm, the mass in g and time in s.
The CHEMKIN pressure is in units of dynes/cm^2 and enthalpy is in units of ergs/mole.
The composition may be specified in terms of mass fractions, mole fractions or concentrations. In PHOENICS the most convenient choice is mass fractions.
The spare C# equations will be activated for the CHEMKIN species. The index for the first species in the PHOENICS dependent variable list will default to 16, but may be set to be greater than 16 by
CHSOB = Index of the first CHEMKIN species.
12 Energy Equation
The temperature, TEM1, is selected as the dependant variable for the energy equation. This equation is written as a conservation equation for the reduced enthalpy, h(T) - h(To), where To is a reference temperature.
For further details of the energy equation refer to POLIS entry 'CHEMKIN, for chemical kinetics' under 'about PHOENICS 2.1',' Description of PHOENICS Modules', 'The add-ons'.
13 Transport Model
The default transport model built into PHOENICS may be used with the thermodynamic and transport properties calculated in the CHEMKIN and TRANLIB routines.
An alternative to the default model is the 'Curtiss-Hirschfelder' model of mass transport and the inclusion of thermal diffusion (the Soret effect).
The enhanced transport model is activated by setting two PIL variables,
ENULA = GRND9 - enhanced mixture-averaged transport model
ENULB = GRND9 - thermal diffusion (Soret effect)
14 Physical Properties I
A call is made to the CHEMKIN routine CKRHOY to calculate the density.
This is activated by setting,
RHO1 = GRND9
in the Q1 file or if PRPS is used by putting GRND9 in the density specification in the PROPS data.
An option exists for calculating the effective specific heat capacities from the species enthalpies calculated by calling the CHEMKIN routine CKHMS.
This option is activated by setting
CP1 = GRND9
in the Q1 file or GRND9 into the PROPS data for Cp.
15 Physical Properties II
An option exists for calculating the mixture-averaged dynamic viscosity by calling the TRANLIB routine MCAVIS. The kinematic viscosity is then calculated by division by density.
This option is activated by setting
ENUL = GRND9 in the Q1 file or
placing GRND9 for ENUL in the PROP data.
The mixture thermal conductivity can be calculated by calling the TRANLIB routine MCACON. This is activated by setting
PRNDTL(TEM1)=-GRND9 in the Q1 file or
by placing GRND9 for thermal conductivity in PROPS data.
16 Physical Properties III
Mixture-averaged diffusivities for all species can be calculated by calling TRANLIB routine MCADIF. This option is activated by setting, for the species variable numbered K,
PRNDTL(K) = -GRND9 in the Q1 file.
Activation through PRPS is not possible.
If the enhanced transport model is activated, then the computed diffusivity is divided by the mean molecular weight.
17 Reaction sources (non-stiff reactions)
For non-stiff or weakly-stiff reactions, the reaction source terms are activated as follows:
PATCH(CHEMK1,VOLUME,1,NX,1,NY,1,NZ,1,LSTEP) DO K=1,KS + COVAL(CHEMK1,C:K:,GRND9,GRND9) ENDDO COVAL(CHEMK1,TEM1,GRND9,GRND9)
where KS is the number of species solved, which is one minus the total number of species involved.
18 Reaction sources (stiff reactions)
When the reaction source terms dominate the terms in the difference equations for advection and diffusion, and when the time-scales for reactions are widely varying, then an implicit method is required which solves all species simultaneously.
The approach employed is to solve all species and enthalpy equations implicitly in a point-by-point manner, using the CHEMKIN's TWOPNT solver.
The implicit solver is activated by setting
CHSOA = GRND9
The mechanisms available for controlling the TWOPNT solver are discussed later on under 'Further Problem Specification II'.
19 Boundary Conditions - Inflow (energy)
It is necessary to specify the properties of gas mixtures flowing into the domain in a way which is consistent with the gas properties to be used within the domain.
For the energy equation the incoming enthalpy must be specified. Hence, when temperature (TEM1) is being solved for, The PATCH name must commence with NOCPCK.
This indicates that the VALUE will not be multiplied by the in-cell specific heat, and that the VALUE is an enthalpy set by the CHEMKIN interface. The gas composition at inlet is obtained by the COVALs for the mass-fractions, and the temperature of the inlet gas is obtained from the PIL variable TMP1A.
20 Boundary Conditions - Inflow (mass-inflow)
Two options for specifying the mass-flux are provided.
(1) Specify the mass-flux from inlet data on flow speed and density calculated from the inflowing gas mixture, or
(2) Specify the inflowing mass-flux and deduce the inlet flow speed from this and the density calculated from the inflowing gas mixture.
A limited facility exists whereby Poiseuille velocity and mass-flux profiles may be calculated for W1 velocity for LOW and HIGH inlets for non-BFC cases.
21 Example of Inflowing Boundary Conditions
An example of the PIL commands for setting an inflow condition based on prescribed inlet velocities is:
PATCH(NOCPCK1,SOUTH,1,NX,1,1,1,NZ,1,LSTEP) COVAL(NOCPCK1,P1, FIXFLU,GRND9) for calculating the flux COVAL(NOCPCK1,U1,ONLYMS,30.) COVAL(NOCPCK1,V1,ONLYMS,40.) gas speed=sqrt(u1**2+v1**2+w1**2) COVAL(NOCPCK1,H2,ONLYMS,0.02) COVAL(NOCPCK1,O2,ONLYMS,0.25) COVAL(NOCPCK1,TEM1,ONLYMS,GRND9) for calculating enthalpy TMP1A = 298.0
22 Example of Inflowing Boundary Conditions II
An example of the PIL commands for setting an inflow condition based on prescribed mass influx and finite normal velocity only, is:
PATCH(NOCPCK1,SOUTH,1,NX,1,1,1,NZ,1,LSTEP) COVAL(NOCPCK1,P1, FIXFLU,1.E-4) COVAL(NOCPCK1,V1,ONLYMS,GRND9) COVAL(NOCPCK1,H2,ONLYMS,0.02) COVAL(NOCPCK1,O2,ONLYMS,0.25) COVAL(NOCPCK1,TEM1,ONLYMS,GRND9) TMP1A = 298.0
23 Boundary Condition - Wall conditions
Discussion here is limited to the required wall conditions for mass fractions. Laminar wall functions should be activated as appropriate.
When the enhanced diffusion model is in use, the laminar wall functions should be setup in the Q1 file, as follows,
PATCH(CKWALIN,SWALL,1,NX,1,1,1,NZ,1,LSTEP) COVAL(CKWALIN,H2,GRND9,0.025)
The first five characters of the PATCH name should be CKWAL.
24 Further Problem Specification Facilities I (PRESS0 & TEMP0)
The reference pressure and temperature, PRESS0 and TEMP0, have both been utilised by the CHEMKIN interface, and all the stored pressures and temperatures are relative to the prevailing values of PRESS0 and TEMP0.
Additionally the value of PRESS0 is set to be a multiple of the standard atmospheric pressure and is obtained from CHEMKIN data-base. The value of the pressure in atmospheres is set via PIL variable CHSOC.
The default value of pressure passed to CHEMKIN and TRANLIB is PRESSO, unless VARPRS T is set in the Q1 file between the line CHKIBEGIN and CHKIEND. If it is, then the local instantaneous pressure will be used.
25 Further Problem Specification Facilities II (Control of TWOPNT)
The action of the TWOPNT solver may be controlled by making settings in the Q1 file which are read directly by EARTH. The settings must be preceeded by a line CHKIBEGIN and must be followed by a line CHKIEND. Starting in column 3 or over, the following settings may be made:
CALL1 - if 1 then for problems for which the energy (TEM1) equation is solved a preliminary solution with fixed temperature should be performed, if 2 or 0 then only the full solution should be performed
NJAC - maximum no. of iterations of N-R solver before the Jacobian is updated without internal false- timestepping (Default: 20) Experience shows that a suitable value may be 50.
26 Further Problem Specification Facilities III (Control of TWOPNT)
ITJAC - maximum no. of iterations of N-R solver before the Jacobian is updated with false-timestepping (Default: 50)
IRETIR - no. of timesteps to be taken before the false- timestep is increased (Default: 50)
NUMDT - the number of time-steps to be taken in the event of a failure of the Newton iteration (Default: 50)
NINIT - no. of false-timesteps taken by the N-R solver prior to attempting a solution without internal false- timestepping (Default: 0)
27 Further Problem Specification Facilities IV (Control of TWOPNT)
UFAC - factor by which the time-steps are increased when IRETIR steps have been performed (Default: 2)
DFAC - factor by which the time-steps are decreased when a time-step fails to converge (Default: 4.1)
DTMIN - the lower limit below which the false-timstep is not permitted to fall (Default: 1.E-10)
DTMAX - the upper limit above which the false-timestep is not permitted to rise (Default: 0.1)
28 Further Problem Specification Facilities IV (Control of TWOPNT)
ATOL - the absolute tolerance used by TWOPNT to determine whether the steady-state problem has converged according to the formula
dVAR < (ATOL + RTOL.VAR)
(Default: 1.E-4)
RTOL - the relative tolerance used by TWOPNT to determine whether the steady-state problem has converged according to the formula given above (Default: 1.E-8)
ATIM - the absolute tolerance used by TWOPNT to determine whether the false-transient problem has converged according to the formula given above (Default: 1.E-4)
29 Further Problem Specification Facilities VI (Control of TWOPNT)
RTIM - the relative tolerance used by TWOPNT to determine whether the false-transient problem has converged according to the formula given above (Default: 1.E-8)
30 Additional Print-out
Facilities exits for storing heat-release rate, rate of production of concentrations, rate of reaction of selected reactions and elemental composition of mass. Mole-fractions are computed and printed at completion. The storage is activated as:
STORE(QDOT) for storage of the heat-release rate
STORE(Cn+) for the rate of production of concentration Cn
STORE(n&) for the rate of reaction of reaction n, eg STORE(4&) for rate of reaction 4
STORE(ELXX) for the mass composition of element XX
31 General Advice on Convergence Control
It is recommended that inertial relaxation should be applied to all species being solved for and the energy equation.
VARMIN and VARMAX should be set to 0.0 and 1.0, respectivly for all species.
ENDIT should be set to at least 1.E-6 for all species and temperature.
RESREF's for species should be of the order of 1.E-6. Suitable RESREF's should be set for mass, momentum and energy. For the TWOPNT solver, RESREF's for species and temperature must be negative, and at present SELREF should be set False.
For large density gradients, DENPCO = T, may assist convergence.
32 Summary and Concluding Remarks
The lecture has presented the facilities provided by the PHOENICS CHEMKIN interface.
Guidance on problem setup, and files required has been provided.
Although the CHEMKIN interface is primaryly for laminar flow, the CHEMKIN reaction rates may be used in turbulent flow if the interaction between turbulence and chemical reaction can be tolerated.
The CHEMKIN interface is supplied with a number of Q1 files for exemplary cases.
FURTHER DEVELOPMENTS include: provision of more library Q1's for practical applications; improved documentation and ease of use; possible allowance for turbulence-chemistry interactions; extension of automatic inlet conditions to BFCs.
33 Sources of Further Information
SANDIA REPORTS by R J Kee, F M Rupley, J A Miller, August 1993
The Twopnt Program for Boundary Value Problems, SANDIA REPORT SAND91-8230 by J F Grear, April 1992
CHEMKIN Interface entry in POLIS under About PHOENICS-2.1 -> Descriptions of the PHOENICS modules -> The add-ons -> CHEMKIN.
Simple Library examples using CHEMKIN in the Chemical Reaction Flow Library: cases C201-C211 inclusive.