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` 函数在 `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 的对象

OF-dev

UML 类图

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

OF-dev 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 的基类之一。

这两个函数很重要,代码如下:
template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
(
const scalar time,
const scalarField& c,
const label li,
scalarField& dcdt
) 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);
//dcdt在这里计算,事实上,还是调用的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_
g2_
g3_
g4_
c21c21c_{21}
c31c31c_{31}
c32c32c_{32}
c41c41c_{41}
c42c42c_{42}
c43c43c_{43}
a21
a31
a32
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

上边是代码追踪,下面来看一下化学 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 成长之路
微信打赏给博主更多动力吧~