OpenFOAM 燃烧化学求解过程解析

分析 OpenFOAM-dev 中的燃烧化学相关类,个人理解有限,希望帮我指正错误!

UML 类图

首先给出 UML 图,右键在新标签页中打开链接,可以获得更好的阅读体验。

Click to expand!

我猜测:ArrheniusReactionRate也可以直接作为ReversibleReaction的模板,而不通过FallOffReactionRate。

Reaction 的模板是 gasHThermoPhysicsReaction 的基类是 Foam::species::thermo,全称是 Foam::species::thermo< Foam::janafThermo<Foam::perfectGas<Foam::specie>>, Foam::sensibleEnthalpy>

ReversibleReaction 的模板是 ReactiongasHThermoPhysicsFallOffReactionRate。而 FallOffReactionRate 的模板是 ArrheniusReactionRateTroeFallOffFunction
ReversibleReaction 的基类是 Reaction<gasHThermoPhysics>

机理文件怎么读取的?
Reaction等相关类的对象在哪里创建的?
在StandardChemistryModel 中申明了反应的 list:

const PtrList<Reaction<ThermoType>>& reactions_;

StandardChemistryModel 的构造函数中:

reactions_
(
dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
)

this->thermo() 返回的是 BasicChemistryModel 中申明的 ReactionThermo& thermo_,这里模板 ReactionThermo 的实例是 psiReactionThermo
thermo_由构造函数中的thermo传入。

StandardChemistryModel的模板ReactionThermo
到底是psiReactionThermo,还是hePsiThermo?
是psiReactionThermom又如何,
是hePsiThermo又如何?
dynamic_cast<const reactingMixture&>(this->thermo())
二者之间,哪个可以转化为reactingMixture?答:hePsiThermo

有这个疑问是因为solver中创建thermo的时候用的是psiReactionThermo,但是后来实际上的对象是hePsiThermo类型的。

autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(mesh));
psiReactionThermo& thermo = pThermo();

const PtrList<Reaction<ThermoType>>& reactions_;
为什么可以这样初始化?没有找到对应的构造函数。
答案:
下边提到了,reactingMixture中实际上包括了从机理文件读取的反应的list,这里的目的是把读取的list转移到这里的reactions_中。

reactions_
(
dynamic_cast<const reactingMixture<ThermoType>&>(this->thermo())
)

StandardChemistryModel::calculateRR中使用了Reaction中的函数omega

const Reaction<ThermoType>& R = reactions_[ri];
const scalar omegai = R.omega(pi, Ti, c_, pf, cf, lRef, pr, cr, rRef);

chemkinReader(chemistryReader的派生类)中申明了另一个反应的list:

ReactionList<gasHThermoPhysics> reactions_;

并且定义了接口 reactions()

读取机理文件后,在chemkinReader::addReactionType中append 所有的reactions_

reactingMixture的构造函数
template<class ThermoType>
Foam::reactingMixture<ThermoType>::reactingMixture
(
const dictionary& thermoDict,
const fvMesh& mesh,
const word& phaseName
)
:
speciesTable(),
autoPtr<chemistryReader<ThermoType>>
(
chemistryReader<ThermoType>::New(thermoDict, *this)
),
multiComponentMixture<ThermoType>
(
thermoDict,
*this,
autoPtr<chemistryReader<ThermoType>>::operator()().speciesThermo(),
mesh,
phaseName
),
PtrList<Reaction<ThermoType>>
(
autoPtr<chemistryReader<ThermoType>>::operator()().reactions()
),
speciesComposition_
(
autoPtr<chemistryReader<ThermoType>>::operator()().specieComposition()
)
{
autoPtr<chemistryReader<ThermoType>>::clear();
}

在 dev 版本的 reactingFoam 中,创建了哪些和燃烧化学相关的对象呢?

