CHEMKIN - Introduction

The CHEMKIN System
The PHOENICS CHEMKIN Interface
Organisation and Data
Specification of Composition
Energy Conservation Equation

1 Introduction

The CHEMKIN System

The CHEMKIN system is supplied by Sandia National Laboratories and is described by Kee et al (1993a). CHEMKIN consists of:

A thermodynamics database
A library of Fortran subroutines which the user may call from his applications programmes to supply thermodynamic data, reaction rate data and various other thermo-chemical data
A "stand-alone" interpreter programme that reads a "plain language" file that specifies the chemical elements, the chemical species, the reactions, and any additional thermodynamic data for the thermo-chemical system under investigation

Associated with CHEMKIN, and also supplied by Sandia National Labs, is a further system that supplies transport data; this is described by Kee et al (1993b). The transport properties system consists of:-

A transport database.
A library of Fortran subroutines which the user may call from his application programmes to supply viscosities, thermal conductivities, diffusion coefficients, and thermal diffusion ratios or coefficients calculated according to two approximations.
A fitting programme that generates polynomial fits to the detailed transport properties in order to make more efficient the calculations performed by the subroutine library.

(b)The PHOENICS CHEMKIN Interface

The PHOENICS CHEMKIN Interface provides a range of facilities, from which the user may choose those that he requires. The facilities provided are:

calculation of thermodynamic and transport properties of gas mixtures;
an improved formulation of the transport equations in the "mixture-averaged" approximation, including thermal-diffusion (the Soret effect);
calculation of source terms for laminar flow for chemical species and enthalpy that result from chemical reactions when the system is to be solved within the usual PHOENICS methodology;
a fully implicit solver for the solution of the chemical species and enthalpy variables in laminar-flow reacting systems with stiff chemical source terms;
the calculation of inlet properties, ie. the inflowing enthalpy of the gas-mixture, and either the mass-flux from the inlet gas velocity, or an inlet gas speed from the inlet mass-flux;

The thermodynamic properties and the reaction rates are obtained by making calls to routines in the CHEMKIN Library, and the transport properties are obtained by making calls to the Sandia National Lab.s Multi-component Transport Library (TRANLIB). In addition the two-point boundary problem solver TWOPNT may be used to solve the difference equations when the reaction rates are known to be large.

1.1 Organisation and Data

The interface to the CHEMKIN and TRANLIB routines takes the form of a single large subroutine GXCHKI which is called from a number of points in GREX3, from GXPRPS (for physical properties), and from GXRHO and GXENUL.

The subroutine is not entirely straightforward because the CHEMKIN and TRANLIB routines calculate properties, or reaction rates, for all species in a cell simultaneously, whereas PHOENICS EARTH requires the properties for a single species for all cells in a slab at each call.

The differences in the order of data access are resolved by making the calls to the CHEMKIN and TRANLIB routines when the first species is accessed; the properties are then stored in dedicated F-array segments, created using GXMAKE calls, which have NX*NY*K elements. Note that K is the number of species specified in the CHEMKIN input file.

The primary data input to PHOENICS is, as usual, through the Q1 file and it is the Q1 file that controls the options in the CHEMKIN interface. However, if the CHEMKIN SATLIT programme CHEMST is run, the SATELLITE will read the CKLINK file and make settings on the basis of its contents.

However, in addition to the EARDAT file generated by the SATELLITE from the Q1 file, the user must supply the following CHEMKIN file:

    CKLINK - the CHEMKIN link file generated by the user from
             
the mechanism file xxxx.CKM by running the
             
CKINTERP program

and if the Sandia Transport Library is used to supply the transport properties;

    TPLINK - the Transport Properties link file generated by
             
the user from the CKLINK file by running the
             
TRANFIT program

The data-flow may be visualised thus:


                 --------                     -------
  ---------     /        /      --------     /       /
  |xxxx.CKM|-->| CKINTERP |---->|CKLINK|--->| TRANFIT |
  ---------     /        /      --------     /       /
                 --------          |          -------
                                  / /            |
                                /    /       --------
                              /       /      |TPLINK|
                            /          /     --------
                          /             /        |
                        /                /       |
                      /                   /      |
                     V                     V     V
              ---------                    -------
  ------     /         /                  /       /
  | Q1 |--->| SATELLITE |--------------->|  EARTH  |
  ------     /         /                  /       /
              ---------                    -------

