OpenFOAM 如何从焓得到温度

OpenFOAM 是如何从焓得到温度的

在 solver 中,求解完 h(单位是 J/Kg)的方程之后,调用 thermo.correct(); 该函数位于 /basic/psiThermo/hePsiThermo.Ccorrect() 函数中,调用了 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::Cp(p, T) 来自 perfectGasI.H


inline Foam::scalar Foam::perfectGas<Specie>::Cp(scalar p, scalar T) const
{
return 0;
}

上边 janafThermoI.H 还用到了机理 therm.dat 文件的系数,但是是经过转换的系数:

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]。

文章作者: Yan Zhang
文章链接: https://openfoam.top/T_calculate/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 OpenFOAM 成长之路
微信打赏给博主更多动力吧~