autoPtr<psiReactionThermo> pThermo(psiReactionThermo::New(mesh));
psiReactionThermo& thermo = pThermo();
basicSpecieMixture& composition = thermo.composition();
PtrList<volScalarField>& Y = composition.Y();
autoPtr<CombustionModel<psiReactionThermo>> reaction
(
CombustionModel<psiReactionThermo>::New(thermo, turbulence())
);

reaction 如何使用呢?在 YEqn.H

reaction->correct();
Qdot = reaction->Qdot();

另外组分方程的源项也通过它调用:reaction->R(Yi)

结合上边的 UML 图,这里的 reaction 的类型实际上是 laminar
函数 correct() 的关键代码是:

this->chemistryPtr_->solve(this->mesh().time().deltaTValue());

其中的 chemistryPtr_ 定义于 laminar 的基类 ChemistryCombustion

autoPtr<BasicChemistryModel<ReactionThermo>> chemistryPtr_;
ChemistryCombustion 构造函数
template<class ReactionThermo>
Foam::ChemistryCombustion<ReactionThermo>::ChemistryCombustion
(
const word& modelType,
ReactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
)
:
CombustionModel<ReactionThermo>
(
modelType,
thermo,
turb,
combustionProperties
),
chemistryPtr_(BasicChemistryModel<ReactionThermo>::New(thermo))
{}
即,通过 `BasicChemistryModel` 调用其派生类 `StandardChemistryModel` 中的 `solve(const scalar deltaT)`。但是 `BasicChemistryModel` 中并没有定义这个 `solve` 函数,事实上,它是定义在基类 `basicChemistryModel` 中的纯虚函数。然后在 `StandardChemistryModel` 中,这个 `solve` 又调用了同一 `class` 中的实际干活的另一个 `solve(const DeltaTType& deltaT)` 函数,在其中,调用了同一 `class` 的 5 参数的 `solve` 函数。5 参数的 `solve` 函数在 `EulerImplicit` 和 `ode` 中重新定义。所以这里的 `reaction->correct()` 就相当于在这里求解化学反应!

laminar 中的函数 Qdot() 最初的纯虚函数定义于combustionModel,在 laminar 实际定义,返回值是 this->chemistryPtr_->Qdot()。因此跳转到了 BasicChemistryModel 类中,和上边的 solve 一样,它也是 basicChemistryModel 的纯虚函数,实际定义于 StandardChemistryModel,返回值是一个 tmp<Foam::volScalarField>,实际计算如下:

forAll(Y_, i)
{
forAll(Qdot, celli)
{
const scalar hi = specieThermo_[i].Hc();
Qdot[celli] -= hi*RR_[i][celli];
}
}

其中,Hc() 是该组分的化学焓(英文是 chemical enthalpy,又叫 formation enthalpy),公式是Δhf,k0\Delta h^0_{f,k},它的计算见 janafThermoI.H。另一篇博文正好给出了燃烧反应能量方程,其中显焓输运方程的反应源项,就是这里的 Qdot

函数 R(volScalarField& Y) 的返回值是 tmp<Foam::fvScalarMatrix>,我对 matrix 不太熟悉,这个矩阵应该是由此组成:this->chemistryPtr_->RR(specieI)。又来到了 BasicChemistryModel 类中,照旧,纯虚函数定义于 basicChemistryModel,实际定义于 StandardChemistryModel(在 StandardChemistryModelI.H 中定义),返回值是 RR_,它的类型是 PtrList<volScalarField::Internal>,代码中的描述是 List of reaction rate per specie [kg/m3/s]。

总结起来,就是在 solver 中,定义一个燃烧的对象:

autoPtr<CombustionModel<psiReactionThermo>> reaction
(
CombustionModel<psiReactionThermo>::New(thermo, turbulence())
);

然后在燃烧模型 ChemistryCombustion 里边,有一个化学模型的指针(Protected data,只能它自己或其派生类能够访问,它的对象不能访问)

