燃烧反应中的能量方程

不同形式的能量方程,关于反应放热源项,绝对焓输运返方程无此源项,显焓输运方程的源项是生成焓乘以反应速率,温度输运方程的源项是绝对焓乘以反应速率。

绝对焓的输运方程:

ρht+ρuihxi=DpDt+J+τijuixj\frac{\partial \rho h}{\partial t}+\frac{\rho u_i h}{\partial x_i}=\frac{D p}{D t}+\nabla\cdot J+\tau_{ij}\frac{\partial u_i}{\partial x_j}

其中J=kTρhsVkYkρΔhf,k0VkYk=JsρΔhf,k0VkYkJ=k\nabla T-\rho\sum h_sV_kY_k-\rho\sum\Delta h_{f,k}^0 V_kY_k=J_s-\rho\sum\Delta h_{f,k}^0 V_kY_k
Js=kTρhsVkYkJ_s=k\nabla T-\rho\sum h_sV_kY_kJJ
显焓hs=hk=1NΔhf,k0Ykh_s=h-\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}} 代入绝对焓的输运方程,得到:

ρhs+k=1NΔhf,k0Ykt+ρui(hs+k=1NΔhf,k0Yk)xi=DpDt+J+τijuixj\partial\rho\frac{ h_s+\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}}{\partial t}+\partial\frac{\rho u_i (h_s+\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}})}{\partial x_i} =\frac{D p}{ D t}+\nabla\cdot J+\tau_{ij}\frac{\partial u_i}{\partial x_j}

移项,得到:

ρhst+ρuihsxi=DpDt+J+τijuixjρk=1NΔhf,k0Yktρuik=1NΔhf,k0Ykxi=DpDt+Js+τijuixjρk=1NΔhf,k0Yktρuik=1NΔhf,k0Ykxi(ρk=1NΔhf,k0Vk,iYk)xi=DpDt+Js+τijuixjΔhf,k0[ρYkt+ρ(ui+Vk,i)Ykxi]=DpDt+Js+τijuixjΔhf,k0ω˙k=DpDt+Js+τijuixj+ω˙T\begin{array}{ll} &\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i}\\ &= \frac{D p}{ D t}+\nabla\cdot J + \tau_{ij}\frac{\partial u_i}{\partial x_j} - \partial\rho\frac{\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}}{\partial t} -\partial\frac{\rho u_i\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}}{\partial x_i}\\ &= \frac{D p}{ D t}+\nabla\cdot J_s+\tau_{ij}\frac{\partial u_i}{\partial x_j}-\partial\rho\frac{\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}}{\partial t}-\partial\frac{\rho u_i\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}}{\partial x_i} -\partial\frac{(\rho\sum_{k = 1}^N{\Delta h_{f,k}^0V_{k,i}{Y_k}})}{\partial x_i}\\ &= \frac{D p}{ D t}+\nabla\cdot J_s+\tau_{ij}\frac{\partial u_i}{\partial x_j}-\sum\Delta h_{f,k}^0[\frac{\partial\rho Y_k}{\partial t}+\partial\frac{\rho(u_i+V_{k,i})Y_k}{\partial x_i}]\\ &= \frac{D p}{ D t}+\nabla\cdot J_s+\tau_{ij}\frac{\partial u_i}{\partial x_j}-\sum\Delta h_{f,k}^0\dot\omega_k\\ &=\frac{D p}{ D t}+\nabla\cdot J_s+\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T \end{array}

其中:

DpDt=pt+??\frac{D p}{D t}=\frac{\partial p }{\partial t}+??

上边推导的最后一步,是由组分的输运方程得到的:

ρYkt+ρ(ui+Vk,i)Ykxi=ω˙k\frac{\partial\rho Y_k}{\partial t}+\partial\frac{\rho(u_i+V_{k,i})Y_k}{\partial x_i}=\dot\omega_k

ω˙T=Δhf,k0ω˙k\dot\omega_T=-\sum\Delta h_{f,k}^0\dot\omega_k 这一项是燃烧热。

Δhf,k0\Delta h_{f,k}^0 是组分 k 的化学焓(chemical),也叫生成焓(formation)。

viscous heating source term:τijuixj\tau_{ij}\frac{\partial u_i}{\partial x_j} 这一项是粘性热。

显焓的输运方程:

ρhst+ρuihsxi=DpDt+Js+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} =\frac{D p}{ D t}+\nabla\cdot J_s+\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

Js=kTρk=1NhsVkYkJ_s=k\nabla T-\rho\sum_{k = 1}^N h_sV_kY_k

ρhst+ρuihsxi=DpDt+(kTρk=1NhsVkYk)+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} =\frac{D p}{ D t}+\nabla\cdot (k\nabla T-\rho\sum_{k = 1}^N h_sV_kY_k)+\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

