openFOAM学习笔记(五)——chemFoam的运行过程

3780

在前面的帖子中已经大概给出了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的主体源码内容。这里罗列大纲，并会不断的进行补充。

展开阅读全文
© 版权声明
文章版权归作者所有，未经允许请勿转载。
THE END
综合
# openFOAM
喜欢就支持一下吧
点赞0 分享
Work hard in silence, let success make the noise.
在沉默中努力，让成功自己发声
ycwang125的头像-卡核
相关推荐
文章目录