OpenFOAM 燃烧化学求解过程解析

2018-06-08| code analysis|字数总计: 9.3k|阅读时长: 46 分钟

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

OF-7

UML 类图

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

OF-7 UML

Reaction

疑问：

机理文件怎么读取的？
Reaction 等相关类的对象在哪里创建的？

尝试解答：

第一个问题：

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();
}

调用 speciesTable 的默认构造函数，它来自 basicMultiComponentMixture::species_。

调用 chemistryReader 的 selector，按照字典构造用户想要的对象（默认是 chemkinReader）。
比如调用 chemkinReader 的构造函数：

Foam::chemkinReader::chemkinReader
(
    const dictionary& thermoDict,
    speciesTable& species
)
:
    lineNo_(1),
    specieNames_(10),
    speciesTable_(species),
    reactions_(speciesTable_, speciesThermo_),
    newFormat_(thermoDict.lookupOrDefault("newFormat", false)),
    imbalanceTol_(thermoDict.lookupOrDefault("imbalanceTolerance", rootSmall))
{
    fileName chemkinFile(fileName(thermoDict.lookup("CHEMKINFile")).expand());
    fileName thermoFile = fileName::null;
    if (thermoDict.found("CHEMKINThermoFile"))
    {
        thermoFile = fileName(thermoDict.lookup("CHEMKINThermoFile")).expand();
    }
    fileName transportFile
    (
        fileName(thermoDict.lookup("CHEMKINTransportFile")).expand()
    );

    read(chemkinFile, thermoFile, transportFile);
}

主要是 read 这个函数：

Click to expand!

void Foam::chemkinReader::read(const fileName& CHEMKINFileName, const fileName& thermoFileName, const fileName& transportFileName)
{
    Reaction<gasHThermoPhysics>::TlowDefault = 0;
    Reaction<gasHThermoPhysics>::ThighDefault = great;

    transportDict_.read(IFstream(transportFileName)());

    if (thermoFileName != fileName::null)
    {
        std::ifstream thermoStream(thermoFileName.c_str());

        if (!thermoStream)
        {
            FatalErrorInFunction
                << "file " << thermoFileName << " not found"
                << exit(FatalError);
        }

        yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
        yy_switch_to_buffer(bufferPtr);

        while (lex() != 0)
        {}

        yy_delete_buffer(bufferPtr);

        lineNo_ = 1;
    }

    std::ifstream CHEMKINStream(CHEMKINFileName.c_str());

    yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
    yy_switch_to_buffer(bufferPtr);

    // 在这里通过关键字将反应分类
    initReactionKeywordTable();
    // 然后 lex 处理时根据关键字来读取不同类别的反应
    // 根据中间符号来判断可逆与不可逆
    // reversibleReactionDelimiter {space}("="|"<=>"){space}
    // irreversibleReactionDelimiter {space}"=>"{space}
    /*
    enum reactionType
    {
        irreversible,
        reversible,
        nonEquilibriumReversible,
        unknownReactionType
    };
    const char* Foam::chemkinReader::reactionTypeNames[4] =
    {
        "irreversible",
        "reversible",
        "nonEquilibriumReversible",
        "unknownReactionType"
    };


   enum reactionRateType
    {
        Arrhenius,
        thirdBodyArrhenius,
        unimolecularFallOff,
        chemicallyActivatedBimolecular,
        LandauTeller,
        Janev,
        powerSeries,
        unknownReactionRateType
    };
    const char* Foam::chemkinReader::reactionRateTypeNames[8] =
    {
        "Arrhenius",
        "thirdBodyArrhenius",
        "unimolecularFallOff",
        "chemicallyActivatedBimolecular",
        "LandauTeller",
        "Janev",
        "powerSeries",
        "unknownReactionRateType"
    };

    enum fallOffFunctionType
    {
        Lindemann,
        Troe,
        SRI,
        unknownFallOffFunctionType
    };
    const char* Foam::chemkinReader::fallOffFunctionNames[4] =
    {
        "Lindemann",
        "Troe",
        "SRI",
        "unknownFallOffFunctionType"
    };
    */

    /*
    reactionKeywordTable_.insert("M", thirdBodyReactionType);
    reactionKeywordTable_.insert("LOW", unimolecularFallOffReactionType);
    reactionKeywordTable_.insert
    (
        "HIGH",
        chemicallyActivatedBimolecularReactionType
    );
    reactionKeywordTable_.insert("TROE", TroeReactionType);
    reactionKeywordTable_.insert("SRI", SRIReactionType);
    reactionKeywordTable_.insert("LT", LandauTellerReactionType); // -->LandauTeller
    reactionKeywordTable_.insert("RLT", reverseLandauTellerReactionType);
    reactionKeywordTable_.insert("JAN", JanevReactionType);
    reactionKeywordTable_.insert("FIT1", powerSeriesReactionRateType);
    reactionKeywordTable_.insert("HV", radiationActivatedReactionType);
    reactionKeywordTable_.insert("TDEP", speciesTempReactionType);
    reactionKeywordTable_.insert("EXCI", energyLossReactionType);
    reactionKeywordTable_.insert("MOME", plasmaMomentumTransfer);
    reactionKeywordTable_.insert("XSMI", collisionCrossSection);
    reactionKeywordTable_.insert("REV", nonEquilibriumReversibleReactionType);
    reactionKeywordTable_.insert("DUPLICATE", duplicateReactionType);
    reactionKeywordTable_.insert("DUP", duplicateReactionType);
    reactionKeywordTable_.insert("FORD", speciesOrderForward);
    reactionKeywordTable_.insert("RORD", speciesOrderReverse);
    reactionKeywordTable_.insert("UNITS", UnitsOfReaction);
    reactionKeywordTable_.insert("END", end);
    */
}

