OpenFOAM 成长之路 首页关于Timeline
OpenFOAM 如何更新混合物的属性
| code analysis|字数总计: 1.4k|阅读时长: 7 分钟

OpenFOAM 如何更新混合物的属性

本文代码全部在 /src/thermophysicalModels 文件夹下。

温度更新于 /basic/psiThermo/hePsiThermo.C 中：

forAll(TCells, celli)
{
    const typename MixtureType::thermoType& mixture_ =
        this->cellMixture(celli);
    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]);
}

  1. 第一步： mixture_ = this->cellMixture(celli); 先通过质量分数加权得到 mixture 的值，mixture 相当于一个特殊的 specie。

  2. 第二步：使用 THE 这个函数用焓更新温度，这个过程比较复杂，见本站另一篇博文。

  3. 第三步：计算 psi，函数见 specie/equationOfState/perfectGasI.H：return 1.0/(this->R()*T);
    其中 R() 在specie\specie 计算: return RR/molWeight_; R() 的单位是 [J/(kg K)]。这里的 molWeight_ 是 mixture 的值，W() 也能返回 molWeight_，单位是[kg/kmol]。
    RR 的位置更加隐蔽。首先给出结果：RR=1000*R=8314.47J/(kmol·k)。详细代码挖掘见本文末尾。

  4. 第四步：计算 mu，见 /specie/transport/sutherland/sutherlandTransportI.H，return As_*::sqrt(T)/(1.0 + Ts_/T);

  5. 第五步：计算 alpha[kg/m/s]，同样见 /specie/transport/sutherland/sutherlandTransportI.H，return kappa(p,T) /Cp(p,T) 这里的 kappa 就是导热系数，我们所说的 lamda。

alphah---Thermal diffusivity of enthalpy [kg/ms]    在Sutherland中计算
kappa----Thermal conductivity [W/mK]              也在Sutherland中计算
template<class Thermo>
inline Foam::scalar Foam::sutherlandTransport<Thermo>::kappa
(
    const scalar p, const scalar T
) const
{
    scalar Cv_ = this->Cv(p, T);
    return mu(p, T)*Cv_*(1.32 + 1.77*this->R()/Cv_);
}


template<class Thermo>
inline Foam::scalar Foam::sutherlandTransport<Thermo>::alphah
(
    const scalar p,
    const scalar T
) const
{
    return kappa(p, T)/this->Cp(p, T);
}

//thermI.H 中：
template<class Thermo, template<class> class Type>
inline Foam::scalar
Foam::species::thermo<Thermo, Type>::Cpv(const scalar p, const scalar T)
const
{
    return Type<thermo<Thermo, Type>>::Cpv(*this, p, T);
}

//单位：
//- Heat capacity at constant volume [J/kg/K]
//virtual tmp<volScalarField> Cv() const = 0;
//大写是[J/kg/K]，小写是[J/kmol/K]???
//Cv = cv(p, T)/this->W()
    
//cpv(或者Cpv)：当能量是焓的形式，cpv = cp，若能量是内能形式，cpv = cv

//那么 mixture 的 Cp(p, T) 又是如何计算的？
//单个组分的 Cp 计算见 [janafThermoI.H](https://cpp.openfoam.org/v3/a10297_source.html)，使用机理文件 therm.dat 中的系数计算得到。
template<class EquationOfState>
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);
}

//SpecieMixture.C 中：
template<class MixtureType>
Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
(
    const label speciei,
    const scalar p,
    const scalar T
) const
{
    return this->getLocalThermo(speciei).Cp(p, T);
}

最后，什么是 mixture 的加权平均？
见 \reactionThermo\mixtures\multiComponentMixture\multiComponentMixture.C

template<class ThermoType>
const ThermoType& Foam::multiComponentMixture<ThermoType>::cellMixture
(
    const label celli
) const
{
    mixture_ = Y_[0][celli]*speciesData_[0];
    for (label n=1; n<Y_.size(); n++)
    {
        mixture_ += Y_[n][celli]*speciesData_[n];
    }
    return mixture_;
}

mixture_ = sum( Y[i] * speciesData[i]) i 从 1 到组分个数

就是说混合物，是每个组分按照质量分数加权得到的。
speciesData_ 包括哪些数据？ PtrList<ThermoType> speciesData_;
speciesData_ 是由所有组分的元素组成的一个 list。
这里的 ThermoType 是指sutherlandTransport<species::thermo<janafThermo<perfectGas<specie>>,sensibleEnthalpy>>。
sutherlandTransport 继承自 species::thermo。
species::thermo 继承自 janafThermo + sensibleEnthalpy。
janafThermo 继承自 perfectGas。
perfectGas 继承自 specie。
此外，sensibleEnthalpy 继承自 species::thermo。
sensibleEnthalpy< species::thermo<janafThermo<perfectGas<specie>>,sensibleEnthalpy> >
species 是命名空间。
另外，它还有一个别名，见 thermoPhysicsTypes.H：

