OpenFOAM 如何更新混合物的属性

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.Hreturn 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.Hreturn As_*::sqrt(T)/(1.0 + Ts_/T);

  5. 第五步:计算 alpha[kg/m/s],同样见 /specie/transport/sutherland/sutherlandTransportI.Hreturn 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 成长之路
微信打赏给博主更多动力吧~