举例看看 chemkinLexer 中发生了什么：
（我们在这里说明的时候，下标都是从0开始的）

Click to expand!

// Read species
 if (!specieIndices_.found(specieName))
 {
     specieNames_.append(specieName);
     specieIndices_.insert(specieName, currentSpecieIndex++);
 }


 // Read thermo
 // 猜测是从第[10]个字符开始读，长度是10个字符
 allCommonT = stringToScalar(temperaturesString(10, 10));


 // 第一个空格之前的字符，当做组分名字，长度限制：18
 // thermoSpecieName      .{18}
 size_t spacePos = specieString.find(' ');
 if (spacePos != string::npos)
 {
     currentSpecieName = specieString(0, spacePos);
 }
 else
 {
     currentSpecieName = specieString;
 }

 // 接下来6个字符随便写什么，会被忽略：[18]-[23]


 // [24]-[43] 列给出组分原子及其原子组成数
 // thermoFormula         .{20}
 // 这 20 列均分为每 5 列一段，共四段
 // 每段又分为 1 个两列和 1 个三列，分别用来写原子名和原子数
 nSpecieElements = 0;
 currentSpecieComposition.setSize(5);

 for (int i=0; i<4; i++)
 {
     word elementName(foamName(thermoFormula(5*i, 2)));
     label nAtoms = atoi(thermoFormula(5*i + 2, 3).c_str());

     if (elementName.size() && nAtoms)
     {
         correctElementName(elementName);
         currentSpecieComposition[nSpecieElements].name() =
             elementName;
         currentSpecieComposition[nSpecieElements++].nAtoms() = nAtoms;
     }
 }


 //第 [44] 列上用（G 、S、 L）指定组分相态
 char phaseChar = YYText()[0]; // 这里应该是表示继续往前读一列

 switch (phaseChar)
 {
     case 'S':
         speciePhase_.insert(currentSpecieName, solid);
     break;

     case 'L':
         speciePhase_.insert(currentSpecieName, liquid);
     break;

     case 'G':
         speciePhase_.insert(currentSpecieName, gas);
     break;
 }


 // [45]-[54] 这 10 列上用来指定低温
 // thermoLowTemp         .{10}
 // [55]-[64] 这 10 列上指定高温
 // thermoHighTemp        .{10}
 // [65]-[72] 这 8 列上指定常温
 // thermoCommonTemp      .{8}
 currentLowT = stringToScalar(temperaturesString(0, 10)); // 这里是相对列数
 currentHighT = stringToScalar(temperaturesString(10, 10));
 currentCommonT = stringToScalar(temperaturesString(20, 8));

 if (currentCommonT < SMALL)
 {
     currentCommonT = allCommonT;
 }

 // 每个数字限制15列
 // thermoCoeff           .{15}
 highCpCoeffs[0] = stringToScalar(thermoCoeffString(0, 15));
 highCpCoeffs[1] = stringToScalar(thermoCoeffString(15, 15));
 highCpCoeffs[2] = stringToScalar(thermoCoeffString(30, 15));
 highCpCoeffs[3] = stringToScalar(thermoCoeffString(45, 15));
 highCpCoeffs[4] = stringToScalar(thermoCoeffString(60, 15));


 highCpCoeffs[5] = stringToScalar(thermoCoeffString(0, 15));
 highCpCoeffs[6] = stringToScalar(thermoCoeffString(15, 15));

 lowCpCoeffs[0] = stringToScalar(thermoCoeffString(30, 15));
 lowCpCoeffs[1] = stringToScalar(thermoCoeffString(45, 15));
 lowCpCoeffs[2] = stringToScalar(thermoCoeffString(60, 15));


 lowCpCoeffs[3] = stringToScalar(thermoCoeffString(0, 15));
 lowCpCoeffs[4] = stringToScalar(thermoCoeffString(15, 15));
 lowCpCoeffs[5] = stringToScalar(thermoCoeffString(30, 15));
 lowCpCoeffs[6] = stringToScalar(thermoCoeffString(45, 15));

 speciesThermo_.insert
 (
     currentSpecieName,
     new gasHThermoPhysics
     (
         janafThermo<perfectGas<specie>>
         (
             specie
             (
                 currentSpecieName,
                 1.0,
                 molecularWeight(currentSpecieComposition)
             ),
             currentLowT,
             currentHighT,
             currentCommonT,
             highCpCoeffs,
             lowCpCoeffs,
             true
         ),
         transportDict_.subDict(currentSpecieName)
     )
 );

 // Read reactions

