在前面的帖子中已经大概给出了chemFoam主程序的结构，这里给出一个比较全面的总结
首先程序结构如下：
添加头文件
//*****************************//
int main(int argc, char *argv[])
{ 初始化
while(runTime.run())
{ 计算化学反应
输出中间结果
}
return 0;
}
我们接下来介绍各个部分的具体内容
while循环中的燃烧计算
其源码如下：
while (runTime.run())
{
//控制时间步长
#include "readControls.H"
#include "setDeltaT.H"
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
//燃烧的求解
#include "solveChemistry.H"
#include "YEqn.H"
#include "hEqn.H"
#include "pEqn.H"
//输出中间过程
#include "output.H"
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}
其中燃烧的过程的头文件源码内容如下：
//#include "solveChemistry.H"
dtChem = chemistry.solve(runTime.deltaT().value()); //求解之后的反应速率存在类中
scalar Qdot = chemistry.Qdot()()[0]/rho[0]; //反应放热速度？直接从类中提取
integratedHeat += Qdot*runTime.deltaT().value(); //反应放热
//#include "YEqn.H" 更新组分
{
forAll(Y, specieI)
{
volScalarField& Yi = Y[specieI];
solve(fvm::ddt(rho, Yi) - chemistry.RR(specieI), "Yi"); //代表 d(rho*Yi)/dt - 速率=0这个ode方程
}
}
//#include "hEqn.H" 更新焓，同时更新温度
if (constProp != "temperature") //非定温的情况下
{
volScalarField& h = thermo.he(); //可以计算焓
if (constProp == "volume") //定容和定压两种
{
h[0] = u0 + p[0]/rho[0] + integratedHeat;
}
else
{
h[0] = h0 + integratedHeat;
}
thermo.correct();
}
//#include "pEqn.H" 更新压力，密度本来就不变
{
rho = thermo.rho();
if (constProp == "volume")
{
scalar invW = 0.0;
forAll(Y, i)
{
invW += Y[i][0]/W[i];
}
Rspecific[0] = 1000.0*constant::physicoChemical::R.value()*invW;
p[0] = rho0*Rspecific[0]*thermo.T()[0]; //p = rho R T
rho[0] = rho0;
}
}
计算过程为，先计算反应速率，然后用反应速率计算组分，然后更新焓和温度，最后更新压力。其中组分和各个热学量的计算不难理解，而反应速率的计算调用了另外一个函数solve。该函数来自类StandardChemistryModel，内容如下：
template<class ReactionThermo, class ThermoType>
template<class DeltaTType>
Foam::scalar Foam::StandardChemistryModel<ReactionThermo, ThermoType>::solve
(
const DeltaTType& deltaT
)
{
BasicChemistryModel<ReactionThermo>::correct();
scalar deltaTMin = great;
if (!this->chemistry_)
{
return deltaTMin;
}
tmp<volScalarField> trho(this->thermo().rho());
const scalarField& rho = trho();
const scalarField& T = this->thermo().T();
const scalarField& p = this->thermo().p();
scalarField c0(nSpecie_);
forAll(rho, celli) //对所有的cell遍历
{
scalar Ti = T[celli];
if (Ti > Treact_) //超过燃点才会燃烧
{
const scalar rhoi = rho[celli];
scalar pi = p[celli];
for (label i=0; i<nSpecie_; i++) //通过密度和质量分数计算摩尔浓度
{
c_[i] = rhoi*Y_[i][celli]/specieThermo_[i].W();
c0[i] = c_[i];
}
// Initialise time progress
scalar timeLeft = deltaT[celli];
// Calculate the chemical source terms
while (timeLeft > small)
{
scalar dt = timeLeft;
this->solve(c_, Ti, pi, dt, this->deltaTChem_[celli]); //这里调用的是basicChemistryModel中纯虚函数的一个实现，更新了c_
timeLeft -= dt;
}
deltaTMin = min(this->deltaTChem_[celli], deltaTMin);
this->deltaTChem_[celli] =
min(this->deltaTChem_[celli], this->deltaTChemMax_);
for (label i=0; i<nSpecie_; i++)
{
RR_[i][celli] =
(c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
} //先计算出更新的摩尔浓度，再除以时间得到反映速率？
}
else
{
for (label i=0; i<nSpecie_; i++)
{
RR_[i][celli] = 0;
}
}
}
return deltaTMin;
}
其中RR_即为计算所得的反应速率。
所以chemFoam计算燃烧的过程如下：
while
{ 更新当前反应速率
利用反应速率计算反应热
更新组分
更新焓和温度
更新压力
}
头文件部分
首先给出一共添加了哪些头文件：
#include "fvCFD.H"
#include "rhoReactionThermo.H"
#include "BasicChemistryModel.H"
#include "reactingMixture.H"
#include "chemistrySolver.H"
#include "OFstream.H"
#include "thermoPhysicsTypes.H"
#include "basicSpecieMixture.H"
#include "cellModeller.H"
#include "thermoTypeFunctions.H"
其中fvCFD.H为有限体积法求解过程中需要使用的代码，它添加了大量的头文件。设计时间相关的类，优先体积法离散格式相关的类，以及单位制转换相关的类。求解过程中的solve(fvm::ddt(rho, Yi) - chemistry.RR(specieI), "Yi");依赖这部分代码。
rhoReactionThermo.H和reactingMixture.H以及baicSpecieMixture.H为多组分和热学相关的类，前面求解过程中焓温度以及压力的更新，依赖这些类中的函数。
BasicChemistryModel.H和chemistrySolver.H为燃烧求解过程使用的类，前面求解过程中的chemistry.solve()依赖这部分含糊，
OFstream.H为文件流和文件操作需要使用的头文件
thermoPhysicsTypes.H和thermoTypeFunctions.H进行了部分类型定义和函数定义，方便后面使用。
初始化部分
源码如下：
argList::noParallel(); //并行操作
#define CREATE_MESH createSingleCellMesh.H
#define NO_CONTROL
#include "postProcess.H" //续算时有用
#include "setRootCaseLists.H" //对终端输入args进行处理
#include "createTime.H" //初始化时间类runTime，同时从控制文件读取控制参数以进行初始化
#include "createSingleCellMesh.H" //创建一个只有一个cell的网格，正常的算例中此处为读取给定的网格
#include "createFields.H" //读取initialConditions，读取基础的热学信息，创建pChemistry
#include "createFieldRefs.H" //利用pChemistry创建chemistry
#include "readInitialConditions.H" //继续读取initialCondtions中信息，constantProperty以及各个组分初始浓度
#include "createControls.H" //继续读取controlDict中的信息，主要为时间步长的控制
主要创建的类为：
Foam::Time runTime
fvMesh mesh
autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh))
rhoReactionThermo& thermo = pThermo()
autoPtr<BasicChemistryModel<rhoReactionThermo>> pChemistry(BasicChemistryModel<rhoReactionThermo>::New(thermo))
OFstream post(args.path()/"chemFoam.out")
BasicChemistryModel<rhoReactionThermo>& chemistry = pChemistry()
分别时间，网格控制，以及热学和化学反应的控制。类对象在创建过程中就通过构造函数读取了控制文件进行了初始化。
执行方法
首先chemFoam是编译好的文件，我们只需要调用即可，将我们希望计算的控制文件配置好。例如我们希望计算氢气燃烧，就将配置好的文件夹添加到路径~.OpenFOAM/dyfluid-7/run/h2中：
键入chemFoam > chemFoam.out把结果保存起来方便后续和CHEMKIN的结果做对比：
之后将计算结果chemFoam.out复制到validation文件夹中，调用预先写好的createGraph脚本作图
脚本内容其实是gnuplot作图的代码
#!/bin/sh
if ! which gnuplot > /dev/null 2>&1
then
echo "gnuplot not found - skipping graph creation" >&2
exit 1
fi
gnuplot<<EOF
set terminal postscript eps color enhanced "Helvetica,20"
set output "OF_vs_CHEMKINII.eps"
set xlabel "Time / [s]" font "Helvetica,24"
set ylabel "Temperature / [K]" font "Helvetica,24"
set grid
set key left top
set xrange [0:0.001]
set yrange [750:2750]
set ytic 250
plot \\
"../chemFoam.out" u 1:2 t "OpenFOAM" with points lt 1 pt 6 ps 1.5,\\
"chemkinII" with lines title "Chemkin II" lt -1
EOF
#------------------------------------------------------------------------------
作图结果如下：
H2计算结果（上）和C8H18计算结果（下）
详细的类见关系解读
这里挖个坑，前面介绍的是程序大概的内容，以及其使用方法。但是代码实现使用的类型其实均经过高度的封装，这些类的支持性文件非常多，反而是openFOAM的主体源码内容。这里罗列大纲，并会不断的进行补充。
- 基础的类，如整型浮点型，字符串，布尔类型
- 基础的数据结构，如向量，数组，链表，域
- 文件读取过程中文件流的获取
- Chemistry一族的类间关系和可进行操作
- therm一族的类间关系和可进行操作
© 版权声明
文章版权归作者所有，未经允许请勿转载。
THE END
![[转]我国CAD软件产业亟待研究现状采取对策-卡核](https://www.caxkernel.com/wp-content/uploads/2024/07/frc-f080b20a9340c1a89c731029cb163f6a-212x300.png)
