OpenFOAM 燃烧化学求解过程解析

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

OF-7

UML 类图

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

OF-7 UML

我猜测: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 函数在 EulerImplicitode 中重新定义。所以这里的 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 的对象

OF-8

UML 类图

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

OF-8 UML!

我猜测: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 ReactionList<ThermoType> reactions_;

StandardChemistryModel 的构造函数中:

reactions_
(
dynamic_cast<const multiComponentMixture<ThermoType>&>(this->thermo()).species(),
specieThermo_,
this->mesh(),
*this
),

调用的是 ReactionList 的构造函数

template<class ThermoType>
Foam::ReactionList<ThermoType>::ReactionList
(
const speciesTable& species,
const PtrList<ThermoType>& speciesThermo,
const objectRegistry& ob,
const dictionary& dict
)
{
// Set general temperature limits from the dictionary
Reaction<ThermoType>::TlowDefault =
dict.lookupOrDefault<scalar>("Tlow", 0);

Reaction<ThermoType>::ThighDefault =
dict.lookupOrDefault<scalar>("Thigh", great);

const dictionary& reactions(dict.subDict("reactions"));

HashPtrTable<ThermoType> thermoDatabase;
forAll(speciesThermo, i)
{
thermoDatabase.insert
(
speciesThermo[i].name(),
speciesThermo[i].clone().ptr()
);
}

this->setSize(reactions.size());
label i = 0;

forAllConstIter(dictionary, reactions, iter)
{
this->set
(
i++,
Reaction<ThermoType>::New
(
species,
thermoDatabase,
ob,
reactions.subDict(iter().keyword())
).ptr()
);
}
}

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

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

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

在代码还没运行的时候,我们不知道 thermo 是 hePsiThermo 类型,那代码中为什么可以这样转换类型呢?

dynamic_cast<const multiComponentMixture<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);

在 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_;

BasicChemistryModel 中没有solve函数,为什么可以这样用呢:

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

原因在此:
ChemistryCombustion 的构造函数中,初始化了 chemistryPtr_(BasicChemistryModel<ReactionThermo>::New(thermo))
而这个 New 函数返回的又是 basicChemistryModel::New,这个 basicChemistryModel::New 函数才是真正的用于 RTS 的 Selector。

template<class ReactionThermo>
Foam::autoPtr<Foam::BasicChemistryModel<ReactionThermo>>
Foam::BasicChemistryModel<ReactionThermo>::New(ReactionThermo& thermo)
{
return basicChemistryModel::New<BasicChemistryModel<ReactionThermo>>
(
thermo
);
}

也就是说,chemistryPtr_ 实际上是指向最底层的 basicChemistryModel 的指针,而 basicChemistryModel 中定义了 solve 等纯虚函数,这样就能解释的通了。

basicChemistryModel::New