autoPtr<BasicChemistryModel<ReactionThermo>> chemistryPtr_

所有关于化学的使用,都是通过它调用。比如求解化学反应方程,就是 reaction->correct()。而这个 correct 函数的内容是

this->chemistryPtr_->solve(this->mesh().time().deltaTValue());

调用了化学模型里边的 solve 函数。

template<class ReactionThermo>
Foam::ChemistryCombustion<ReactionThermo>::ChemistryCombustion
(
const word& modelType,
ReactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
)
:
CombustionModel<ReactionThermo>
(
modelType,
thermo,
turb,
combustionProperties
),
chemistryPtr_(BasicChemistryModel<ReactionThermo>::New(thermo))
{}

燃烧模型的构造函数里边,创建了化学模块的对象,传入了 thermo

三个函数完整定义
template<class ReactionThermo>
void Foam::combustionModels::laminar<ReactionThermo>::correct()
{
if (this->active())
{
if (integrateReactionRate_)
{
if (fv::localEulerDdt::enabled(this->mesh()))
{
const scalarField& rDeltaT =
fv::localEulerDdt::localRDeltaT(this->mesh());

if (this->coeffs().found("maxIntegrationTime"))
{
scalar maxIntegrationTime
(
readScalar(this->coeffs().lookup("maxIntegrationTime"))
);

this->chemistryPtr_->solve
(
min(1.0/rDeltaT, maxIntegrationTime)()
);
}
else
{
this->chemistryPtr_->solve((1.0/rDeltaT)());
}
}
else
{
this->chemistryPtr_->solve(this->mesh().time().deltaTValue());
}
}
else
{
this->chemistryPtr_->calculate();
}
}
}


template<class ReactionThermo>
Foam::tmp<Foam::fvScalarMatrix>
Foam::combustionModels::laminar<ReactionThermo>::R(volScalarField& Y) const
{
tmp<fvScalarMatrix> tSu(new fvScalarMatrix(Y, dimMass/dimTime));

fvScalarMatrix& Su = tSu.ref();

if (this->active())
{
const label specieI =
this->thermo().composition().species()[Y.member()];

Su += this->chemistryPtr_->RR(specieI);
}

return tSu;
}