构造 multiComponentMixture 的部分，调用 chemistryReader::speciesThermo() 这个纯虚函数。
它定义于 chemkinReader，返回 HashPtrTable<gasHThermoPhysics> speciesThermo_;。
构造 PtrList<Reaction<ThermoType>>，通过复制 chemkinReader::reactions_。
构造 speciesComposition_，通过复制 chemkinReader::speciesComposition_。

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

ReactionList<gasHThermoPhysics> reactions_;

并且定义了接口 reactions()。

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

第二个问题：

在 StandardChemistryModel 中申明了反应的 list：

const ReactionList<ThermoType> reactions_;

StandardChemistryModel 的构造函数中：

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

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

reactingMixture 中实际上包括了从机理文件读取的反应的 list，这里 reactions_ 获取了反应 list 的引用。

StandardChemistryModel::calculateRR 中使用了 Reaction<ReactionThermo>::omega

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

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 等纯虚函数，这样就能解释的通了。

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() 就相当于在这里求解化学反应！

因此，实际调用的是派生类 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(c_, Ti, pi, 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;
}

在其中，调用了 5 参数的 solve 函数。5 参数的 solve 函数在继承自 StandardChemistryModel 的 chemistrySolver 中被申明为纯虚函数，实际上在 EulerImplicit 和 ode 中重新定义。所以这里的 reaction->correct() 就相当于在这里求解化学反应！

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

template<class ChemistryModel>
void Foam::ode<ChemistryModel>::solve
(
    scalarField& c, // 通过该函数后，被更新
    scalar& T, // 通过该函数后，被更新
    scalar& p, // 通过该函数后，被更新
    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_, subDeltaT);

    for (int i=0; i<nSpecie; i++)
    {
        c[i] = max(0.0, cTp_[i]); // 计算得到新的 c
    }
    T = cTp_[nSpecie]; // 计算得到新的 T
    p = cTp_[nSpecie+1]; // 计算得到新的 p
}

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

