燃烧反应中的能量方程

2018-03-19| combustion|字数总计: 3k|阅读时长: 17 分钟

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

绝对焓的输运方程：

$\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=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$
即 $J_s=k\nabla T-\rho\sum h_sV_kY_k$
显焓 $h_s=h-\sum_{k = 1}^N {\Delta h_{f,k}^0{Y_k}}$ 代入绝对焓的输运方程，得到：

$\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}$

移项，得到：

$\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}$

其中：

$\dfrac{D p}{D t}=\dfrac{\partial p}{\partial t}+u_i\dfrac{\partial p}{x_i}$

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

$\frac{\partial\rho Y_k}{\partial t}+\partial\frac{\rho(u_i+V_{k,i})Y_k}{\partial x_i}=\dot\omega_k$

$\dot\omega_T=-\sum\Delta h_{f,k}^0\dot\omega_k$ 这一项是燃烧热。

$\Delta h_{f,k}^0$ 是组分 k 的化学焓（chemical），也叫生成焓（formation）。

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

显焓的输运方程：

$\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$

$J_s=k\nabla T-\rho\sum_{k = 1}^N h_sV_kY_k$

$\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$

$\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$

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

$\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$ ，导热系数。
由于温度和焓的关系： $\alpha=\frac{\lambda}{C_p}$ ，代入上式，得：

$\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$

温度的输运方程:

显焓的输运方程：

$\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$

$h_s=\sum_{k = 1}^N Y_k{h_{s,k}}$

代入显焓输运方程得到：

$\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$

方程第一项：

$\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}$

将 $h_{s,k}=\int_{T0}^TC_{p,k}dT$ 代入上式，得到：

$=\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}$

又有 $C_p=\sum Y_k{C_{p,k}}$ ，因此：

$=\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}$

接下来呢？

$\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$

把组分输运方程： $\rho\frac{dY_k}{dt}=-\frac{\partial}{\partial x_i}\rho V_{k,i}Y_k+\dot\omega_k$ 代入上式左边第二项，得到：

$\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$

移项，得：

$\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$

$\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$

$\dot\omega_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 terms $h_{s,k}$ . Section 1.2.2 shows that they are equal when the heat capacities $C_{p,k}$ are supposed equal for all species.

$\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\\$

because $h_{s,k}=\int_{T_0}^T C_{p,k}dT$ .

if $C_{p,k}$ are supposed equal for all species, then so do $h_{s,k}$ , then $\sum_{k = 1}^N h_{s,k}\dot\omega_k=h_s\sum_{k = 1}^N \dot\omega_k$ and $\sum_{k = 1}^N\dot\omega_k =0$ , therefore $h_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;
}

动能输运方程：

$\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=\dfrac{1}{2}u_ku_k$ 。

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

$u_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}$

显焓方程中的 $\dfrac{D p}{D t}=\dfrac{\partial p}{\partial t}+u_i\dfrac{\partial p}{x_i}$ ，这里的 $u_i\dfrac{\partial p}{x_i}$ 一正一负抵消了。

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

$\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 代码中被忽略的部分是 $\frac{\partial \tau_{ij} u_i}{\partial x_j}$ 。