template
Foam::autoPtrFoam::basicChemistryModel::New
(
typename ChemistryModel::reactionThermo& thermo
)
{
IOdictionary chemistryDict
(
IOobject
(
thermo.phasePropertyName(“chemistryProperties”),
thermo.db().time().constant(),
thermo.db(),
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);

if (!chemistryDict.isDict("chemistryType"))
{
    FatalErrorInFunction
        << "Template parameter based chemistry solver selection is no "
        << "longer supported. Please create a chemistryType dictionary"
        << "instead." << endl << endl << "For example, the entry:" << endl
        << "    chemistrySolver ode<StandardChemistryModel<"
        << "rhoChemistryModel,sutherlandspecie<janaf<perfectGas>,"
        << "sensibleInternalEnergy>>" << endl << endl << "becomes:" << endl
        << "    chemistryType" << endl << "    {" << endl
        << "        solver ode;" << endl << "        method standard;"
        << endl << "    }" << exit(FatalError);
}

const dictionary& chemistryTypeDict =
    chemistryDict.subDict("chemistryType");

const word& solverName
(
    chemistryTypeDict.found("solver")
  ? chemistryTypeDict.lookup("solver")
  : chemistryTypeDict.found("chemistrySolver")
  ? chemistryTypeDict.lookup("chemistrySolver")
  : chemistryTypeDict.lookup("solver") // error if neither entry is found
);

const word& methodName
(
    chemistryTypeDict.lookupOrDefault<word>
    (
        "method",
        chemistryTypeDict.lookupOrDefault<bool>("TDAC", false)
      ? "TDAC"
      : "standard"
    )
);

dictionary chemistryTypeDictNew;
chemistryTypeDictNew.add("solver", solverName);
chemistryTypeDictNew.add("method", methodName);

Info<< "Selecting chemistry solver " << chemistryTypeDictNew << endl;

typedef typename ChemistryModel::thermoConstructorTable cstrTableType;
cstrTableType* cstrTable = ChemistryModel::thermoConstructorTablePtr_;

const word chemSolverCompThermoName =
    solverName + '<' + methodName + '<'
  + ChemistryModel::reactionThermo::typeName + ','
  + thermo.thermoName() + ">>";

typename cstrTableType::iterator cstrIter =
    cstrTable->find(chemSolverCompThermoName);

if (cstrIter == cstrTable->end())
{
    FatalErrorInFunction
        << "Unknown " << typeName_() << " type " << solverName << endl
        << endl;

    const wordList names(cstrTable->toc());

    wordList thisCmpts;
    thisCmpts.append(word::null);
    thisCmpts.append(word::null);
    thisCmpts.append(ChemistryModel::reactionThermo::typeName);
    thisCmpts.append(basicThermo::splitThermoName(thermo.thermoName(), 5));

    List<wordList> validNames;
    validNames.append(wordList(2, word::null));
    validNames[0][0] = "solver";
    validNames[0][1] = "method";
    forAll(names, i)
    {
        const wordList cmpts(basicThermo::splitThermoName(names[i], 8));

        bool isValid = true;
        for (label i = 2; i < cmpts.size() && isValid; ++ i)
        {
            isValid = isValid && cmpts[i] == thisCmpts[i];
        }

        if (isValid)
        {
            validNames.append(SubList<word>(cmpts, 2));
        }
    }

    FatalErrorInFunction
        << "All " << validNames[0][0] << '/' << validNames[0][1]
        << " combinations for this thermodynamic model are:"
        << endl << endl;
    printTable(validNames, FatalErrorInFunction);

    FatalErrorInFunction << endl;

    List<wordList> validCmpts;
    validCmpts.append(wordList(8, word::null));
    validCmpts[0][0] = "solver";
    validCmpts[0][1] = "method";
    validCmpts[0][2] = "reactionThermo";
    validCmpts[0][3] = "transport";
    validCmpts[0][4] = "thermo";
    validCmpts[0][5] = "equationOfState";
    validCmpts[0][6] = "specie";
    validCmpts[0][7] = "energy";
    forAll(names, i)
    {
        validCmpts.append(basicThermo::splitThermoName(names[i], 8));
    }

    FatalErrorInFunction
        << "All " << validCmpts[0][0] << '/' << validCmpts[0][1] << '/'
        << validCmpts[0][2] << "/thermoPhysics combinations are:"
        << endl << endl;
    printTable(validCmpts, FatalErrorInFunction);

    FatalErrorInFunction << exit(FatalError);
}

return autoPtr<ChemistryModel>(cstrIter()(thermo));

}

因此,实际调用的是派生类 StandardChemistryModel::solve(const scalarField& deltaT)。它返回的是实际干活的 solve 函数。

实际干活的 solve 函数:

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)
{
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)
//这里应该体现了 time split 算法,不对啊,这里只循环了一次啊。因为dt进入solve之后并没有修改!
//这样的话,运行一次之后,timeLeft就等于0了!
//上述表述适用于ode求解器,但是选用EulerImplicit的话,确实是在这里time split
//因为EulerImplicit中会更新dt,这样的话,这个while循环就会不止执行一次
//ode的time split体现在别处
{
scalar dt = timeLeft;
this->solve(pi, Ti, c_, celli, dt, this->deltaTChem_[celli]);
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;
}

在其中,调用了 6 参数的 solve 函数。6 参数的 solve 函数在继承自 StandardChemistryModelchemistrySolver 中被申明为纯虚函数,实际上在 EulerImplicitode 中重新定义。所以这里的 reaction->correct() 就相当于在这里求解化学反应!

但是化学反应怎么求解的,本文还没有涉及到。
下面以 ode 为例分析一下:
ode::solve 传进来的数据有 压力、温度、摩尔浓度、celli、dt、deltaTChem_。

ode::solve

在这个函数中,只有一个 const 参数,其余的传进来的数,都被该函数修改了!不对啊,deltaT 没有被修改啊!deltaT 传进了 5 参数的 ODESolver::solve 的第二个参数,但这个参数是 const 啊。
subDeltaT 被修改了。
p T c 肯定是被修改了的。