ode 的 time split 体现在这里！！！
//这里的 x 是指时间步进
void Foam::ODESolver::solve
(
    const scalar xStart, //实参=0
    const scalar xEnd, // 实参=deltaT
    scalarField& y, // 实参=cTp_，通过该函数后，被更新
    scalar& dxTry // 化学时间步长
) const
{
    stepState step(dxTry);
    scalar x = xStart; // x 从 0 开始

    for (label nStep=0; nStep<maxSteps_; nStep++)
    // 这里是最大循环次数，但是中间可能退出循环
    // maxSteps_(dict.lookupOrDefault<scalar>("maxSteps", 10000))
    // maxSteps 设置于 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, step); // 每调用一次，x 步进或步退

        // Check if reached xEnd
        //即 x >= 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, stepState& step)，这个函数在 ODESolver和它的派生类 seulex 中都定义了，那么这里应该是调用的是派生类 seulex 的。

void Foam::seulex::solve
(
    scalar& x, // 时间
    scalarField& y, // 实参=cTp_
    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_)
    {
        // 使用 x y 更新 dfdx_ 和 dfdy_， dfdx_ 是 dcdt，dfdy_是雅克比矩阵
        odes_.jacobian(x, y, 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++)
        {
            // seul 中用到了上边计算得到的 dfdy_，dfdy_ 就是 jacobian 雅克比矩阵
            bool success = seul(x, y0_, 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, dfdx_, dfdy_); // 使用 x y 更新 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, dfdx_, dfdy_) ， seulex::seul 和 seulex::extrapolate。

seulex::seul 中调用了 ODESystem::derivatives。

但是 ODESystem 中只定义了三个纯虚函数 nEqns， derivatives， jacobian。
实际定义都在 StandardChemistryModel 中，因为 ODESystem 是 StandardChemistryModel 的基类之一。

derivatives 和 jacobian 这两个函数很重要，代码如下：

template<class ReactionThermo, class ThermoType>
void Foam::StandardChemistryModel<ReactionThermo, ThermoType>::derivatives
(
    const scalar time, //这个参数不起任何作用（很合理，因为是瞬态项，和时间无关）
    const scalarField& c, //cTp
    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(c_, T, p, dcdt);
    // 给定 c_ T p，更新 dcdt，就是瞬时的化学反应速率
    // 事实上，还是调用的 Reaction::omega，定义于 Reaction.C

    // 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;

    // dT/dt
    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,  // 实参=cTp_
    scalarField& dcdt, // 它也是，但是这个量计算出来，没有任何地方调用！
    scalarSquareMatrix& J //雅克比矩阵 ddc(dcdt)
) 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;
        // 下边两个函数用于计算 Jacobian-J 和 derivative-dcdt，定义于 Reaction.C
        R.dwdc(p, T, c_, J, dcdt, omegaI, kfwd, kbwd, false, dummy);
        R.dwdT(p, T, c_, omegaI, kfwd, kbwd, J, false, dummy, nSpecie_);
    }

    // 以下计算 dT'/dYi
    // 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)
    }


    // 以下计算 dT'/dT
    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 功能：
给定当前的 c T p，计算得到 deltaT 之后新的浓度 c。
然后据此求解组分方程的源项：

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

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

6 阶 RK

$y_{n+1}=y_n+h(c_1k_1+c_2k_2+c_3k_3+c_4k_4+c_5k_5+c_6k_6)$

$k_1=f(t_n,y_n)=dydt$

$k_2=f(t_n+a_2h,y_n+hb_{21}k_1)$

$k_3=f(t_n+a_3h,y_n+hb_{31}k_1+hb_{32}k_2)$

$k_4=f(t_n+a_4h,y_n+hb_{41}k_1+hb_{42}k_2+hb_{43}k_3)$

$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)$

$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_`	$k_2$
`ak3_`	$k_3$
`ak4_`	$k_4$
`ak5_`	$k_5$
`ak6_`	$k_6$
`y`	$y_n$
`yout`	$y_{n+1}$
`dydx`	$k_1$
第一个 `yTemp_`	$y_n+hb_{21}k_1$
第二个 `yTemp_`	$y_n+hb_{31}k_1+hb_{32}k_2$
第三个 `yTemp_`	$y_n+hb_{41}k_1+hb_{42}k_2+hb_{43}k_3$
第四个 `yTemp_`	$y_n+hb_{51}k_1+hb_{52}k_2+hb_{53}k_3+hb_{54}k_4$
第五个 `yTemp_`	$y_n+hb_{61}k_1+hb_{62}k_2+hb_{63}k_3+hb_{64}k_4+hb_{65}k_5$

代码中返回的是： $y_{n+1}=y_n+h(c_1k_1+c_3k_3+c_4k_4+c_6k_6)$ ，只需要保证 $\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_{n+1}=y_n+hf(y_{n+})$

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

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

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

$(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$

$g_i=\sum_{j=1}^{i-1}\gamma_{ij}k_j+\gamma k_i$

$(1/\gamma h -f')\cdot g_1=f(y_0)$

$(1/\gamma h -f')\cdot g_2=f(y_0+a_{21}g_1)+c_{21}g_1/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/\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_`	$f'$
`dfdx_`	$\dfrac{\partial f}{\partial x}$
`g1_`	$g_1$
`g2_`	$g_2$
`g3_`	$g_3$
`g4_`	$g_4$
`c21`	$c_{21}$
`c31`	$c_{31}$
`c32`	$c_{32}$
`c41`	$c_{41}$
`c42`	$c_{42}$
`c43`	$c_{43}$
`a21`	$a_{21}$
`a31`	$a_{31}$
`a32`	$a_{32}$
`b1`	$b1$
`b2`	$b2$
`b3`	$b3$
`b4`	$b4$

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 求解理论

OpenFOAM 代码对应的

derivatives, size= nSpecie_ + 2