typedef
sutherlandTransport
<
    species::thermo
    <
        janafThermo
        <
            perfectGas<specie>
        >,
        sensibleEnthalpy
    >
> gasHThermoPhysics;
//即 typedef sutherlandTransport<species::thermo<janafThermo<perfectGas<specie>>,sensibleEnthalpy>> gasHThermoPhysics;

同样，mixture_ 也是这个类型，说明混合物在这里相当于一个组分，代表所有组分的一个组分。质量分数加权，就是把单个组分的所有信息都加权一遍！如下图所示的所有信息。包括：Y_ moleWeight_ name_ Tcommon_ Thigh_ Tlow_ highCpCoeffs_ lowCpCoeffs_ As_ Ts_。

  • As_ Ts_ 是由 suterlandTransport 提供的；
  • Tcommon_ Thigh_ Tlow_ highCpCoeffs_ lowCpCoeffs_ 是由 janafThermo 提供的；therm.dat 文件中，每个组分有 14（或者 15，最后一个没用）参数，依次是高温的 7 个系数和低温的 7 个系数。但并不是 highCpCoeffs_ ，还需要 *R()，即 *RR/moleWeight_ 才成为 highCpCoeffs_ 。RR=8314.47 J/(kmol·k)，moleWeight 单位是 g/mol（亦即 kg/kmol）。
  • Y_ moleWeight_ name_ 是由 Foam::specie 提供的。

template<class EquationOfState>
inline void Foam::janafThermo<EquationOfState>::operator+=
(
    const janafThermo<EquationOfState>& jt
)
{
    scalar Y1 = this->Y();

    EquationOfState::operator+=(jt);

    if (mag(this->Y()) > small)
    {
        Y1 /= this->Y();
        const scalar Y2 = jt.Y()/this->Y();

        Tlow_ = max(Tlow_, jt.Tlow_);
        Thigh_ = min(Thigh_, jt.Thigh_);

        if
        (
            janafThermo<EquationOfState>::debug
         && notEqual(Tcommon_, jt.Tcommon_)
        )
        {
            FatalErrorInFunction
                << "Tcommon " << Tcommon_ << " for "
                << (this->name().size() ? this->name() : "others")
                << " != " << jt.Tcommon_ << " for "
                << (jt.name().size() ? jt.name() : "others")
                << exit(FatalError);
        }

        for
        (
            label coefLabel=0;
            coefLabel<janafThermo<EquationOfState>::nCoeffs_;
            coefLabel++
        )
        {
            highCpCoeffs_[coefLabel] =
                Y1*highCpCoeffs_[coefLabel]
              + Y2*jt.highCpCoeffs_[coefLabel];

            lowCpCoeffs_[coefLabel] =
                Y1*lowCpCoeffs_[coefLabel]
              + Y2*jt.lowCpCoeffs_[coefLabel];
        }
    }
}

上边这段代码说明，机理文件中的thermo文件，中间那个温度一定要保持一致，不然会有问题。不过这个只有在debug模式才会提示。

RR 代码挖掘：

定义在\src\OpenFOAM\global\constants\thermodynamic\thermodynamicConstants.C
const scalar RR = 1e3*physicoChemical::R.value();

R 的定义在 \src\OpenFOAM\global\constants\physicoChemical\physicoChemicalConstants.C

defineDimensionedConstantWithDefault
(
    physicoChemical::group,
    physicoChemical::R,
    dimensionedScalar
    (
        "R",
        physicoChemical::NA*physicoChemical::k
    ),
    constantphysicoChemicalR,
    "R"
);

k 见 \etc\controlDict:

physicoChemical
{
    mu              mu [1 0 0 0 0 0 0] 1.66054e-27;
    k               k [1 2 -2 -1 0 0 0] 1.38065e-23;
}

另外这个文件还定义了 Pstd=1e5，Tstd=298.15。
NA 见 \src\OpenFOAM\global\constants\fundamental\fundamentalConstants.C

defineDimensionedConstantWithDefault
(
    physicoChemical::group,
    physicoChemical::NA,
    Foam::dimensionedScalar
    (
        "NA",
        dimensionSet(0, 0, 0, 0, -1), //Foam::dimless/Foam::dimMoles,
        6.0221417930e+23
    ),
    constantphysicoChemicalNA,
    "NA"
);

至此，R=NA*k=6.0221417930e+23*1.38065e-23=8.31447J/(mol·k)，RR=1000*R=8314.47J/(kmol·k)

文章作者: Yan Zhang
文章链接: https://openfoam.top/mixture/
版权声明: 本博客所有文章除特别声明外，均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路！
thermophysicalModels
您的肯定会给我更大动力~