OpenFOAM中热物理量的计算（一）

对于做反应流体的Foamer，温度是重点关注的物理量之一，但是在一般的OpenFoam求解器中能量方程输运的都是he, enthalpy/Internal energy [J/kg]，比如在coalChemistryFoam中：

{
    volScalarField& he = thermo.he(); 
    // 声明he 为thermo类中的he，即显焓
    // thermo类的设置在<constant>/thermophysicalProperties中指定

    // 求解he的输运方程
    fvScalarMatrix EEqn
    (
        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
      + fvc::ddt(rho, K) + fvc::div(phi, K)
      + (
            he.name() == "e"
          ? fvc::div
            (
                fvc::absolute(phi/fvc::interpolate(rho), U),
                p,
                "div(phiv,p)"
            )
          : -dpdt
        )
      - fvm::laplacian(turbulence->alphaEff(), he)
     ==
        rho*(U&g)
      + Qdot
      + coalParcels.Sh(he)
      + limestoneParcels.Sh(he)
      + radiation->Sh(thermo, he)
      + fvOptions(rho, he)
    );

    EEqn.relax();

    fvOptions.constrain(EEqn);

    EEqn.solve();

    fvOptions.correct(he);

    thermo.correct(); // 计算热物理量
    radiation->correct(); // 计算辐射传热

    Info<< "T gas min/max   = " << min(T).value() << ", "
        << max(T).value() << endl;
}

求解完he的输运方程之后，利用thermo.correct()更新热物理量，以coalChemistryFoam算例中的hePsiThermo为例，其correct()函数调用了calculate()函数，主要计算部分为：

template<class BasicPsiThermo, class MixtureType>
void Foam::hePsiThermo<BasicPsiThermo, MixtureType>::calculate
(
    const volScalarField& p,
    volScalarField& T,
    volScalarField& he,
    volScalarField& psi,
    volScalarField& mu,
    volScalarField& alpha,
    const bool doOldTimes
)
{
    // Note: update oldTimes before current time so that if T.oldTime() is
    // created from T, it starts from the unconverted T
......
      forAll(TCells, celli)
    {
        const typename MixtureType::thermoType& mixture_ =
            this->cellMixture(celli);

        if (this->updateT())
        {
            TCells[celli] = mixture_.THE
            (
                hCells[celli],
                pCells[celli],
                TCells[celli]
            );
        }

        psiCells[celli] = mixture_.psi(pCells[celli], TCells[celli]);

        muCells[celli] = mixture_.mu(pCells[celli], TCells[celli]);
        alphaCells[celli] = mixture_.alphah(pCells[celli], TCells[celli]);
    }
......

在correct()函数中对温度、粘性系数、导热率等参数进行了更新，其中温度计算用到了
THE(const scalar he, const scalar p, const scalar T0)

// $FOAM_SRC/thermophysicalModels/specie/thermo/thermo/thermoI.H
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THE
(
    const scalar he,
    const scalar p,
    const scalar T0
) const
{
    return Type<thermo<Thermo, Type>>::THE(*this, he, p, T0);
}

根据<constant>/thermophysicalProperties中指定的焓的类型，OF会调用不同的THE函数进一步计算温度，以显焓为例：

// $FOAM_SRC/thermophysicalModels/specie/thermo/sensibleEnthalpy/sensibleEnthalpy.H
 //- Temperature from sensible enthalpy
//  given an initial temperature T0
            scalar THE
            (
                const Thermo& thermo,
                const scalar h,
                const scalar p,
                const scalar T0
            ) const
            {
                #ifdef __clang__
                // Using volatile to prevent compiler optimisations leading to
                // a sigfpe
                volatile const scalar ths = thermo.THs(h, p, T0);
                return ths;
                #else
                return thermo.THs(h, p, T0);
                #endif
            }

进一步调用THs函数

// $FOAM_SRC/thermophysicalModels/specie/thermo/thermoI.H
template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::THs
(
    const scalar hs,
    const scalar p,
    const scalar T0
) const
{
    return T
    (
        hs,
        p,
        T0,
        &thermo<Thermo, Type>::Hs,
        &thermo<Thermo, Type>::Cp,
        &thermo<Thermo, Type>::limit
    );
}

T函数也在thermoI.H中

template<class Thermo, template<class> class Type>
inline Foam::scalar Foam::species::thermo<Thermo, Type>::T
(
    scalar f,
    scalar p,
    scalar T0, // Told 上一时刻温度
    scalar (thermo<Thermo, Type>::*F)(const scalar, const scalar) const, // 求解出的显焓，Hs_new
    scalar (thermo<Thermo, Type>::*dFdT)(const scalar, const scalar)
        const,
    scalar (thermo<Thermo, Type>::*limit)(const scalar) const
) const
{
    if (T0 < 0)
    {
        FatalErrorInFunction
            << "Negative initial temperature T0: " << T0
            << abort(FatalError);
    }

    scalar Test = T0;
    scalar Tnew = T0;
    scalar Ttol = T0*tol_;
    int    iter = 0;

    do
    {
        Test = Tnew; // Told
        Tnew =
            (this->*limit)
            (Test - ((this->*F)(p, Test) - f)/(this->*dFdT)(p, Test));

        if (iter++ > maxIter_)
        {
            FatalErrorInFunction
                << "Maximum number of iterations exceeded: " << maxIter_
                << abort(FatalError);
        }

    } while (mag(Tnew - Test) > Ttol);

    return Tnew;
}

可以看出OpenFOAM采用迭代法求解出温度，即
Tnew-Told > Ttol 时，计算 Tnew = Told - (Hs_old - Hs_new)/cp
在以上计算过程中还涉及到Cp, mu等参数的更新，下次继续介绍。