$f = {\Phi_i'}=\frac{\partial \Phi}{\partial t} = \left \lbrace \frac{\partial c_1}{\partial t}, \frac{\partial c_2}{\partial t}, \dotsc, \frac{\partial c_N}{\partial t}, \frac{\partial T}{\partial t}, \frac{\partial p}{\partial t} \right \rbrace^{\text{T}}$

$\frac{\partial c_k}{\partial t} = \dot{\omega}_k \quad k = 1, \dotsc, N$

$\frac{\partial T}{\partial t} = \frac{-1}{\rho c_p} \sum_{k=1}^{N} h_k \frac{\partial c_k}{\partial t}$

$\frac{\partial p}{\partial t} = 0$

Jacobian, size= (nSpecie_+1)*(nSpecie+1)

$\mathcal{J}_{i,j} = \frac{\partial f_i}{\partial \Phi_j} = \frac{\partial \dot{\Phi_i}}{\partial \Phi_j}$

展开为：

$\left[ \begin{array}{cccc} \dfrac{\partial {c_1'}}{\partial c_1} & \dfrac{\partial {c_1'}}{\partial c_2} & \ldots &\ldots& \dfrac{\partial {c_1'}}{\partial c_N} & \dfrac{\partial {c_1'}}{\partial T}\\ \dfrac{\partial {c_2'}}{\partial c_1} & \dfrac{\partial {c_2'}}{\partial c_2} & \ldots &\ldots& \dfrac{\partial {c_2'}}{\partial c_N} & \dfrac{\partial {c_2'}}{\partial T}\\ \vdots & \vdots & \ddots && \vdots & \vdots \\ \vdots & \vdots & &\ddots & \vdots & \vdots \\ \dfrac{\partial c_N'}{\partial c_1} & \dfrac{\partial c_N'}{\partial c_2} & \ldots&\ldots & \dfrac{\partial c_N'}{\partial c_N} & \dfrac{\partial {c_N'}}{\partial T}\\ \dfrac{\partial {T'}}{\partial c_1} & \dfrac{\partial {T'}}{\partial c_2} & \ldots&\ldots & \dfrac{\partial {T'}}{\partial c_N} & \dfrac{\partial {T'}}{\partial T} \end{array} \right]$

pyJac 中描述的：

官方文档：http://slackha.github.io/pyJac/overview.html

待求解的量有：
$\Phi = \left \lbrace T, Y_1, Y_2, \dotsc, Y_{N_{\text{sp}} - 1} \right \rbrace^{\text{T}}$

$Y_{N_{\text{sp}}} = 1 - \sum_{k=1}^{N_{\text{sp}} - 1} Y_k$

derivatives, size= nSpecie_

$f = {\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}}$

$\frac{\partial T}{\partial t} = \frac{-1}{\rho c_p} \sum_{k=1}^{N_{\text{sp}}} h_k W_k \dot{\omega}_k$

$\frac{\partial Y_k}{\partial t} = \frac{1}{\rho} W_k \dot{\omega}_k \quad k = 1, \dotsc, N_{\text{sp}} - 1$

Jacobian, size= nSpecie_*nSpecie_

$\mathcal{J}_{i,j} = \frac{\partial f_i}{\partial \Phi_j} = \frac{\partial \dot{\Phi_i}}{\partial \Phi_j}$

展开为：

$\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），公式是 $\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。

OF-8

UML 类图

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

OF-8 UML

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

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

ReversibleReaction 的模板是 Reaction，gasHThermoPhysics 和 FallOffReactionRate。而 FallOffReactionRate 的模板是 ArrheniusReactionRate 和 TroeFallOffFunction。
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()
        );
    }
}

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 分别为 word 和 dictionaryConstructorPtr 的 HashTable 重命名为 argNames##ConstructorTable，即 dictionaryConstructorTable，并创建了一个以它为类型的指针 argNames##ConstructorTablePtr_，即 dictionaryConstructorTablePtr_。
申明一个类：add##argNames##ConstructorToTable，即 adddictionaryConstructorToTable，它需要模板 baseType##Type，其中的 baseType 是 CombustionModel。
类中定义了 New 函数，argList 是其形参列表，parList 是返回的 baseType##Type 的构造函数的形参列表。注意，这里调用了基类的构造函数！！！
但是如果我们在派生类中创建这个 add 类的对象的话，它的模板则是派生类，就是说 New 函数的返回值是指向派生类对象基类指针。
在 add##argNames##ConstructorToTable 的构造函数中，向 dictionaryConstructorTablePtr_ 中 insert 了一组 key and value：
key 是 typeName，vaule 是返回 基类 autoPtr 的 New 函数。

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)

这里的 baseType 是 CombustionModel##psiReactionThermo，argNames 是 dictionary。
包括三个宏函数，第一个函数定义了一个类型为 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 就是 CombustionModel，reactionThermo 是它的模板 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::typeName 和 New 函数。
说起这个 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 的对象