ρhst+ρuihsxi=DpDt+xi(kTxi)ρhs,kVk,iYkxi+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} = \frac{D p}{ D t} + \frac{\partial}{\partial x_i}(k\frac{\partial T}{\partial x_i})-\partial\frac{\rho\sum h_{s,k}V_{k,i}Y_k}{\partial x_i} +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

ρhs,kVk,iYkxi\partial\frac{\rho\sum h_{s,k}V_{k,i}Y_k}{\partial x_i} 经常被视为0
所以显焓的输运方程:

ρhst+ρuihsxi=DpDt+xi(κTxi)+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} = \frac{D p}{ Dt} + \frac{\partial}{\partial x_i}(\kappa\frac{\partial T}{\partial x_i}) +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

κ\kappa 就是λ\lambda,导热系数。
由于温度和焓的关系:α=λCp\alpha=\frac{\lambda}{C_p}, 代入上式,得:

ρhst+ρuihsxi=DpDt+xi(αhsxi)+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} = \frac{D p}{ D t} + \frac{\partial}{\partial x_i}(\alpha\frac{\partial h_s}{\partial x_i}) +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

温度的输运方程:

显焓的输运方程:

ρdhsdt=ρhst+ρuihsxi=DpDt+xi(kTxi)ρhs,kVk,iYkxi+τijuixj+ω˙T\rho\frac{d h_s}{dt}=\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} = \frac{D p}{ D t} + \frac{\partial}{\partial x_i}(k\frac{\partial T}{\partial x_i})-\partial\frac{\rho\sum h_{s,k}V_{k,i}Y_k}{\partial x_i} +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

hs=k=1NYkhs,kh_s=\sum_{k = 1}^N Y_k{h_{s,k}}

代入显焓输运方程得到:

ρddtk=1NYkhs,k=DpDt+xi(kTxi)ρk=1Nhs,kVk,iYkxi+τijuixj+ω˙T\rho\frac{d }{dt}\sum_{k = 1}^N Y_k{h_{s,k}} = \frac{D p}{ D t} + \frac{\partial}{\partial x_i}(k\frac{\partial T}{\partial x_i})-\partial\frac{\rho\sum_{k = 1}^N h_{s,k}V_{k,i}Y_k}{\partial x_i} +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

方程第一项:

ρddtk=1N(Ykhs,k)=ρk=1Nddt(hs,kYk)=ρk=1N(Ykdhs,kdt)+ρk=1Nhs,kdYkdt\rho\frac{d }{dt}\sum_{k = 1}^N (Y_k{h_{s,k}})=\rho\sum_{k = 1}^N \frac{d}{dt}(h_{s,k}Y_k)=\rho\sum_{k = 1}^N (Y_k\frac{dh_{s,k}}{dt})+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}

hs,k=T0TCp,kdTh_{s,k}=\int_{T0}^TC_{p,k}dT 代入上式,得到:

=ρk=1NYkddtT0TCp,kdT+ρk=1Nhs,kdYkdt=ρdT0Tk=1NYkCp,kdTdt+ρk=1Nhs,kdYkdt=\rho\sum_{k = 1}^N Y_k\frac{d}{dt}\int_{T0}^TC_{p,k}dT+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}=\rho\frac{d\int_{T0}^T\sum_{k = 1}^N Y_kC_{p,k}dT}{dt}+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}

又有Cp=YkCp,kC_p=\sum Y_k{C_{p,k}}, 因此:

=ρdT0TCpdTdt+ρk=1Nhs,kdYkdt=ρCpdTdt+ρk=1Nhs,kdYkdt=\rho\frac{d\int_{T0}^TC_pdT}{dt}+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}=\rho C_p\frac{dT}{dt}+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}

接下来呢?

ρCpdTdt+ρk=1Nhs,kdYkdt=pt+xi(kTxi)ρk=1Nhs,kVk,iYkxi+τijuixj+ω˙T\rho C_p\frac{dT}{dt}+\rho\sum_{k = 1}^N h_{s,k}\frac{dY_k}{dt}=\frac{\partial p}{ \partial t} + \frac{\partial}{\partial x_i}(k\frac{\partial T}{\partial x_i})-\partial\frac{\rho\sum_{k = 1}^N h_{s,k}V_{k,i}Y_k}{\partial x_i} +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

把组分输运方程:ρdYkdt=xiρVk,iYk+ω˙k\rho\frac{dY_k}{dt}=-\frac{\partial}{\partial x_i}\rho V_{k,i}Y_k+\dot\omega_k 代入上式左边第二项,得到:

ρCpdTdt+k=1Nhs,kω˙k=pt+xi(kTxi)+τijuixj+ω˙T\rho C_p\frac{dT}{dt}+\sum_{k = 1}^N h_{s,k}\dot\omega_k=\frac{\partial p}{ \partial t} + \frac{\partial}{\partial x_i}(k\frac{\partial T}{\partial x_i}) +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega_T

移项,得:

