A Tutorial on Chemkin
Niket Kaisare (11-12-2007)
Indian Institute of Technology Madras Original Presentation: Dr. Ashish Mhadeshwar (June 2004)
Why use Chemkin?
1. Efficient handling of large reaction mechanisms 2. Generalized framework, simple and standard inputs 3. Minimization of common mistakes 4. Common platform for distributing new products
Example: Ammonia Decomposition on Ruthenium
NH3 NH3* NH2* NH* H*
H2
Ruthenium
Requirements
We intend to solve a combined flow-diffusion-reaction problem 1. 2. 3. 4. Reaction stoichiometry and rate constants Reaction thermo-chemistry Thermodynamic properties Transport properties (diffusion, viscosity and thermal conductivity)
Example: Ammonia Decomposition on Ruthenium Matrix of stoichiometric coefficients
R1 1. 2. 3. 4. 5. 6. N2 + 2* H2 +2* NH* + * NH2* + * NH3* + * NH3 + * 2N* 2H* N* + H* NH* + H* NH2* + H* NH3* N* H* NH* NH2* NH3* 2 2 R2 R3 1 1 -1 1 1 -1 1 -1 1 1 R4 R5 R6
Kinetic parameters for each reaction
T r j = A0 j T ref
j E j m e RT akj
Reaction thermo-chemistry H , G, S
Thermodynamic and Transport properties
cp , H , ; k g , Dik ,
Example: Ammonia Decomposition on Ruthenium
Matrix of stoichiometric coefficients
R1 R2 R3 R4 R5 R6 N* 2 1 1. N2 + 2* 2N* 2. H2 +2* 2H* H* 2 1 1 1 [Link] NH* NH* + * N* + H* 1. -1 1 4. NH2* + * NH* + H* -1 1 5. NH3* + * NH2* + H* NH2* [Link] NH3 + * NH3* NH3* -1 1
straightforward
T j E j Kinetic parameters for 3. Adding/removingreactionsfromkj e RT a m r j = A0 j each reaction T ref mechanismiscumbersome
Reaction thermo-chemistry H , G, S Thermodynamic and Transport properties cp , H , ; k g , Dik ,
Utilities
AURORA/SENKIN AURORA predicts the timeevolution or steady state of a well mixed reactor.
CRESLAF simulates laminar, chemically reacting, boundary-layer flow in cylindrical or planar channels.
Utilities
PLUG simulates a plug-flow reactor with gas-phase and surface chemistry.
EQUIL calculates the equilibrium state of systems containing ideal gas mixtures or ideal solutions.
Utilities
OPPDIF predicts the structure of opposed-flow diffusion flames in cylindrical or planar coordinate systems.
PREMIX predicts the structure of steady, laminar, 1D, pre-mixed flames.
Utilities
SHOCK: A program for predicting chemical behavior behind incident and reflected shock waves.
SPIN: A program for modeling one-dimensional, rotating-disk or stagnation-flow chemicalvapor-deposition reactors.
Components of Chemkin Suite
1. Chemkin 2. Surface Chemkin 3. Transport
Components of Chemkin Suite
Chemistry Input Properties Input Link File
ChemKin Utilities
Interpreter
Library Application
CHEMKIN Interpreter Library Chemistry Input Physical Properties Input Link File cklink ckinterp.x cklib.f fort.15
SURFACE TRANSPORT skinterp.x sklib.f [Link] tranfit.x tranlib.f --trandat tplink
thermdat sklink
Note: Surface Chemkin and Transport Libraries require cklink for species information
Example: Ammonia Decomposition on Ruthenium
NH3 NH3* NH2* NH* H*
H2
Ruthenium
1. CHEMKIN: To define gas phase species and thermodynamics 2. SURFACE CHEMKIN: To define surface species and reactions 3. TRANSPORT: To define gas phase transport properties (diffusion, viscosity and thermal conductivity) if applicable
Defining Gas Phase Conditions in Chemkin
fort.15 thermdat
ChemKin Utilities
ckinterp.x
cklink
cklib.f Application
Steps: Create fort.15 (gas phase reaction mechanism)
Input of elements, species, and reactions in fort.15
Input of thermodynamic database: thermdat Compile ckinterp.f to create ckinterp.x Run ckinterp.x to create fort.16 and cklink (MUST READ fort.16) Compile cklib.f to create cklib.o and link with your code
1. Chemkin inputs: Sample fort.15
ELEMENTS N H AR RH END SPECIES N2 NH3 END REACTIONS END H2 AR
5 elements
6 species
0 reactions
Chemkin peculiarity: Catalyst should be declared with gas phase species
2. Chemkin inputs: Thermodynamic database thermdat
Low
THERMO 300.000 1000.000 5000.000 NH3 L 7/88N 1H 3
Mid High
G 200.000 3500.000 1000.000
1 2
3.85746029E+00 4.41437026E-03-2.21481404E-06 5.23490188E-10-4.72084164E-14
-4.87591660E+04 2.27163806E+00 2.35677352E+00 8.98459677E-03-7.12356269E-06 3 2.45919022E-09-1.43699548E-13-4.83719697E+04 9.90105222E+00 4
Chemkin peculiarity: Fixed format that must be followed for every new entry
3. Check output file fort.16 for errors Run ckinterp.x. Two output files are created: cklink and fort.16
Defining Surface Conditions in Chemkin
[Link] cklink thermdat
Surface ChemKin
skinterp.x
sklink
sklib.f Application
Steps: Create [Link] (surface reaction mechanism) Input surface species and reactions Input of thermodynamic database: thermdat Compile skinterp.f to create skinterp.x
Run skinterp.x to create [Link] and sklink (MUST READ [Link]) Compile cklib.f to create cklib.o and link with your code
1. Surface Chemkin inputs: Sample [Link]
SITE/SURFACE/ H(S) NH2(S) END N(S) NH3(S) SDEN/1.66058e-9/ NH(S) RU(S)
site density (mol/cm2s)
surface species units reaction
BULK RU(B)/12.37/ END
REACTIONS H2 + 2RU(S) STICK
kcal/mole => 2H(S) + 2RU(B) 1.0 0.0 1.9
2H(S) + 2RU(B) => H2 + 2RU(S)
6.022e21
0.0
23.7
N2 + 2RU(S) STICK
=> 2N(S) + 2RU(B)
1.0
0.0
6.2
sticking keyword
6.022e21 0.0 50.3 => N(S) + H(S) + RU(B) 6.022e19 0.0 5.8
2N(S) + 2RU(B) => N2 + 2RU(S)
NH(S) + RU(S)
COV/N(S) 0.0 0.0 15.4/ N(S) + H(S) + RU(B) => NH(S) + RU(S) 6.022e19 0.0 37.3
Coverage dependent parameters
T E k = A0 exp Tref RT
COV/N(S) 0.0 0.0 -19.7/
...
1. Surface Chemkin inputs: Peculiarities in [Link] Site density is a 3 element vector (gas, surface, bulk). Surface site density, i.e., sden(2) is typically used in the calculations, other two are zero. For every consumed vacant surface site, a bulk site is produced. However, for reactions in the category of A(S) + B(S) C(S) + D(S), no vacant surface sites or bulk sites are involved. Reactions are written in irreversible manner, due to unavailability of thermodynamic database (specific heats) for surface species. Units of pre-exponential factors are dependent on the type of reaction in the mechanism. The pre-exponential factor should be defined such that the reaction rate is calculated in mol/cm2sec.
2. Surface Chemkin inputs: Thermodynamic database thermdat
RU(S) 81292RU 1 S 300.00 2000.00 1000.00 1 2 0.02318546E+02 0.01929426E-01-0.01682130E-05 0.04131402E-09-0.16692935E-14
-0.07634172E+04-0.09802473E+02 0.15056368E+01 0.09217769E-01-0.01955590E-03 3 0.02020330E-06-0.07251403E-10-0.07226661E+04-0.06689818E+02 N(S) TPIS79N 1 S 300.000 3000.000 1000.000 4 1 2 3 4
0.00000E+00 0.00000E-00 -0.00000E-00 -0.00000E-00 0.00000E-00 -0.00000E+00 -0.00000E+00 -0.00000E+00 0.00000E-00 0.00000E-00 -0.00000E-00 0.00000E-00 -0.00000E+00 0.00000E+00
Follows same format as earlier Can be provided in [Link] or thermdat Thermodynamic inputs for surface species are usually not available. Put all specific heats as zero, since all heats of reactions are calculated using BOC Thermodynamic consistency is therefore important
Defining Surface Conditions in Chemkin
trandat cklink
Transport Utilities
tranfit.x
tplink
tranlib.f Application
Steps: Input of transport database: thermdat Compile tranfit.f and dmath.f to create tranfit.x Run tranfit.x to create tplink Compile tranlib.f to create tranlib.o and link with your code
10
1. Transport inputs: Transport database trandat (optional)
AR C C2 C2O CN2
0 136.500 0 1 71.400 97.530
3.330 3.298 3.621 3.828 3.828
0.000 0.000 0.000 0.000 0.000
0.000 0.000 1.760 0.000 0.000
0.000 0.000 ! * 4.000 1.000 ! * 1.000 ! OIS
1 232.400 1 232.400
species molecular index
L-J collision diameter
polariz ability rotational relaxation collision number
L-J potential well depth
dipole moment
Interfacing with Application codes Typical structure of Chemkin subroutines:
subroutine ckcpbs(t, yk, ickwrk, rckwrk, cpmean)
Name Input(s)
Work arrays
Output(s)
Interfacing with Application codes: Chemkin Mnemonics ck = Chemkin sk = Surface Chemkin mc = Transport iwk = integer work array rwk = real work array cwk = character work array
11
Interfacing with Application codes: Dimensionalize Work Arrays
! *** File identifiers integer, parameter :: linc = 25, linksk = 26, linkmc = 35, lout = 6
! *** Dimensionalizing work arrays integer, parameter :: leniwk = 6000, lenrwk = 6000, lencwk = 6000 integer, parameter :: lsiwk = 6000, lsrwk = 6000, lscwk = 6000 integer, parameter :: lenimc = 6000, lenrmc = 6000 integer :: ickwrk(leniwrk), iskwrk(lsiwk), imcwrk(lenimc) double precision :: rckwrk(lenrwrk), rskwrk(lsrwk) rmcwrk(lenrmc) character (len=16) :: cckwrk(lencwk), cskwrk(lscwk)
Interfacing with Application codes: Initialization
! *** Open chemkin link files open (linc, file='cklink', form='unformatted') open (linkmc, file='tplink', form='unformatted') open (linksk, file='sklink', form='unformatted') ! *** Initialize chemkin call ckinit (leniwk, lenrwk, lencwk, linc, lout, ickwrk, rckwrk, cckwrk) call mcinit (linkmc, lout, lenimc, lenrmc, imcwrk, rmcwrk) call skinit (lenisk,lenrsk,lencsk,linksk,lout,iskwrk, rskwrk,cskwrk)
Interfacing with Application codes: Using subroutines
!.......tutorial open (unit=394, file='[Link]',status='unknown') write(394,*)'Inlet mole fractions' write(394,39)xin call ckxty(xin,ickwrk,rckwrk,yin) write(394,*)'Inlet mass fractions' write(394,39)yin
Inlet mole fractions 0.300000 0.600000 0.100000 0.000000
Inlet mass fractions 0.743363 0.106195 0.150442 0.000000
12
Interfacing with Application codes: Using subroutines
write(394,*)'Initial coverages' write(394,39)(act(ii),ii=kgmax+1,ksmax-1) call skrat (p,temper,act,sden,iskwrk,rskwrk,sdot,sitdot) write(394,*)'Rates of species production' write(394,39)sdot
pressure activity Gas: mole fractions Surface: coverage Bulk: 1.0
Site density sden(2) is the site density of the catalyst surface
Initial coverages 0.10000E-05 0.99990E+00 0.10000E-05 0.98000E-04
Rates of species production -0.53416E-05 -0.21331E-05 0.45925E-05 0.00000E+00 0.21331E-05 -0.45925E-05 0.24594E-05 -0.24594E-05 0.10678E-10
Recent modifications in Surface Chemkin 1. UBI-QEP implementation using a new subroutine (subroutine ubi) internally called from subroutine skrrop 2. Modification of sklib.f to incorporate Molecular beam effects (subroutines skatcz, skden) Large reaction mechanisms (subroutine skinu) Replacement of temperature exponents using a new subroutine (subroutine skrbex) Modification of Arrhenius temperature exponent formula (subroutine skrrop)
13
Useful links
Reaction Design: [Link] History of chemkin: [Link] Overview of chemkin: [Link] Equilibrium calculations: [Link] GRI-Mech: [Link] Sensitivity, reaction path, and principal component analysis modules: [Link] [Link] [Link]
Useful links
Cantera - MATLAB based Chemkin: [Link] Burcats extended thermodynamic database: [Link] Chemical reaction network toolbox: [Link] EGLIB A multicomponent transport software: [Link] Caltech gas phase mechanism library: [Link] Software aid for developing mechanisms (unit conversion, deleting species, etc.): [Link]
14
Useful links
Information about Equil, Aurora, Senkin, Premix, Spin, Creslaf, Surftherm: [Link] A general link for chemical engineers: [Link]
15
