CHEMKIN III 学习笔记
首先简单说明一下，CHEMKIN由美国Sandia国家实验室开发，CHEMKIN II即以下版本为开源版本，CHEMKIN III开始科研使用需要购买版权，或使用旧版本实现相同功能。CHEMKIN III之后的版本开始带有简单的图形界面。这个软件包是一个比较经典的数值燃烧的代码包，主要提供和化学反应相关的翻译器，利用Arrhenius公式计算反应速率的计算器，以及一些经典的算例。此处说几句题外话.
文章共4,679字 · 阅读需要大约16分钟
首先简单说明一下，CHEMKIN由美国Sandia国家实验室开发，CHEMKIN II即以下版本为开源版本，CHEMKIN III开始科研使用需要购买版权，或使用旧版本实现相同功能。CHEMKIN III之后的版本开始带有简单的图形界面。这个软件包是一个比较经典的数值燃烧的代码包，主要提供和化学反应相关的翻译器，利用Arrhenius公式计算反应速率的计算器，以及一些经典的算例。
此处说几句题外话，CHEMKIN III的完成时间大概是九十年代，当时的编程语言和编程模式发展还非常不成熟，程序本身并不具备面向对象的特性，函数名和变量名也由于时代的限制没有完整的体现功能和存储内容，程序包中的各个fortran文件都是极其冗长的工具箱。这导致CHEMKIN III的代码本身其实结构并不方便阅读和维护。如果对燃烧的基本知识并不熟悉，建议首先学习李新亮openCFD-comb模块的使用手册，模块本身其实是对这个代码库的简单实现，更加方便阅读和使用。
CHEMKIN III是有使用手册的，使用手册中主要介绍的是CKLIB1,2两个代码库，最后介绍了CONP.f函数，函数调用前面的代码库中的一系列函数，完成了一个0维燃烧的算例。另外一个实现的功能是Premix.f中的一维预混燃烧，这里计算的是一个定常火焰燃烧模型，但是这个代码的计算过程使用的是牛顿迭代求解非线性系统，计算过程并非现在常见的时间推进过程。
CHEMKIN的基本结构
代码的基本结构如下：
其中chem.inp为化学反应机理的输入文件，therm.dat为一个热学参数的数据库。将他们分别用翻译器ckinterp.f进行读取，生成输出文件chem.out和linking file文件名chem.bin。linking file可以在Gas-Phase Library即文件cklib.f中使用。具体使用方法如下：
这里的unit 5是指fortran代码中read(5,*)中的文件编号。另外需要注意的是，sample.f中如果想要调用cklib.f中的子函数，必须首先调用subroutine ckinit创建化学反应所使用的数组，并读入linking file中的数据。首先我们给出chem.inp的一个例子并进行说明：
ELEMENTS
H O AR
END
SPECIES
H2 O2 H O OH HO2 H2O2 H2O AR
END
REACTIONS
H2+O2=2OH 1.7E13 0.0 47780.
OH+H2=H2O+H 1.17E9 1.3 3626. !D-L$W
H+O2=OH+O 5.13E16 -0.816 16507. !JAM,JCP 1981
O+H2=OH+H 1.8E10 1.0 8826.
H+O2+M=HO2+M 2.1E18 -1.0 0. !SLACK
H2O/21./ H2/3.3/ O2/0.0/
H+O2+O2=HO2+O2 6.7E19 -1.42 0. !SLACK,JAN
OH+HO2=H2O+O2 5.0E13 0.0 1000.
H+HO2=2OH 2.5E14 0.0 1900.
O+HO2=O2+OH 4.8E13 0.0 1000.
2OH=O+H2O 6.0E+8 1.3 0. !COHEN-WEST.
H2+M=H+H+M 2.23E12 0.5 92600.
H2O/6/ H/2/ H2/3/
O2+M=O+O+M 1.85E11 0.5 95560.
H+OH+M=H2O+M 7.5E23 -2.6 0.
H2O/20/
H+HO2=H2+O2 2.5E13 0.0 700.
HO2+HO2=H2O2+O2 2.0E12 0.0 0.
H2O2+M=OH+OH+M 1.3E17 0.0 45500.
H2O2+H=HO2+H2 1.6E12 0.0 3800.
H2O2+OH=H2O+HO2 1.0E13 0.0 1800.
END
基本结构是首先声明元素ELEMENTS，然后声明组分SPECIES，最后声明化学反应REACTIONS。三个部分的声明都必须以这三个关键字开头，然后以END结尾。对于声明的详细声明在说明文件里有特别详细的解释，这里不做进一步说明。除去化学反应相关的信息外，还需要热学参数的数据，即文件therm.dat，这部分非常长，只罗列出其中一部分：
THERMO
300.000 1000.000 5000.000
! GRI-Mech Version 3.0 Thermodynamics released 7/30/99
! NASA Polynomial format for CHEMKIN-II
! see README file for disclaimer
O L 1/90O 1 G 200.000 3500.000 1000.000 1
2.56942078E+00-8.59741137E-05 4.19484589E-08-1.00177799E-11 1.22833691E-15 2
2.92175791E+04 4.78433864E+00 3.16826710E+00-3.27931884E-03 6.64306396E-06 3
-6.12806624E-09 2.11265971E-12 2.91222592E+04 2.05193346E+00 4
O2 TPIS89O 2 G 200.000 3500.000 1000.000 1
3.28253784E+00 1.48308754E-03-7.57966669E-07 2.09470555E-10-2.16717794E-14 2
-1.08845772E+03 5.45323129E+00 3.78245636E+00-2.99673416E-03 9.84730201E-06 3
-9.68129509E-09 3.24372837E-12-1.06394356E+03 3.65767573E+00 4
END
热学参数是以THERMO开头，以END结尾的，另外可以通过THERMO ALL关键字控制是否调用LTHRM中内置的热学参数。同样说明文件中有非常详细的说明，而主要的热学参数来源可以参考这个文章
Kee, R. J., Rupley, F. M. and Miller, J. A.: “The CHEMKIN Thermodynamic Data Base,” Sandia
National Laboratories Report SAND87-8215B (1990).
将。经过翻译器翻译过后的输出文件chem.out如下：
**************************************************************
* CHEMKIN Collection Release 3.7 *
* CHEM Application *
* GAS-PHASE MECHANISM INTERPRETER *
* Copyright 1997-2002 Reaction Design. All Rights Reserved. *
**************************************************************
--------------------
ELEMENTS ATOMIC
CONSIDERED WEIGHT
--------------------
1. H 1.00797
2. O 15.9994
3. AR 39.9480
--------------------
-------------------------------------------------------------------------------
C
P H
H A
A R
SPECIES S G MOLECULAR TEMPERATURE ELEMENT COUNT
CONSIDERED E E WEIGHT LOW HIGH H O AR
-------------------------------------------------------------------------------
1. H2 G 0 2.01594 200 3500 2 0 0
2. O2 G 0 31.99880 200 3500 0 2 0
3. H G 0 1.00797 200 3500 1 0 0
4. O G 0 15.99940 200 3500 0 1 0
5. OH G 0 17.00737 200 3500 1 1 0
6. HO2 G 0 33.00677 200 3500 1 2 0
7. H2O2 G 0 34.01474 200 3500 2 2 0
8. H2O G 0 18.01534 200 3500 2 1 0
9. AR G 0 39.94800 300 5000 0 0 1
-------------------------------------------------------------------------------
(k = A T**b exp(-E/RT))
REACTIONS CONSIDERED A b E
1. H2+O2=2OH 1.70E+13 0.0 47780.0
2. OH+H2=H2O+H 1.17E+09 1.3 3626.0
3. H+O2=OH+O 5.13E+16 -0.8 16507.0
4. O+H2=OH+H 1.80E+10 1.0 8826.0
5. H+O2+M=HO2+M 2.10E+18 -1.0 0.0
H2O Enhanced by 2.100E+01
H2 Enhanced by 3.300E+00
O2 Enhanced by 0.000E+00
6. H+O2+O2=HO2+O2 6.70E+19 -1.4 0.0
7. OH+HO2=H2O+O2 5.00E+13 0.0 1000.0
8. H+HO2=2OH 2.50E+14 0.0 1900.0
9. O+HO2=O2+OH 4.80E+13 0.0 1000.0
10. 2OH=O+H2O 6.00E+08 1.3 0.0
11. H2+M=H+H+M 2.23E+12 0.5 92600.0
H2O Enhanced by 6.000E+00
H Enhanced by 2.000E+00
H2 Enhanced by 3.000E+00
12. O2+M=O+O+M 1.85E+11 0.5 95560.0
13. H+OH+M=H2O+M 7.50E+23 -2.6 0.0
H2O Enhanced by 2.000E+01
14. H+HO2=H2+O2 2.50E+13 0.0 700.0
15. HO2+HO2=H2O2+O2 2.00E+12 0.0 0.0
16. H2O2+M=OH+OH+M 1.30E+17 0.0 45500.0
17. H2O2+H=HO2+H2 1.60E+12 0.0 3800.0
18. H2O2+OH=H2O+HO2 1.00E+13 0.0 1800.0
NO ERRORS FOUND ON INPUT:
ASCII Vers. 1.2 CHEMKIN linkfile chem.asc written.
WORKING SPACE REQUIREMENTS ARE
INTEGER: 593
REAL: 441
CHARACTER: 12
CKLIB的主要功能
CKLIB.f中实现了和化学反应相关的几乎所有功能，库本身的函数变量此处进行简单的说明。首先是所有的相关函数均以CK开头避免和用户的子函数混淆，例如suboutine ckinit。
变量也分别给出了相应的命名，比如说压力P温度T质量分数Y摩尔分数X摩尔浓度C定压比热CP定容比热CV焓H熵S内能U等等。
子函数中均以CK开头，但是紧跟字母W的代表是计算反应速率相关的，CK,CD分别代表产生和消耗的速率的关键字，Q代表和反应平衡常数相关的子函数。
手册中详细的罗列了CKLIB中所有的子函数极其功能，此处不全部例举，但是要罗列出该模块能够实现的功能：
- 初始化：读取文件，输出数组长度，输出反应机理，查看是否为等离子体反应
- 元素相关：输出原子质量，输出元素对应的字符串数组，输出元素名的字符
- 组分相关：返回组分的电荷、字符串、元素种类数、组分标签、分子质量
- 反应机理相关：返回Arrhenius公式系数、摩尔反应热、组分和对应的系数、描述反应的字符串、某个反应的反应系数。输出是否为某个特殊反应。
- 气态常数和单位：返回通用气态常数和标准大气压
- 浓度和摩尔质量相关：计算平均分子质量、压力、密度、摩尔分数、质量分数、摩尔浓度或者摩尔分数。
- 无量纲热学性质相关：返回比热的拟合参数、定压比热、无量纲焓、可以计算的最高温。
- 带量纲热学性质相关：计算Cp,Cv，焓、熵、内能
- 摩尔单位下的热学性质：同上
- 质量单位下的热学性质：同上
- 摩尔单位下的平均热学性质：同上
- 化学反应速率相关：计算反应速率、组分的摩尔变化、生成速率、正逆反应速率、组分生成速率
- 平衡常数相关：计算平衡常数
算例conp.f
这是一个constant presssure的0维燃烧算例。依赖的是ode方程，形式如下：
给初始时刻各个组分的浓度，利用CKLIB计算反应速率，然后更新组分的质量分数，并利用反应放热和比热计算温度变化。这是一个具有一定刚性的ode问题，对于ode的求解，CHEMKIN中也提供了相应的计算模块vode.f。而这个ode方程的描述，写在subroutine fun中。这部分直接打印手册中的代码阅读即可，可读性较强。需要注意的是程序本身存在文件名大小写读取的bug，在运行之前需要进行调试。运行程序后，在控制台依次输入所需的参数
1 1000
H2 1
O2 3
N2 .1
END
3.0e-4 3.0e-5
相当于给定p为一个标准大气压，温度为1000K，组分和摩尔浓度初始值和步长和终止时间。之后输出的计算结果为：
T(SEC) TMP(K) H2 H O2 O OH
HO2 H2O2 H2O N N2
NO
0.000E+00 0.100E+04 0.244E+00 0.000E+00 0.732E+00 0.000E+00 0.000E+00
0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.244E-01
0.000E+00
0.300E-04 0.100E+04 0.244E+00 0.817E-05 0.732E+00 0.425E-05 0.144E-05
0.129E-04 0.103E-07 0.259E-04 0.181E-20 0.244E-01
0.375E-19
0.600E-04 0.196E+04 0.890E-02 0.169E-01 0.625E+00 0.570E-01 0.411E-01
0.174E-03 0.355E-04 0.224E+00 0.229E-09 0.262E-01
0.167E-07
0.900E-04 0.235E+04 0.367E-02 0.331E-02 0.658E+00 0.235E-01 0.392E-01
0.845E-04 0.445E-05 0.246E+00 0.193E-08 0.271E-01
0.163E-05
0.120E-03 0.243E+04 0.258E-02 0.185E-02 0.665E+00 0.165E-01 0.352E-01
0.693E-04 0.254E-05 0.251E+00 0.229E-08 0.272E-01
0.438E-05
0.150E-03 0.246E+04 0.216E-02 0.139E-02 0.669E+00 0.138E-01 0.330E-01
0.641E-04 0.197E-05 0.254E+00 0.236E-08 0.273E-01
0.730E-05
0.180E-03 0.248E+04 0.197E-02 0.120E-02 0.670E+00 0.125E-01 0.319E-01
0.619E-04 0.173E-05 0.255E+00 0.237E-08 0.273E-01
0.102E-04
0.210E-03 0.248E+04 0.188E-02 0.111E-02 0.671E+00 0.119E-01 0.313E-01
0.609E-04 0.162E-05 0.255E+00 0.238E-08 0.273E-01
0.131E-04
0.240E-03 0.249E+04 0.183E-02 0.106E-02 0.671E+00 0.116E-01 0.310E-01
0.604E-04 0.157E-05 0.256E+00 0.239E-08 0.273E-01
对应各个组分在随时间推进时的摩尔浓度变化。其中使用手册也对ode求解器进行了说明，subroutine vode中使用了计算刚性问题的算法，有需求可以进一步了解。
算例Premix.f
这是一个比较特别的算例，计算的是一个一维的定常燃烧问题，建议查看一下《层流预混燃烧的数值研究》秦凤华，东北大学。使用的就是这个代码模块里面的算法，并进行了更多的数值试验。算法计算的并非NS方程，而是针对当前定解，且一维的情形进行简化的火焰模型：
计算过程中假设工况如图：
预混气体从左侧以恒定的流量导入，经过预热区的加热后，在反应区燃烧。燃烧产物再从右侧排出。燃烧消耗的反应物速度和气体流入的速度一致，所以变成一个定常问题。这类定常问题在现在求解时，常常使用随时间推进的迭代算法。而当前程序为了满足计算结果定常，采用的方法是给定计算区域某点的温度为定值。而计算方法也并非随时间迭代，而是直接牛顿法迭代求解离散的代数方程组。求解的结果大体图像如下：
下面介绍Premix.f代码中的基本结构。参考另外一篇文献《PREMIX: a fortran program for modeling steady laminar one-dimensional premixed flames》Printed April 1998。程序需要依赖cklib，具体的程序依赖关系如图：
首先我们看到和conp.f一样，使用到了CHEMKIN INTERPRETER和CHEMKIN subroutine library，而另外引入的模块TRANSPORT PROPERTY FITTING和TRANSPORT subroutine library在这个文献中可以找到详细说明：
R. J. Kee, G. Dixon-Lewis, J. Warnatz, M. E. Coltrin, and J. A. Miller, “A Fortran Computer Code Package for the Evaluation of Gas-Phase Multicomponent Transport Properties,” Sandia National Laboratories Report SAND86-8246 (1986).
在组分的守恒方程中，存在组分之间的扩散系数。这一系数的计算需要借助图中右端的库函数。此处给出各个模块和CHEMKIN各个模块之间的对应关系：
Gas-phase Species and Kinetics Description ----> chem.inp
Thermodynamics Data Base ----> therm.dat
Transport Data Base ----> tran.data
CHEMKIN Interpreter ----> ckinterp.f
Transport Property Fitting ----> tranfit.f
ASCII CHEMKIN Output File ----> chem.out
CHEMKIN Linking file ----> chem.asc 注意此处是.asc不是.bin
Transport Linking file ----> tran.bin
ASCII Transport Output file ----> tran.out
CHEMKIN subroutine library ----> cklib1.f cklib2.f
Transport subroutine library ----> tranlib.f
Restart File ----> rest.bin 读取文件继续算之前没算完的部分
Save File ----> save.bin
Recover File ----> recov.bin 没算完的话存在这里
函数的调用过程比较简单：
PROGRAM PRDRIV
C*****precision > double
IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
C*****END precision > double
C*****precision > single
C IMPLICIT REAL (A-H,O-Z), INTEGER (I-N)
C*****END precision > single
PARAMETER (LENLWK=100, LENIWK=5000, LENRWK=200000, LENCWK=100,
1 LENSYM=16, LIN=5, LOUT=6, LRIN=14, LROUT=15,
2 LRCVR=16, LINCK=25, LINMC=35, NMAX=65, ZERO=0.0)
C LENLWK allocates the logical working space
C LENIWK allocates the integer working space
C LENRWK allocates the real working space
C LENCWK allocates the character working space
C LENSYM is the length of a character string
C LIN is the unit from which user input is read
C LOUT is the unit to which printed output is written
C LRIN is the unit from which the restart file is read
C LROUT is the unit to which the save file is written
C LRCVR is the unit to which the scratch save file is written
C LINCK is unit from which the Chemkin binary file is read
C LINTP is unit from which the Transport binary file is read
C NMAX is the total number of grid points allowed
DIMENSION IWORK(LENIWK), RWORK(LENRWK)
CHARACTER CWORK(LENCWK)*(LENSYM)
LOGICAL LWORK(LENLWK)
EXTERNAL CKTIME
TSTART = CKTIME(ZERO)
C
OPEN(LRIN,FORM='UNFORMATTED',STATUS='UNKNOWN',FILE='./rest.bin')
OPEN(LROUT,FORM='UNFORMATTED',STATUS='UNKNOWN',FILE='./save.bin')
OPEN(LRCVR,FORM='UNFORMATTED',STATUS='UNKNOWN',FILE='./recov.bin')
OPEN(LINCK,STATUS='UNKNOWN',FORM='FORMATTED',FILE='./chem.asc')
OPEN(LINMC,STATUS='UNKNOWN',FORM='FORMATTED', FILE='./tran.asc')
C
CALL PREMIX (NMAX, LIN, LOUT, LINCK, LINMC, LRIN, LROUT, LRCVR,
1 LENLWK, LWORK, LENIWK, IWORK, LENRWK, RWORK, LENCWK,
2 CWORK)
TEND = CKTIME (TSTART)
IF (TEND .GT. 60.) THEN
WRITE (LOUT, '(A,1PE15.2)') ' Total CPUtime (min): ',TEND/60.
ELSE
WRITE (LOUT, '(A,1PE15.2)') ' Total CPUtime (sec): ',TEND
ENDIF
CLOSE(LRIN)
CLOSE(LROUT)
CLOSE(LRCVR)
CLOSE(LINCK)
CLOSE(LINMC)
STOP
END
就是说另外写函数然后，之间在函数中调用这个子函数即可这个代码段对应文件Premxdem.f。计算过程中的求解器是一个Newton-Raphson方法，对应着文件Twopoints.f。变量说明在文献中有详细说明，这里不展开。文章介绍了两个算例，分别为氢氧燃烧和甲烷燃烧。对应的输入输出文件数据都很长，这里均不给出。对于当前算例，首先需要根据问题给定合适的输入文件chem.inp和therm.dat，以及tran.dat。通过两个翻译器分别生成两个linking file名为chem.asc和tarn.bin之外。还需要在运行后手动输入如下的代码段，作为控制关键字：
/ flame configuration, burner stabilized with specified temperature
BURN
TGIV
/ in the event of a Newton failure, take 100 timesteps of 1.E-6
TIME 100 1.00E-6
/ begin on a uniform mesh of 6 points
NPTS 6
/ definition of the computational interval
XEND 10.0
XCEN 5.0
WMIX 10.0
/ pressure and inlet mass flow rate
PRES 0.0329 (atmospheres)
FLRT 4.63E-3 (g/cm**2-sec)
/ adaptive mesh criteria
GRAD 0.2
CURV 0.5
/ unreacted mole fractions
MOLE
REAC O2 0.09
REAC AR .63
REAC H2 0.28
/ estimated products
PROD AR 0.68
PROD H2O 0.12
PROD H2 0.15
PROD OH 0.02
PROD O 0.02
PROD H 0.01
/ estimated intermediate mole fractions
INTM H2O2 .00001
INTM HO2 .001
INTM H2 .01
/ tolerances for the Newton iteration
ATOL 1.E-10
RTOL 1.E-4
/ tolerances for the time step Newton iteration
ATIM 1.E-5
RTIM 1.E-5
/ print control
PRNT 1
/ given temperature profile
TEMP 0. 373.7
TEMP .1250 484.5
TEMP .250 583.7
TEMP .375 672.2
TEMP .5 753.5
TEMP .75 901.4
TEMP 1.0 1027.
TEMP 1.25 1120.
TEMP 1.5 1184.
TEMP 2.0 1260.
TEMP 3.0 1348.
TEMP 6.0 1475.
TEMP 10.0 1524.
END
注意/中的说明文字也必须敲进去程序才能识别，空格多少个并不影响。当前调试好的Premix文件夹中，在控制台分别敲入如下命令：
make interp
make tranfit
make premix
./chem.exe
./tranfit.exe
./premix.exe
然后按顺序键入前面提到的控制命令，就可以输出结果，为文献中Page56的计算结果，如图：
MAKEFILE文件
CHEMKIN III中自带了一个MAKEFILE文件，内容如下：
# Make file for compile and link CHEMKIN Demo
# Compatible to Mircosoft Fortran Powerstation 4.0 or Digital Fortran 5.0
# Note: PCMACH.F is used
FC = DF
# Mircosoft Fortran Powerstation 4.0
FFLAGS = -4R8 -W0 -G5
# Compaq fortran 6.0
# FFLAGS = /align:rec8byte /integer_size:32 /real_size:64 /architecture:p5
LD = fl32
.F.o :
$(FC) $(FFLAGS) $<
INTPOBJ=CKINTERP.obj CKLIB1.OBJ CKLIB2.OBJ
EQUIOBJ=equi.obj eqlib.obj cklib1.obj cklib2.obj stanlib.obj
CONPOBJ=conp.obj cklib1.obj CKLIB2.OBJ vode.obj math.obj pcmach.obj
PREMIXOBJ=premix.obj premxdem.obj cklib1.obj cklib2.obj tranlib.obj twopnt.obj cputim.obj math.obj pcmach.obj
FITOBJ=fitdat.obj lsei.obj math.obj pcmach.obj
SENKINOBJ=senkin.obj senkdem.obj cklib1.obj cklib2.obj dasac.obj
PSROBJ=psr.obj psrdem.obj eqlib.obj stanlib.obj cklib1.obj cklib2.obj twopnt.obj cputim.obj math.obj pcmach.obj
TRANFITOBJ=tranfit.obj cklib1.obj cklib2.obj math.obj pcmach.obj
SHOCKOBJ=shock.obj cklib1.obj cklib2.obj vode.obj math.obj pcmach.obj
all: ckinterp EQUI CONP PREMIX FIT SENKIN PSR TRANFIT SHOCK
ckinterp: $(INTPOBJ)
$(LD) $(INTPOBJ)
EQUI: $(EQUIOBJ)
$(LD) $(EQUIOBJ)
CONP: $(CONPOBJ)
$(LD) $(CONPOBJ)
PSR: $(PSROBJ)
$(LD) $(PSROBJ)
PREMIX: $(PREMIXOBJ)
$(LD) $(PREMIXOBJ)
FIT: $(FITOBJ)
$(LD) $(FITOBJ)
SENKIN: $(SENKINOBJ)
$(LD) $(SENKINOBJ)
TRANFIT: $(TRANFITOBJ)
$(LD) $(TRANFITOBJ)
SHOCK: $(SHOCKOBJ)
$(LD) $(SHOCKOBJ)
说明我们当前的算例只是CHEMKIN中的其中两个模块。另外经过测试，当前这个Makefile文件并不能直接使用。文件中不同.f文件的依赖关系存在缺失，并且当前的编译器对输入输出文件以及.f文件的大小写是区分的，导致程序中有很多bug，而且代码段中部分format定义也存在bug，需要经过调试才能使用。
旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。
更多推荐
14
0- 0
已为社区贡献3条内容
所有评论(0)