template<class ReactionThermo>
Foam::tmp<Foam::volScalarField>
Foam::combustionModels::laminar<ReactionThermo>::Qdot() const
{
tmp<volScalarField> tQdot
(
new volScalarField
(
IOobject
(
this->thermo().phasePropertyName(typeName + ":Qdot"),
this->mesh().time().timeName(),
this->mesh(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
this->mesh(),
dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
)
);

if (this->active())
{
tQdot.ref() = this->chemistryPtr_->Qdot();
}

return tQdot;
}

RTS

declareRunTimeSelectionTable

首先在源代码中全局搜索发现 declareRunTimeSelectionTable 函数在 CombustionModel.H 中被调用。

declareRunTimeSelectionTable
(
autoPtr,
CombustionModel,
dictionary,
(
const word& modelType,
ReactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
),
(modelType, thermo, turb, combustionProperties)
);

而这个函数定义于 \src\OpenFOAM\db\runTimeSelection\construction\runTimeSelectionTables.H 中:

#define declareRunTimeSelectionTable(autoPtr,baseType,argNames,argList,parList)\
typedef autoPtr<baseType> (*argNames##ConstructorPtr)argList; \
typedef HashTable<argNames##ConstructorPtr, word, string::hash> \
argNames##ConstructorTable; \
static argNames##ConstructorTable* argNames##ConstructorTablePtr_; \
static void construct##argNames##ConstructorTables(); \
template<class baseType##Type> \
class add##argNames##ConstructorToTable \
{ \
public: \
static autoPtr<baseType> New argList \
{ \
return autoPtr<baseType>(new baseType##Type parList); \
} \
\
add##argNames##ConstructorToTable \
( \
const word& lookup = baseType##Type::typeName \
) \
{ \
construct##argNames##ConstructorTables(); \
argNames##ConstructorTablePtr_->insert(lookup, New); \
} \
}; \
  • 定义一个函数指针,这样定义的结果是, dictionaryConstructorPtr 代表一个指向参数为 argList,返回类型为 autoPtr< baseType> 的函数指针。
  • 将一个 key 和 value 分别为 worddictionaryConstructorPtrHashTable 重命名为 argNames##ConstructorTable,即 dictionaryConstructorTable,并创建了一个以它为类型的指针 argNames##ConstructorTablePtr_,即 dictionaryConstructorTablePtr_
  • 申明一个类:add##argNames##ConstructorToTable,即 adddictionaryConstructorToTable,它需要模板 baseType##Type,其中的 baseTypeCombustionModel
  • 类中定义了 New 函数,argList 是其形参列表,parList 是返回的 baseType##Type 的构造函数的形参列表。注意,这里调用了基类的构造函数!!!
  • 但是如果我们在派生类中创建这个 add 类的对象的话,它的模板则是派生类,就是说 New 函数的返回值是指向派生类对象基类指针。
  • add##argNames##ConstructorToTable 的构造函数中,向 dictionaryConstructorTablePtr_insert 了一组 key and value:
    key 是 typeName,vaule 是 返回 基类 autoPtrNew 函数。

defineTemplateRunTimeSelectionTable

defineTemplateRunTimeSelectionTable 函数被发现调用于 makeCombustionTypes.H 中的 makeCombustion 函数中。
makeCombustion 函数在 CombustionModels.C 中被调用:

makeCombustion(psiReactionThermo);

具体看看这个函数:

#define makeCombustion(Comp)                                                   \
\
typedef CombustionModel<Comp> CombustionModel##Comp; \
\
defineTemplateTypeNameAndDebugWithName \
( \
CombustionModel##Comp, \
( \
word(CombustionModel##Comp::typeName_()) + "<" + Comp::typeName \
+ ">" \
).c_str(), \
0 \
); \
\
defineTemplateRunTimeSelectionTable \
( \
CombustionModel##Comp, \
dictionary \
);

首先 defineTemplateTypeNameAndDebugWithName 函数创建了一个 word 名字是 typeName 值是 Name,即 laminar<psiReactionThermo>

其次是调用 defineTemplateRunTimeSelectionTable 函数,它定义于 \src\OpenFOAM\db\runTimeSelection\construction\runTimeSelectionTables.H 中。

#define defineTemplateRunTimeSelectionTable(baseType,argNames)                 \
defineRunTimeSelectionTablePtr(baseType,argNames); \
defineRunTimeSelectionTableConstructor(baseType,argNames); \
defineRunTimeSelectionTableDestructor(baseType,argNames)

这里的 baseTypeCombustionModel##psiReactionThermoargNamesdictionary
包括三个宏函数,第一个函数定义了一个类型为 dictionary##ConstructorTable (即 HashTable 的别名)的名字为 dictionary##ConstructorTablePtr_ 的空指针。

#define defineRunTimeSelectionTablePtr(baseType,argNames)                      \
baseType::argNames##ConstructorTable* \
baseType::argNames##ConstructorTablePtr_ = nullptr

第二个函数定义了一个函数,名为 constructdictionaryConstructorTables(),作用是创建了一个类型为 dictionaryConstructorTable,名字为 dictionaryConstructorTablePtr_ 的对象。

#define defineRunTimeSelectionTableConstructor(baseType,argNames)              \
void baseType::construct##argNames##ConstructorTables() \
{ \
baseType::argNames##ConstructorTablePtr_ \
= new baseType::argNames##ConstructorTable; \
}

第三个函数定义了一个函数(后来在 add##argNames##ConstructorToTable 的析构函数中被调用),名字是 destroydictionaryConstructorTables()

创建对象

我们在创建燃烧模型的对象时(createFields.H)是这样的:

Info<< "Creating reaction model\n" << endl;
autoPtr<CombustionModel<psiReactionThermo>> reaction
(
CombustionModel<psiReactionThermo>::New(thermo, turbulence())
);

调用的是 CombustionModel 中的 New

//- Selector
static autoPtr<CombustionModel> New
(
ReactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties=combustionPropertiesName //默认值是 combustionPropertiesName
);

定义于它的 C 文件中,返回类型是 autoPtr<CombustionModel>CombustionModel 需要一个模板,我们给它的是 psiReactionThermo

template<class ReactionThermo>
Foam::autoPtr<Foam::CombustionModel<ReactionThermo>>
Foam::CombustionModel<ReactionThermo>::New
(
ReactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
)
{
return
combustionModel::New<CombustionModel<ReactionThermo>>
// 这里发现这个 New 调用时需要给它模板,我们给的就是 CombustionModel<psiReactionThermo>
(
thermo,
turb,
combustionProperties
);
}

调用了它的基类 combustionModel(但没有在 C 文件中定义)的 New 函数:

// Selectors

//- Generic New for each of the related chemistry model
template<class CombustionModel>
static autoPtr<CombustionModel> New
(
typename CombustionModel::reactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
);

上边这个 New 函数发现是在 combustionModelTemplates.C 中定义的。
这里的 CombustionModel 就是 CombustionModelreactionThermo 是它的模板 ReactionThermo,我们在这里指定它为 psiReactionThermo

template<class CombustionModel>
Foam::autoPtr<CombustionModel> Foam::combustionModel::New
(
typename CombustionModel::reactionThermo& thermo,
const compressibleTurbulenceModel& turb,
const word& combustionProperties
)
{
word combModelName("none");
IOdictionary(combIO).lookup("combustionModel") >> combModelName;
Info<< "Selecting combustion model " << combModelName << endl;
// 字典 combustionProperties 里写的那个名字,如 laminar

typedef typename CombustionModel::dictionaryConstructorTable cstrTableType;
//这是在 declareRunTimeSelectionTable 中的那个 HashTable
cstrTableType* cstrTable = CombustionModel::dictionaryConstructorTablePtr_;
//dictionaryConstructorTablePtr_ 在 declareRunTimeSelectionTable 申明过
// 并且在 constructdictionaryConstructorTables() 这个函数中分派了实际指向的对象。
//这个函数申明于 declareRunTimeSelectionTable,定义于 defineTemplateRunTimeSelectionTable
//而这个函数又在 adddictionaryConstructorToTable 的构造函数中被调用
//adddictionaryConstructorToTable,即 add##argNames##ConstructorToTable 在 declareRunTimeSelectionTable 中定义

const word compCombModelName =
combModelName + '<' + CombustionModel::reactionThermo::typeName + '>';

const word thermoCombModelName =
combModelName + '<' + CombustionModel::reactionThermo::typeName + ','
+ thermo.thermoName() + '>';

typename cstrTableType::iterator compCstrIter =
cstrTable->find(compCombModelName);
//在哈希表中查询,key 是 laminar<psiReactionThermo> ,返回值是 返回指向基类autoPtr的New函数
typename cstrTableType::iterator thermoCstrIter =
cstrTable->find(thermoCombModelName);

return autoPtr<CombustionModel>
(
thermoCstrIter != cstrTable->end()
? thermoCstrIter()(combModelName, thermo, turb, combustionProperties)
: compCstrIter()(combModelName, thermo, turb, combustionProperties)
//modeltype 从这里由 combModelName 传入!!!如 laminar
);
}

里边用到了一个类:CombustionModel::dictionaryConstructorTable,它是在 declareRunTimeSelectionTable 中定义的,就是那个哈希表的别名。

typedef HashTable<argNames##ConstructorPtr, word, string::hash>
argNames##ConstructorTable; //argNames 是 dictionary

因为是在 CombustionModel.H 中调用的 declareRunTimeSelectionTable,于是 CombustionModel::dictionaryConstructorTable 就被定义好了。

laminars.C

makeCombustionTypes(laminar, psiReactionThermo);

而这个函数定义于 makeCombustionTypes.H:

#define makeCombustionTypes(CombModel, Comp)                                   \
\
typedef combustionModels::CombModel<Comp> CombModel##Comp; \
\
defineTemplateTypeNameAndDebugWithName \
( \
CombModel##Comp, \
( \
word(CombModel##Comp::typeName_()) + "<" + Comp::typeName + ">" \
).c_str(), \
0 \
); \
\
CombustionModel<Comp>:: \
add##dictionary##ConstructorToTable<CombModel##Comp> \
add##CombModel##Comp##dictionary##ConstructorTo##CombustionModel##Comp\
##Table_;
//这句话的作用和 addToRunTimeSelectionTable 一致,都是创建一个 add 类的对象,创建它的同时,将返回 指向派生类对象的基类指针 的 New 函数插进了 HashTable,即 dictionaryConstructorTable。后续就是在这个里边 find 我们在字典文件中指定的派生类。

CombModel 是 laminar,Comp 是 psiReactionThermo
combustionModels::laminar<psiReactionThermo> 重命名为 laminar##psiReactionThermo
CombustionModel<psiReactionThermo> 的命名空间下,
创建了类型为 adddictionaryConstructorToTable<laminar##psiReactionThermo> 的对象。于是在这里调用了 adddictionaryConstructorToTable 的构造函数!!!倒退至前边所述,这个构造函数中调用了 constructdictionaryConstructorTables(),作用是给 dictionaryConstructorTablePtr_ new 一个对象,类型是 dictionaryConstructorTable ,即 HashTable。并且构造函数中还向该dictionaryConstructorTablePtr_ insert 了 New 函数:
这个函数 new 了一个对象,并将其绑定至基类的指针,该对象的类型是 baseType##Type,即 add 这个类的模板的类型,而我们在这里给它的模板是 laminar##psiReactionThermo。就是说,我们在派生类中创建add这个类的对象时,同时向哈希表中插入了一组key and value,即派生类的typeName和返回指向派生类对象的基类指针的New函数!

总结

declareRunTimeSelectionTable 创建了一个哈希表的指针 argNames##ConstructorTablePtr_,但是没有给它分配对象。只有当创建第一个 add##argNames##ConstructorToTable<baseType##Type> 对象时,才会 new 一个baseType::argNames##ConstructorTable 对象给它指向。这是这个函数的功能:construct##argNames##ConstructorTables(申明于 declareRunTimeSelectionTable(autoPtr,baseType,argNames,argList,parList),定义于defineRunTimeSelectionTableConstructor(baseType,argNames))。
此外,创建每个 add##argNames##ConstructorToTable<baseType##Type> 对象时,都会调用它的构造函数,该构造函数数向 argNames##ConstructorTablePtr_ 插入一组 key and value:baseType##Type::typeNameNew 函数。
说起这个 New 函数,那就大有来头,它返回基类指针,但是指向的对象类型是 baseType##Type,即add##argNames##ConstructorToTable<baseType##Type> 这个类的模板。
我们经常在派生类中创建一个 add 的对象,如:laminars.C 中的 makeCombustionTypes 中的:

CombustionModel<Comp>::add##dictionary##ConstructorToTable<CombModel##Comp>
add##CombModel##Comp##dictionary##ConstructorTo##CombustionModel##Comp##Table_;

模板是 laminar##psiReactionThermo,即往哈希表中插入的 New 函数,返回值是基类指针,指向的是 laminar##psiReactionThermo 的对象

文章作者: Yan Zhang
文章链接: https://openfoam.top/chemistry/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路
微信打赏给博主更多动力吧~