OpenFOAM 是如何从焓得到温度的

在 solver 中，求解完 h（单位是 J/Kg）的方程之后，调用 thermo.correct(); 该函数位于 hePsiThermo.C的 correct() 函数中，调用了 calculate()，calculate() 大致内容如下：

const scalarField& hCells = this->he_;
const scalarField& pCells = this->p_;
scalarField& TCells = this->T_.primitiveFieldRef();

    TCells[celli] = mixture_.THE
    (
        hCells[celli],
        pCells[celli],
        TCells[celli]
    );

THE(h,p,T) 函数在 thermoI.H

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

THE(*this, he, p, T0) 函数在 sensibleEnthalpy.H

scalar THE
(
    const Thermo& thermo,
    const scalar h,
    const scalar p,
    const scalar T0
) const
{
    return thermo.THs(h, p, T0);
}

THs(h, p, T0) 函数在 thermoI.H

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

inline Foam::scalar Foam::species::thermo<Thermo, Type>::T
(
    scalar f,
    scalar p,
    scalar T0,
    scalar (thermo<Thermo, Type>::*F)(const scalar, const scalar) const,
    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;
        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;
}

上边的 F 是 Hs 函数（单位是 J/Kg），dFdT 是 Cp 函数，f 是输运方程求解得到的每个网格的 hs（单位是 J/Kg），Test 是 Told。

为了通过焓求温度，利用牛顿迭代法求解 Hs(p_old,T) - hs =0 ，令 f(x) = Hs(p_old,T) - hs，则：
T_new = T_old - [ Hs(p_old,T_old) - hs（输运方程得到的）]/Cp(p_old,T_old)

由 Ha 求 T 呢？也是
T_new = T_old - [ Ha(p_old,T_old) - ha（输运方程得到的）]/Cp(p_old,T_old)

上边用到了 Hs(p,T) 和 Cp(p,T) 函数，位于 janafThermoI.H

inline Foam::scalar Foam::janafThermo<EquationOfState>::Ha
(
    const scalar p,
    const scalar T
) const
{
    const coeffArray& a = coeffs(T);
    return
    (
        ((((a[4]/5.0*T + a[3]/4.0)*T + a[2]/3.0)*T + a[1]/2.0)*T + a[0])*T
      + a[5]
    ) + EquationOfState::H(p, T);
}

inline Foam::scalar Foam::janafThermo<EquationOfState>::Hs
(
    const scalar p,
    const scalar T
) const
{
    return Ha(p, T) - Hc();
}

inline Foam::scalar Foam::janafThermo<EquationOfState>::Hc() const
{
    const coeffArray& a = lowCpCoeffs_;
    return
    (
        (
            (((a[4]/5.0*Tstd + a[3]/4.0)*Tstd + a[2]/3.0)*Tstd + a[1]/2.0)*Tstd
          + a[0]
        )*Tstd + a[5]
    );
}


inline Foam::scalar Foam::janafThermo<EquationOfState>::Cp
(
    const scalar p,
    const scalar T
) const
{
    const coeffArray& a = coeffs(T);
    return
        ((((a[4]*T + a[3])*T + a[2])*T + a[1])*T + a[0])
      + EquationOfState::Cp(p, T);
}

其中 EquationOfState::H(p, T) EquationOfState::Cp(p, T) 来自EOS，比如： [perfectGas(https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/thermophysicalModels/specie/equationOfState/perfectGas/perfectGasI.H) 中这两个函数都返回 0。

上边 [janafThermoI.H](https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/thermophysicalModels/specie/thermo/janaf/janafThermoI.H) 还用到了机理 therm.dat 文件的系数，但是是经过转换的系数：

```c++
inline const typename Foam::janafThermo<EquationOfState>::coeffArray&
Foam::janafThermo<EquationOfState>::coeffs
(
    const scalar T
) const
{
    if (T < Tcommon_)
    {
        return lowCpCoeffs_;
    }
    else
    {
        return highCpCoeffs_;
    }
}

// Convert coefficients to mass-basis
for (label coefLabel=0; coefLabel<nCoeffs_; coefLabel++)
{
    highCpCoeffs_[coefLabel] *= this->R();
    lowCpCoeffs_[coefLabel] *= this->R();
}

意思就是，系数从 therm.dat 读取进来，先乘一个 R()， 再拿来使用的。

还有一个 HE 函数，忘了有什么用了。。。位于 sensibleEnthalpy.H

scalar HE
(
    const Thermo& thermo,
    const scalar p,
    const scalar T
) const
{
    return thermo.Hs(p, T);
}

另外 thermoI.H 中：

inline Foam::scalar
Foam::species::thermo<Thermo, Type>::HE(const scalar p, const scalar T) const
{
    return Type<thermo<Thermo, Type>>::HE(*this, p, T);
}
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::he(const scalar p, const scalar T) const
{
    return this->HE(p, T)*this->W();
}

可能是一般调用的是抽象出来的 HE(p,T) 或 he(p,T) 函数，最终实际调用的是 Hs(p, T) 或 Es 或 Ha 或 Ea（根据能量形式自动选择）。
小写单位是 [J/kmol]，大写单位是 [J/kg]。输运方程中求解的量也是 [J/kg]。