ρCpdTdt=DpDt+xi(κTxi)+τijuixj+ω˙T\rho C_p\frac{dT}{dt}=\frac{D p}{ D t}+\frac{\partial}{\partial x_i}(\kappa\frac{\partial T}{\partial x_i}) +\tau_{ij}\frac{\partial u_i}{\partial x_j}+\dot\omega'_T

ω˙T=ω˙Tk=1Nhs,kω˙k=k=1NΔhf,kω˙kk=1Nhs,kω˙k=k=1Nhkω˙k\dot\omega^\prime_T =\dot\omega_T-\sum_{k = 1}^N h_{s,k}\dot\omega_k =-\sum_{k = 1}^N\Delta h_{f,k}\dot\omega_k-\sum_{k = 1}^N h_{s,k}\dot\omega_k =-\sum_{k = 1}^N h_k\dot \omega_k

ω˙T\dot\omega_Tω˙T\dot\omega^\prime_T 的关系,来自 Poinsot 的书:
These two terms are both called heat release so that different authors use the same term for different quantities. They differ by a small amount due to the contribution of sensible enthalpy termshs,kh_{s,k} . Section 1.2.2 shows that they are equal when the heat capacitiesCp,kC_{p,k} are supposed equal for all species.

ω˙T=k=1Nhkω˙k=k=1Nhs,kω˙kk=1NΔhf,k0ω˙k=hsk=1Nω˙kk=1NΔhf,k0ω˙k=ω˙T\dot\omega^\prime_T=-\sum_{k = 1}^N h_k\dot\omega_k=-\sum_{k = 1}^N h_{s,k}\dot\omega_k-\sum_{k = 1}^N\Delta h_{f,k}^0\dot\omega_k=-h_s\sum_{k = 1}^N\dot\omega_k-\sum_{k = 1}^N \Delta h_{f,k}^0\dot\omega_k=\dot\omega_T\\

becausehs,k=T0TCp,kdTh_{s,k}=\int_{T_0}^T C_{p,k}dT.

ifCp,kC_{p,k} are supposed equal for all species, then so dohs,kh_{s,k}, thenk=1Nhs,kω˙k=hsk=1Nω˙k\sum_{k = 1}^N h_{s,k}\dot\omega_k=h_s\sum_{k = 1}^N \dot\omega_k andk=1Nω˙k=0\sum_{k = 1}^N\dot\omega_k =0, thereforehsk=1Nω˙k =0h_s\sum_{k = 1}^N \dot\omega_k\ =0

reactingFoam 中的能量方程

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)
==
Qdot
+ fvOptions(rho, he)
);

其中 Qdot = reaction->Qdot();,然后在燃烧模型中:Qdot.ref() = this->chemistryPtr_->Qdot();

template<class ReactionThermo, class ThermoType>
Foam::tmp<Foam::volScalarField>
Foam::StandardChemistryModel<ReactionThermo, ThermoType>::Qdot() const
{
tmp<volScalarField> tQdot
(
new volScalarField
(
IOobject
(
......
),
this->mesh_,
dimensionedScalar("zero", dimEnergy/dimVolume/dimTime, 0)
)
);

if (this->chemistry_)
{
scalarField& Qdot = tQdot.ref();

forAll(Y_, i)
{
forAll(Qdot, celli)
{
const scalar hi = specieThermo_[i].Hc();
Qdot[celli] -= hi*RR_[i][celli];
}
}
}

return tQdot;
}

动能输运方程:

ρKt+ρuiKxi=ujσijxi\frac{\partial \rho K}{\partial t}+\partial\frac{\rho u_i K}{\partial x_i}=u_j\frac{\partial \sigma_{ij}}{\partial x_i}

动能K=12ukukK=\dfrac{1}{2}u_ku_k

显焓+动能的输运方程推导如下,动能方程右边这一项可以表示为

ujσijxi=ujτijxiuipxiu_j\dfrac{\partial \sigma_{ij}}{\partial x_i}=u_j\dfrac{\partial \tau_{ij}}{\partial x_i}-u_i\dfrac{\partial p}{\partial x_i}

显焓方程中的DpDt=pt+uipxi\dfrac{D p}{D t}=\dfrac{\partial p}{\partial t}+u_i\dfrac{\partial p}{x_i},这里的uipxiu_i\dfrac{\partial p}{x_i} 一正一负抵消了。

因此显焓+动能的输运方程:

ρhst+ρuihsxi+ρKt+ρuiKxi=pt+xi(αhsxi)+τijuixj+ω˙T\frac{\partial \rho h_s}{\partial t}+\partial\frac{\rho u_i h_s}{\partial x_i} +\frac{\partial \rho K}{\partial t}+\partial\frac{\rho u_i K}{\partial x_i} = \frac{\partial p}{ \partial t} + \frac{\partial}{\partial x_i}(\alpha\frac{\partial h_s}{\partial x_i}) +\frac{\partial \tau_{ij} u_i}{\partial x_j}+\dot\omega_T

OF 代码中被忽略的部分是τijuixj\frac{\partial \tau_{ij} u_i}{\partial x_j}

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