template<class ChemistryModel>
void Foam::ode<ChemistryModel>::solve
(
scalar& p,
scalar& T,
scalarField& c,
const label li,
scalar& deltaT,
scalar& subDeltaT
) const
{
// Reset the size of the ODE system to the simplified size when mechanism
// reduction is active
if (odeSolver_->resize())
{
odeSolver_->resizeField(cTp_);
}

const label nSpecie = this->nSpecie();

// Copy the concentration, T and P to the total solve-vector
for (int i=0; i<nSpecie; i++)
{
cTp_[i] = c[i];
}
cTp_[nSpecie] = T;
cTp_[nSpecie+1] = p;

odeSolver_->solve(0, deltaT, cTp_, li, subDeltaT);

for (int i=0; i<nSpecie; i++)
{
c[i] = max(0.0, cTp_[i]);
}
T = cTp_[nSpecie];
p = cTp_[nSpecie+1];
}

ode::solve 中又调用了 ODESolver::solve

ode 的 time split 体现在这里!!!

void Foam::ODESolver::solve
(
const scalar xStart,
const scalar xEnd,
scalarField& y,
const label li,
scalar& dxTry
) const
{
stepState step(dxTry);
scalar x = xStart;

for (label nStep=0; nStep<maxSteps_; nStep++)
// 这里是最大循环次数,但是中间可能退出循环
// maxSteps_(dict.lookupOrDefault<scalar>("maxSteps", 10000))
// 设置于chemistryProperties/odeCoeffs中,与relTol相同的位置
{
// Store previous iteration dxTry
scalar dxTry0 = step.dxTry;

step.reject = false;

// Check if this is a truncated step and set dxTry to integrate to xEnd
if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0)
{
step.last = true;
step.dxTry = xEnd - x;
}

// Integrate as far as possible up to step.dxTry
solve(x, y, li, step);

// Check if reached xEnd
if ((x - xEnd)*(xEnd - xStart) >= 0)
{
if (nStep > 0 && step.last)
{
step.dxTry = dxTry0;
}

dxTry = step.dxTry;

return;
}

step.first = false;

// If the step.dxTry was reject set step.prevReject
if (step.reject)
{
step.prevReject = true;
}
}

FatalErrorInFunction
<< "Integration steps greater than maximum " << maxSteps_ << nl
<< " xStart = " << xStart << ", xEnd = " << xEnd
<< ", x = " << x << ", dxDid = " << step.dxDid << nl
<< " y = " << y
<< exit(FatalError);
}

其中又调用了 solve(scalar& x,scalarField& y,const label li,stepState& step),这个函数在 ODESolver和它的派生类 seulex 中都定义了,那么这里应该是调用的是派生类的。不信的话,用 debug 模式验证一下。我验证过了,确实调用的 seulex::solve

seulex::solve