If variable CSG4 is set in the Q1 file, then the non-blank characters of CSG4 are used to construct a file-name for the CKLINK and TPLINK files. For example, if the setting

    CSG4 = 'ho11'

is made, the CKLINK data will be read from the file

    ho11ckln

and the TPLINK data will be read from the file

    ho11mcln

The PHOENICS-EARTH code will attempt to read the CKLINK and TPLINK files from the users private directory, however, if one of the files is not found, then an attempt will be made read the file from the directory;

    PHOENICS/d_earth/d_spefor/d_chmkin

After a successful read will report the files used in the RESULT file and in the screen output. The user may change the directories searched by modifying the appropriate lines of the EARCON file. The same procedure is followed by the SATELLITE when the CHEMST coding is activated (see below for further details).

1.2 Basics of the Model

a) Specification of Composition

There is a choice to be made regarding the variables used to specify the composition of the gas mixtures that will be used. Compositions may be specified in terms of mass fractions, Y, mole-fractions, X, and concentrations, C. From the point of view of the conservation equations to be solved the most convenient variables to use are the mass fractions, Y.

For laminar cases, the last species of the set is not solved, but instead is derived from the mass-fractions of the other species and the knowledge that the sum of all mass fractions must be 1.

The index of the first CHEMKIN species in the PHOENICS dependent variable list defaults to 16, but may be set by the user to be greater than 16, by setting the variable

    CHSOB = 'index of 1st CHEMKIN species'

b) Energy Conservation Equation

Choices must also be made regarding the form of the energy conservation equation. The first choice is whether to solve for temperature or for an enthalpy. If enthalpy is solved, then the temperature must be recovered by means of an iterative procedure from the enthalpy and composition, whereas if temperature is solved the enthalpy may be calculated directly from the temperature and composition. In addition, the specification of conductive heat fluxes through the PHOENICS COVAL mechanism is much simpler when temperature solution has been chosen. So, we solve for temperature (TEM1). A further choice must be made regarding the formulation of the energy conservation equation.

Essentially, the equation is written as a conservation equation for enthalpy and so it must be decided whether to include the bulk kinetic energy, and whether to include the chemical energy in the enthalpy; the exclusion of either necessitates a source term for each contribution excluded. It has been decided to exclude the kinetic energy from the enthalpy, and to solve for a reduced enthalpy, h', defined

       h'(T) = h(T) - h(To)

where To is a convenient reference temperature usually chosen to be 273K.

In PHOENICS when the temperature, TEM1, is solved the energy conservation equation takes, for one direction, the form

   
  M.(<Cp>.T - <Cp>".T")/dt
   
+ m(+).(<Cp>(+).T(+) - <Cp>.T)
   
+ m(-).(<Cp>(-).T(-) - <Cp>.T)
   
- (A(+).k(+)/dX(+))(T-T(+)) - (A(-).k(-)/dX(-))(T-T(-))
   
- So(h) = 0

wherein +/- refer to the + and - direction neighbour cells, " refers to values at the previous time-step, M is the mass of fluid in the cell, m is a mass-flow rate which is either positive or zero, A is a cell-face area, dX is the distance between cell centres, k is the conductivity, So is the energy source term, and <Cp> is an effective specific heat capacity. A, dX and k are defined on the cell-faces, whilst all other quantities are defined at cell-centres.

The form of the effective specific heat-capacity is dictated by the need to recover the correct energy conservation equation for which we require

       h'(T) = <Cp>.T

ie.

       <Cp> = h'(T)/T

The source term in the energy conservation equation arising from chemical reactions is easily derived, and takes the form:


                    K
       So(h) = -M. Sum(w(k).W(k).ho(k))
                   k=1

wherein w(k) is the molar production rate of species k, W(k) is the molecular-mass, and ho(k) is the zero-point enthalpy defined as:

       ho(k) = h(To;k)

so that

       h'(T;k) = h(T;k) - ho(k)