void Foam::seulex::solve
(
scalar& x,
scalarField& y,
const label li,
stepState& step
) const
{
temp_[0] = great;
scalar dx = step.dxTry;
y0_ = y;
dxOpt_[0] = mag(0.1*dx);

if (step.first || step.prevReject)
{
theta_ = 2*jacRedo_;
}

if (step.first)
{
// NOTE: the first element of relTol_ and absTol_ are used here.
scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
}

forAll(scale_, i)
{
scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
}

bool jacUpdated = false;

if (theta_ > jacRedo_)
{
odes_.jacobian(x, y, li, dfdx_, dfdy_);
jacUpdated = true;
}

int k;
scalar dxNew = mag(dx);
bool firstk = true;

while (firstk || step.reject)
{
dx = step.forward ? dxNew : -dxNew;
firstk = false;
step.reject = false;

if (mag(dx) <= mag(x)*sqr(small))
{
WarningInFunction
<< "step size underflow :" << dx << endl;
}

scalar errOld = 0;

for (k=0; k<=kTarg_+1; k++)
{
bool success = seul(x, y0_, li, dx, k, ySequence_, scale_);

if (!success)
{
step.reject = true;
dxNew = mag(dx)*stepFactor5_;
break;
}

if (k == 0)
{
y = ySequence_;
}
else
{
forAll(ySequence_, i)
{
table_[k-1][i] = ySequence_[i];
}
}

if (k != 0)
{
extrapolate(k, table_, y);
scalar err = 0;
forAll(scale_, i)
{
scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
err += sqr((y[i] - table_(0, i))/scale_[i]);
}
err = sqrt(err/n_);
if (err > 1/small || (k > 1 && err >= errOld))
{
step.reject = true;
dxNew = mag(dx)*stepFactor5_;
break;
}
errOld = min(4*err, 1);
scalar expo = 1.0/(k + 1);
scalar facmin = pow(stepFactor3_, expo);
scalar fac;
if (err == 0)
{
fac = 1/facmin;
}
else
{
fac = stepFactor2_/pow(err/stepFactor1_, expo);
fac = max(facmin/stepFactor4_, min(1/facmin, fac));
}
dxOpt_[k] = mag(dx*fac);
temp_[k] = cpu_[k]/dxOpt_[k];

if ((step.first || step.last) && err <= 1)
{
break;
}

if
(
k == kTarg_ - 1
&& !step.prevReject
&& !step.first && !step.last
)
{
if (err <= 1)
{
break;
}
else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4)
{
step.reject = true;
kTarg_ = k;
if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{
kTarg_--;
}
dxNew = dxOpt_[kTarg_];
break;
}
}

if (k == kTarg_)
{
if (err <= 1)
{
break;
}
else if (err > nSeq_[k + 1]*2)
{
step.reject = true;
if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{
kTarg_--;
}
dxNew = dxOpt_[kTarg_];
break;
}
}

if (k == kTarg_+1)
{
if (err > 1)
{
step.reject = true;
if
(
kTarg_ > 1
&& temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
)
{
kTarg_--;
}
dxNew = dxOpt_[kTarg_];
}
break;
}
}
}
if (step.reject)
{
step.prevReject = true;
if (!jacUpdated)
{
theta_ = 2*jacRedo_;

if (theta_ > jacRedo_ && !jacUpdated)
{
odes_.jacobian(x, y, li, dfdx_, dfdy_);
jacUpdated = true;
}
}
}
}

jacUpdated = false;

step.dxDid = dx;
x += dx;

label kopt;
if (k == 1)
{
kopt = 2;
}
else if (k <= kTarg_)
{
kopt=k;
if (temp_[k-1] < kFactor1_*temp_[k])
{
kopt = k - 1;
}
else if (temp_[k] < kFactor2_*temp_[k - 1])
{
kopt = min(k + 1, kMaxx_ - 1);
}
}
else
{
kopt = k - 1;
if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
{
kopt = k - 2;
}
if (temp_[k] < kFactor2_*temp_[kopt])
{
kopt = min(k, kMaxx_ - 1);
}
}

if (step.prevReject)
{
kTarg_ = min(kopt, k);
dxNew = min(mag(dx), dxOpt_[kTarg_]);
step.prevReject = false;
}
else
{
if (kopt <= k)
{
dxNew = dxOpt_[kopt];
}
else
{
if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
{
dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
}
else
{
dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
}
}
kTarg_ = kopt;
}

step.dxTry = step.forward ? dxNew : -dxNew;
}

其中调用了 ODESystem::jacobian(x, y, li, dfdx_, dfdy_)seulex::seul,还调用了 seulex::extrapolate

seulex::seul 中调用了 ODESystem::derivatives(xnew, y0, li, dy_)

但是 ODESystem 中只定义了三个纯虚函数 nEqnsderivativesjacobiannEqns 定义于 ODESolver,因为在 ODESolver 中,申明了 friend class ODESystem;。另外两个函数都定义在 StandardChemistryModel 中,因为 ODESystemStandardChemistryModel 的基类之一。

celli 传给了 li 这个参数,最终哪里用到了?
可能出现在
ODE/ 排除
thermophysicalModels/ 没找到。。。

官方说明
chemistryModel, reactions, ODESolver: Added “li” argument to all functions which evaluate the reaction rate

In chemistryModel “li” is set to the current cell index but for other reacting
systems it should be set to the current index of the element for which the
reaction system is being evaluated.

In the ODESolver “li” is the current index of the element for which the ODE
system is being solved if there is a list of related systems being solved,
otherwise it can be set to 0.

这两个函数很重要,代码如下:

template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
(
const scalar time,
const scalarField& c,
const label li,
scalarField& dcdt // 该函数用于计算 dcdt,包括 dc/dt dT/dt dp/dt 这里的c是摩尔浓度
) const
{
const scalar T = c[nSpecie_];
const scalar p = c[nSpecie_ + 1];

forAll(c_, i)
{
c_[i] = max(c[i], 0);
}

omega(p, T, c_, li, dcdt);
// dY/dt在这里计算,事实上,还是调用的Reaction::omega

// Constant pressure
// dT/dt = ...
scalar rho = 0;
scalar cSum = 0;
for (label i = 0; i < nSpecie_; i++)
{
const scalar W = specieThermo_[i].W();
cSum += c_[i];
rho += W*c_[i];
}
scalar cp = 0;
for (label i=0; i<nSpecie_; i++)
{
cp += c_[i]*specieThermo_[i].cp(p, T);
}
cp /= rho;

scalar dT = 0;
for (label i = 0; i < nSpecie_; i++)
{
const scalar hi = specieThermo_[i].ha(p, T);
dT += hi*dcdt[i];
}
dT /= rho*cp;

dcdt[nSpecie_] = -dT;

// dp/dt = ...
dcdt[nSpecie_ + 1] = 0;
}


template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::jacobian
(
const scalar t,
const scalarField& c,
const label li,
scalarField& dcdt,
scalarSquareMatrix& J
) const
{
const scalar T = c[nSpecie_];
const scalar p = c[nSpecie_ + 1];

forAll(c_, i)
{
c_[i] = max(c[i], 0);
}

J = Zero;
dcdt = Zero;

// To compute the species derivatives of the temperature term,
// the enthalpies of the individual species is needed
scalarField hi(nSpecie_);
scalarField cpi(nSpecie_);
for (label i = 0; i < nSpecie_; i++)
{
hi[i] = specieThermo_[i].ha(p, T);
cpi[i] = specieThermo_[i].cp(p, T);
}
scalar omegaI = 0;
List<label> dummy;
forAll(reactions_, ri)
{
const Reaction<ThermoType>& R = reactions_[ri];
scalar kfwd, kbwd;
R.dwdc(p, T, c_, li, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
//原来真正的Jacobian-J和derivative-dcdt计算在src/thermophysicalModels/specie/reaction/Reactions/Reaction/Reaction.C里边!
//这里dcdt其实和derivatives中计算的是一样的。不确定,需要验证一下。
R.dwdT(p, T, c_, li, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
}

// The species derivatives of the temperature term are partially computed
// while computing dwdc, they are completed hereunder:
scalar cpMean = 0;
scalar dcpdTMean = 0;
for (label i=0; i<nSpecie_; i++)
{
cpMean += c_[i]*cpi[i]; // J/(m^3 K)
dcpdTMean += c_[i]*specieThermo_[i].dcpdT(p, T);
}
scalar dTdt = 0.0;
for (label i=0; i<nSpecie_; i++)
{
dTdt += hi[i]*dcdt[i]; // J/(m^3 s)
}
dTdt /= -cpMean; // K/s

for (label i = 0; i < nSpecie_; i++)
{
J(nSpecie_, i) = 0;
for (label j = 0; j < nSpecie_; j++)
{
J(nSpecie_, i) += hi[j]*J(j, i);
}
J(nSpecie_, i) += cpi[i]*dTdt; // J/(mol s)
J(nSpecie_, i) /= -cpMean; // K/s/(mol/m3)
}

// ddT of dTdt
J(nSpecie_, nSpecie_) = 0;
for (label i = 0; i < nSpecie_; i++)
{
J(nSpecie_, nSpecie_) += cpi[i]*dcdt[i] + hi[i]*J(i, nSpecie_);
}
J(nSpecie_, nSpecie_) += dTdt*dcpdTMean;
J(nSpecie_, nSpecie_) /= -cpMean;
J(nSpecie_, nSpecie_) += dTdt/T;
}

OF 里的 ODE 是为了得到 deltaT 之后新的浓度 c,在给定当前的 c T p 的前提下。然后据此求解组分方程的源项:

for (label i=0; i<nSpecie_; i++)
{
RR_[i][celli] =
(c_[i] - c0[i])*specieThermo_[i].W()/deltaT[celli];
}

把 c T P 设为Φ\Phi,现在是已知Φn\Phi_n,要求Φn+1\Phi_{n+1}
我们可以根据化学反应来得到 dcdt 和 dTdt,即现在已知Φ=Φt=f(Φ)\Phi'=\dfrac{\partial \Phi}{\partial t}=f(\Phi),然后寻找合适的算法来求解Φn+1\Phi_{n+1}
最简单的方法就是Φn+1=Φn+hf(Φ)\Phi_{n+1}=\Phi_n + hf(\Phi)f(Φ)f(\Phi) 是斜率,h 是时间间隔。
那么现有先进的求解方法,应该也是基于这个思想。
比如 RK 方法,就是把时间间隔 h 分为若干小段。

6 阶 RK

yn+1=yn+h(c1k1+c2k2+c3k3+c4k4+c5k5+c6k6)y_{n+1}=y_n+h(c_1k_1+c_2k_2+c_3k_3+c_4k_4+c_5k_5+c_6k_6)

k1=f(tn,yn)=dydtk_1=f(t_n,y_n)=dydt

k2=f(tn+a2h,yn+hb21k1)k_2=f(t_n+a_2h,y_n+hb_{21}k_1)

k3=f(tn+a3h,yn+hb31k1+hb32k2)k_3=f(t_n+a_3h,y_n+hb_{31}k_1+hb_{32}k_2)

k4=f(tn+a4h,yn+hb41k1+hb42k2+hb43k3)k_4=f(t_n+a_4h,y_n+hb_{41}k_1+hb_{42}k_2+hb_{43}k_3)

k5=f(tn+a5h,yn+hb51k1+hb52k2+hb53k3+hb54k4)k_5=f(t_n+a_5h,y_n+hb_{51}k_1+hb_{52}k_2+hb_{53}k_3+hb_{54}k_4)

k6=f(tn+a6h,yn+hb6k1+hb62k2+hb63k3+hb64k4+hb65k5)k_6=f(t_n+a_6h,y_n+hb_{6}k_1+hb_{62}k_2+hb_{63}k_3+hb_{64}k_4+hb_{65}k_5)

OF-2.0.x 的代码,对应公式中的量:

代码中的变量公式中的变量
ak2_k2k_2
ak3_k3k_3
ak4_k4k_4
ak5_k5k_5
ak6_k6k_6
yyny_n
youtyn+1y_{n+1}
dydxk1k_1
第一个 yTemp_yn+hb21k1y_n+hb_{21}k_1
第二个 yTemp_yn+hb31k1+hb32k2y_n+hb_{31}k_1+hb_{32}k_2
第三个 yTemp_yn+hb41k1+hb42k2+hb43k3y_n+hb_{41}k_1+hb_{42}k_2+hb_{43}k_3
第四个 yTemp_yn+hb51k1+hb52k2+hb53k3+hb54k4y_n+hb_{51}k_1+hb_{52}k_2+hb_{53}k_3+hb_{54}k_4
第五个 yTemp_yn+hb61k1+hb62k2+hb63k3+hb64k4+hb65k5y_n+hb_{61}k_1+hb_{62}k_2+hb_{63}k_3+hb_{64}k_4+hb_{65}k_5

代码中返回的是:yn+1=yn+h(c1k1+c3k3+c4k4+c6k6)y_{n+1}=y_n+h(c_1k_1+c_3k_3+c_4k_4+c_6k_6),只需要保证ck=1\sum c_k =1 就行了。代码中刚好是相加为 1:

RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0,
RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0,

其中的 dydx 即 derivatives,但是没有用到 Jacobian?这个RK确实没有用到 Jacobian,但是其他的 ODE 求解器用到了。

KRR4

已知y=f(y)y'=f(y),我们这样得到yn+1=yn+hf(yn+)y_{n+1}=y_n+hf(y_{n+})

yn+1=yn+h[f(yn)+fyyn(yn+1yn)]y_{n+1}=y_n+h[f(y_n)+\dfrac{\partial f}{\partial y}|_{y_n}(y_{n+1}-y_n)]

yn+1=yh+h[1hfy]1f(yn)y_{n+1}=y_h+h[1-h\dfrac{\partial f}{\partial y}]^{-1} \cdot f(y_n)

y(x0+h)=y0+i=1sbikiy(x_0+h)=y_0+\sum_{i=1}^s b_ik_i

(1γhf)ki=hf(y0+j=1i1αijkj)+hfj=1i1γijkj,i=1,,s(1-\gamma h f) k_i = h f(y_0+\sum_{j=1}^{i-1}\alpha_{ij}k_j)+h f' \cdot \sum_{j=1}^{i-1}\gamma_{ij}k_j, \quad i=1,\cdots,s

gi=j=1i1γijkj+γkig_i=\sum_{j=1}^{i-1}\gamma_{ij}k_j+\gamma k_i

(1/γhf)g1=f(y0)(1/\gamma h -f')\cdot g_1=f(y_0)

(1/γhf)g2=f(y0+a21g1)+c21g1/h(1/\gamma h -f')\cdot g_2=f(y_0+a_{21}g_1)+c_{21}g_1/h

(1/γhf)g3=f(y0+a31g1+a32g2)+(c31g1+c32g2)/h(1/\gamma h -f')\cdot g_3=f(y_0+a_{31}g_1+a_{32}g_2)+(c_{31}g_1+c_{32}g_2)/h

(1/γhf)g4=f(y0+a41g1+a42g2+a43g3)+(c41g1+c42g2+c43g3)/h(1/\gamma h -f')\cdot g_4=f(y_0+a_{41}g_1+a_{42}g_2+a_{43}g_3)+(c_{41}g_1+c_{42}g_2+c_{43}g_3)/h

OF-2.0.x 的代码,对应公式中的量:

代码中的变量公式中的变量
dfdy_ff'
dfdx_fx\dfrac{\partial f}{\partial x}
g1_g1g_1
g2_g2g_2
g3_g3g_3
g4_g4g_4
c21c21c_{21}
c31c31c_{31}
c32c32c_{32}
c41c41c_{41}
c42c42c_{42}
c43c43c_{43}
a21a21a_{21}
a31a31a_{31}
a32a32a_{32}
b1b1b1
b2b2b2
b3b3b3
b4b4b4
KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0,
KRR4::b4 = 125.0/108.0,

b1+b2+b3+b4=4

seulex

OpenFOAM-dev

化学 ODE 求解理论

上边是代码追踪,下面来看一下化学 ODE 求解的理论知识。
待求解的量有:

Φ={T,Y1,Y2,,YNsp1}T\Phi = \left \lbrace T, Y_1, Y_2, \dotsc, Y_{N_{\text{sp}} - 1} \right \rbrace^{\text{T}}

YNsp=1k=1Nsp1YkY_{N_{\text{sp}}} = 1 - \sum_{k=1}^{N_{\text{sp}} - 1} Y_k

derivatives 即:

f=Φi=Φt={Tt,Y1t,Y2t,,YNsp1t}Tf = {\Phi_i'}=\frac{\partial \Phi}{\partial t} = \left \lbrace \frac{\partial T}{\partial t}, \frac{\partial Y_1}{\partial t}, \frac{\partial Y_2}{\partial t}, \dotsc, \frac{\partial Y_{N_{\text{sp}} - 1}}{\partial t} \right \rbrace^{\text{T}}

Tt=1ρcpk=1NsphkWkω˙k\frac{\partial T}{\partial t} = \frac{-1}{\rho c_p} \sum_{k=1}^{N_{\text{sp}}} h_k W_k \dot{\omega}_k

Ykt=1ρWkω˙kk=1,,Nsp1\frac{\partial Y_k}{\partial t} = \frac{1}{\rho} W_k \dot{\omega}_k \quad k = 1, \dotsc, N_{\text{sp}} - 1

Jacobian 矩阵:

Ji,j=fiΦj=Φi˙Φj\mathcal{J}_{i,j} = \frac{\partial f_i}{\partial \Phi_j} = \frac{\partial \dot{\Phi_i}}{\partial \Phi_j}

展开为:

[TTTY1TYNsp1Y1TY1Y1Y1YNsp1YNsp1TYNsp1Y1YNsp1YNsp1]\left[ \begin{array}{cccc} \dfrac{\partial {T'}}{\partial T} & \dfrac{\partial {T'}}{\partial Y_1} & \ldots & \dfrac{\partial {T'}}{\partial Y_{N_{\text{sp}} - 1}} \\ \dfrac{\partial {Y_1'}}{\partial T} & \dfrac{\partial {Y_1'}}{\partial Y_1} & \ldots & \dfrac{\partial {Y_1'}}{\partial Y_{N_{\text{sp}} - 1}} \\ \vdots & & \ddots & \vdots \\ \dfrac{\partial {Y}'_{N_{\text{sp}} - 1}}{\partial T} & \dfrac{\partial {Y}'_{N_{\text{sp}} - 1}}{\partial Y_1} & \ldots & \dfrac{\partial {Y}'_{N_{\text{sp}} - 1}}{\partial Y_{N_{\text{sp}} - 1}} \end{array} \right]

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 成长之路
您的肯定会给我更大动力~