一、积分形式的能量方程

1、能量守恒定律

运动流体除满足质量守恒与动量定律外，还须满足能量守恒定律，即热力学第一定律：对于由确定流体质点组成的流体系统，其总能量的时间变化率等于单位时间内传给系统的热量与外力对系统所做的功之和。

记系统总能量为 $E$，传入系统的热流率为 $\sum Q_h$，外力对系统做功的功率为 $\sum N$，则

$$ \frac{\mathrm{D}E}{\mathrm{D}t}=\sum Q_h+\sum N $$

采用流体吸热为正、外力对流体做功为正的约定。

记密度为 $\rho$，速度为 $\boldsymbol{V}=u\boldsymbol{i}+v\boldsymbol{j}+w\boldsymbol{k}$，并令 $V^2=\boldsymbol{V}\cdot\boldsymbol{V}$。设单位质量流体的内能为 $e$，同一流体系统在 $t$ 时刻占据的体积为 $\tau(t)$，则系统总能量为

$$ E(t)=\int_{\tau(t)}\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau $$

其中 $V^2/2$ 为单位质量流体的动能。这里将体积力的作用作为外力功计入右端，总能量中包含内能与动能。

2、热量输入与外力功率

记 $\tau$ 为在所考察时刻与流体系统重合的固定控制体，$A$ 为其封闭控制面，$\boldsymbol{n}$ 为外法向单位矢量。设 $q$ 为单位时间传给单位质量流体的热量，可包含辐射吸热及燃烧、化学反应等效应的等效热源，则体积热源对系统的供热率为

$$ \int_\tau\rho q\,\mathrm{d}\tau $$

记温度为 $T$，导热系数为 $k$。由傅里叶热传导定律，导热热流密度矢量为

$$ \boldsymbol{q}_h=-k\nabla T $$

$\boldsymbol{q}_h\cdot\boldsymbol{n}$ 表示沿外法向流出的热流密度，因此经控制面传入系统的导热量为

$$ -\oiint_A\boldsymbol{q}_h\cdot\boldsymbol{n}\,\mathrm{d}A =\oiint_A k\nabla T\cdot\boldsymbol{n}\,\mathrm{d}A $$

总供热率为

$$ \sum Q_h =\int_\tau\rho q\,\mathrm{d}\tau +\oiint_A k\nabla T\cdot\boldsymbol{n}\,\mathrm{d}A $$

记 $\boldsymbol{f}$ 为单位质量流体所受的体积力，$\boldsymbol{p}_{\boldsymbol{n}}$ 为控制面上的应力矢量。由功率等于力与速度的点积，体积力和表面力对流体做功的总功率为

$$ \sum N =\int_\tau\rho\boldsymbol{f}\cdot\boldsymbol{V}\,\mathrm{d}\tau +\oiint_A\boldsymbol{p}_{\boldsymbol{n}}\cdot\boldsymbol{V}\,\mathrm{d}A $$

3、固定控制体的能量方程

对系统总能量应用体积分的随体导数公式，得到

$$ \begin{aligned} \frac{\mathrm{D}}{\mathrm{D}t}\int_{\tau(t)}\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau &=\frac{\partial}{\partial t}\int_\tau\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau\\ &\quad+\oiint_A\rho\left(e+\frac{V^2}{2}\right)(\boldsymbol{V}\cdot\boldsymbol{n})\,\mathrm{d}A \end{aligned} $$

代入能量守恒定律，得到积分形式的能量方程

$$ \begin{aligned} &\frac{\partial}{\partial t}\int_\tau\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau +\oiint_A\rho\left(e+\frac{V^2}{2}\right)(\boldsymbol{V}\cdot\boldsymbol{n})\,\mathrm{d}A\\ &\qquad=\int_\tau\rho q\,\mathrm{d}\tau +\oiint_A k\nabla T\cdot\boldsymbol{n}\,\mathrm{d}A +\int_\tau\rho\boldsymbol{f}\cdot\boldsymbol{V}\,\mathrm{d}\tau +\oiint_A\boldsymbol{p}_{\boldsymbol{n}}\cdot\boldsymbol{V}\,\mathrm{d}A \end{aligned} $$

左端为控制体内能量的局部积累率与通过控制面的能量净流出率之和；右端依次为体积热源、导热、体积力做功和表面力做功的贡献。

积分形式用于描述有限控制体的总体能量收支，不要求控制体内所有物理量处处连续。

二、微分形式的能量方程

1、总能量方程

在物理量满足相应光滑性条件的区域内，由雷诺第二输运方程，得

$$ \frac{\mathrm{D}}{\mathrm{D}t}\int_{\tau(t)}\rho\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau =\int_\tau\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau $$

由高斯公式，导热面积分可化为

$$ \oiint_A k\nabla T\cdot\boldsymbol{n}\,\mathrm{d}A =\int_\tau\nabla\cdot(k\nabla T)\,\mathrm{d}\tau $$

记应力张量为 $\boldsymbol{P}=(p_{ij})$，第一个下标表示受力面的法向，第二个下标表示力的分量方向。应力矢量的分量满足

$$ (\boldsymbol{p}_{\boldsymbol{n}})_j=n_i p_{ij} $$

因此，表面力功率的被积函数为

$$ \boldsymbol{p}_{\boldsymbol{n}}\cdot\boldsymbol{V} =n_i p_{ij}u_j =(\boldsymbol{P}\cdot\boldsymbol{V})\cdot\boldsymbol{n} $$

其中重复指标对 $1,2,3$ 求和。再次应用高斯公式，得到

$$ \oiint_A\boldsymbol{p}_{\boldsymbol{n}}\cdot\boldsymbol{V}\,\mathrm{d}A =\int_\tau\nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V})\,\mathrm{d}\tau $$

将各项写为体积分，得到

$$ \begin{aligned} \int_\tau\rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{V^2}{2}\right)\mathrm{d}\tau =\int_\tau\left[\rho q+\nabla\cdot(k\nabla T) +\rho\boldsymbol{f}\cdot\boldsymbol{V} +\nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V})\right]\mathrm{d}\tau \end{aligned} $$

由被积函数的连续性及控制体的任意性，得到一般的微分形式能量方程

$$ \rho\frac{\mathrm{D}}{\mathrm{D}t}\left(e+\frac{V^2}{2}\right) =\rho q+\nabla\cdot(k\nabla T) +\rho\boldsymbol{f}\cdot\boldsymbol{V} +\nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V}) $$

2、动能方程

运动方程为

$$ \rho\frac{\mathrm{D}\boldsymbol{V}}{\mathrm{D}t} =\rho\boldsymbol{f}+\nabla\cdot\boldsymbol{P} $$

两端与速度 $\boldsymbol{V}$ 点乘，得到

$$ \rho\boldsymbol{V}\cdot\frac{\mathrm{D}\boldsymbol{V}}{\mathrm{D}t} =\rho\boldsymbol{f}\cdot\boldsymbol{V} +\boldsymbol{V}\cdot(\nabla\cdot\boldsymbol{P}) $$

由乘积求导公式，

$$ \begin{aligned} \boldsymbol{V}\cdot\frac{\mathrm{D}\boldsymbol{V}}{\mathrm{D}t} &=u\frac{\mathrm{D}u}{\mathrm{D}t} +v\frac{\mathrm{D}v}{\mathrm{D}t} +w\frac{\mathrm{D}w}{\mathrm{D}t}\\ &=\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{u^2+v^2+w^2}{2}\right) =\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{V^2}{2}\right) \end{aligned} $$

因此，动能方程为

$$ \rho\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{V^2}{2}\right) =\rho\boldsymbol{f}\cdot\boldsymbol{V} +\boldsymbol{V}\cdot(\nabla\cdot\boldsymbol{P}) $$

动能方程由运动方程导出，不是独立的基本方程。

3、内能方程

总能量方程减去动能方程，体积力功率项相消，得到

$$ \rho\frac{\mathrm{D}e}{\mathrm{D}t} =\rho q+\nabla\cdot(k\nabla T) +\nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V}) -\boldsymbol{V}\cdot(\nabla\cdot\boldsymbol{P}) $$

对最后两项作指标展开，有

$$ \begin{aligned} \nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V}) &=\frac{\partial(p_{ij}u_j)}{\partial x_i}\\ &=u_j\frac{\partial p_{ij}}{\partial x_i} +p_{ij}\frac{\partial u_j}{\partial x_i} \end{aligned} $$

而

$$ \boldsymbol{V}\cdot(\nabla\cdot\boldsymbol{P}) =u_j\frac{\partial p_{ij}}{\partial x_i} $$

两式相减，得到

$$ \nabla\cdot(\boldsymbol{P}\cdot\boldsymbol{V}) -\boldsymbol{V}\cdot(\nabla\cdot\boldsymbol{P}) =p_{ij}\frac{\partial u_j}{\partial x_i} $$

对不考虑内部偶力的通常流体，应力张量对称，即 $p_{ij}=p_{ji}$。记应变率张量为 $\boldsymbol{E}=(\varepsilon_{ij})$，其中

$$ \varepsilon_{ij} =\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right) $$

利用应力的对称性，得

$$ \begin{aligned} p_{ij}\frac{\partial u_j}{\partial x_i} &=\frac{1}{2}\left(p_{ij}\frac{\partial u_j}{\partial x_i} +p_{ji}\frac{\partial u_i}{\partial x_j}\right)\\ &=p_{ij}\frac{1}{2}\left(\frac{\partial u_j}{\partial x_i} +\frac{\partial u_i}{\partial x_j}\right)\\ &=p_{ij}\varepsilon_{ij} =\boldsymbol{P}:\boldsymbol{E} \end{aligned} $$

其中冒号表示张量的双点积，即 $\boldsymbol{P}:\boldsymbol{E}=\sum_{i=1}^3\sum_{j=1}^3p_{ij}\varepsilon_{ij}$。因此，内能方程为

$$ \rho\frac{\mathrm{D}e}{\mathrm{D}t} =\rho q+\nabla\cdot(k\nabla T)+\boldsymbol{P}:\boldsymbol{E} $$

$\boldsymbol{P}:\boldsymbol{E}$ 表示单位体积流体因变形而接受的应力功率。体积力功率在总能量与动能方程中相消，故不在内能方程中直接出现。

三、牛顿流体的内能方程

1、有粘性可压缩流动

对牛顿流体，采用平均正应力的负值表示压强 $p$，应力张量写为

$$ \boldsymbol{P} =-p\boldsymbol{I}+\boldsymbol{\tau} =-p\boldsymbol{I}+2\mu\boldsymbol{E} -\frac{2}{3}\mu(\nabla\cdot\boldsymbol{V})\boldsymbol{I} $$

其中 $\boldsymbol{I}$ 为单位张量，$\boldsymbol{\tau}$ 为黏性应力张量，$\mu$ 为动力黏度。将 $p$ 同时视为热力学压强时，采用斯托克斯假设。

由单位张量与应变率张量的双点积，得

$$ \boldsymbol{I}:\boldsymbol{E} =\delta_{ij}\varepsilon_{ij} =\varepsilon_{ii} =\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y} +\frac{\partial w}{\partial z} =\nabla\cdot\boldsymbol{V} $$

因此，应力功率可分解为

$$ \begin{aligned} \boldsymbol{P}:\boldsymbol{E} &=(-p\boldsymbol{I}+\boldsymbol{\tau}):\boldsymbol{E}\\ &=-p\nabla\cdot\boldsymbol{V}+\boldsymbol{\tau}:\boldsymbol{E} \end{aligned} $$

定义黏性耗散函数为

$$ \Phi=\boldsymbol{\tau}:\boldsymbol{E} $$

代入内能方程，得到

$$ \rho\frac{\mathrm{D}e}{\mathrm{D}t} =\rho q+\nabla\cdot(k\nabla T) -p\nabla\cdot\boldsymbol{V}+\Phi $$

其中 $-p\nabla\cdot\boldsymbol{V}$ 为压缩或膨胀产生的内能变化率。对于 $p>0$ 的流体，压缩时 $\nabla\cdot\boldsymbol{V}<0$，该项为正；膨胀时 $\nabla\cdot\boldsymbol{V}>0$，该项为负。$\Phi$ 表示机械能因黏性作用转化为内能的速率，单位为单位体积的功率。由黏性应力张量的表达式，得

$$ \begin{aligned} \Phi &=\left[2\mu\boldsymbol{E} -\frac{2}{3}\mu(\nabla\cdot\boldsymbol{V})\boldsymbol{I}\right]:\boldsymbol{E}\\ &=2\mu\boldsymbol{E}:\boldsymbol{E} -\frac{2}{3}\mu(\nabla\cdot\boldsymbol{V})^2 \end{aligned} $$

其中

$$ \begin{aligned} \boldsymbol{E}:\boldsymbol{E} &=\varepsilon_{xx}^2+\varepsilon_{yy}^2+\varepsilon_{zz}^2 +2\varepsilon_{xy}^2+2\varepsilon_{yz}^2+2\varepsilon_{zx}^2\\ &=\left(\frac{\partial u}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial z}\right)^2\\ &\quad+\frac{1}{2}\left[ \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\right] \end{aligned} $$

代入耗散函数，得到

$$ \begin{aligned} \Phi &=\mu\left[ \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\right]\\ &\quad+2\mu\left[ \left(\frac{\partial u}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial z}\right)^2 -\frac{1}{3}\left(\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y} +\frac{\partial w}{\partial z}\right)^2\right] \end{aligned} $$

利用恒等式

$$ a^2+b^2+c^2-\frac{1}{3}(a+b+c)^2 =\frac{1}{3}\left[(a-b)^2+(b-c)^2+(c-a)^2\right] $$

可将后三个法向应变率的贡献整理为平方和，得到

$$ \begin{aligned} \Phi &=\mu\left[ \left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\right]\\ &\quad+\frac{2}{3}\mu\left[ \left(\frac{\partial u}{\partial x}-\frac{\partial v}{\partial y}\right)^2 +\left(\frac{\partial v}{\partial y}-\frac{\partial w}{\partial z}\right)^2 +\left(\frac{\partial w}{\partial z}-\frac{\partial u}{\partial x}\right)^2\right] \end{aligned} $$

当 $\mu\ge0$ 时，各项均非负，故

$$ \Phi\ge0 $$

黏性耗散使机械能不可逆地转化为内能。上面的推导适用于有粘性的可压缩流动，具体的说，适用于那些粘性影响不能忽略的高速气流。利用气体的热力学变量之间的关系式，还可以把式（32）改写为由其他的热力学量表示的能量方程。

2、定比热完全气体

对于定容比热容 $c_v$ 为常数的完全气体，单位质量内能可写为

$$ e=c_vT $$

于是

$$ \frac{\mathrm{D}e}{\mathrm{D}t} =c_v\frac{\mathrm{D}T}{\mathrm{D}t} $$

代入牛顿流体的内能方程，得到温度形式的能量方程

$$ \rho c_v\frac{\mathrm{D}T}{\mathrm{D}t} =\rho q+\nabla\cdot(k\nabla T) -p\nabla\cdot\boldsymbol{V}+\Phi $$

温度的随体导数为

$$ \frac{\mathrm{D}T}{\mathrm{D}t} =\frac{\partial T}{\partial t} +u\frac{\partial T}{\partial x} +v\frac{\partial T}{\partial y} +w\frac{\partial T}{\partial z} $$

在直角坐标系中，导热项展开为

$$ \nabla\cdot(k\nabla T) =\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) +\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) +\frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right) $$

因此，完整的分量表达式为

$$ \begin{aligned} \rho c_v\left(\frac{\partial T}{\partial t} +u\frac{\partial T}{\partial x} +v\frac{\partial T}{\partial y} +w\frac{\partial T}{\partial z}\right) &=\rho q +\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) +\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) +\frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right)\\ &\quad-p\left(\frac{\partial u}{\partial x} +\frac{\partial v}{\partial y} +\frac{\partial w}{\partial z}\right)+\Phi \end{aligned} $$

当导热系数 $k$ 在空间中均匀时，导热项可简化为

$$ \nabla\cdot(k\nabla T)=k\nabla^2T $$

3、不可压缩流动

对不可压缩流动，由连续性方程得

$$ \nabla\cdot\boldsymbol{V}=0 $$

故压强的体积变形功为零，内能方程化为

$$ \rho\frac{\mathrm{D}e}{\mathrm{D}t} =\rho q+\nabla\cdot(k\nabla T)+\Phi $$

若采用 $\mathrm{d}e=c_v\,\mathrm{d}T$ 的热物性关系，并取 $\rho,c_v,k$ 为常数，则温度方程为

$$ \rho c_v\frac{\mathrm{D}T}{\mathrm{D}t} =\rho q+k\nabla^2T+\Phi $$

展开为直角坐标形式，得

$$ \rho c_v\left(\frac{\partial T}{\partial t} +u\frac{\partial T}{\partial x} +v\frac{\partial T}{\partial y} +w\frac{\partial T}{\partial z}\right) =\rho q+k\left(\frac{\partial^2T}{\partial x^2} +\frac{\partial^2T}{\partial y^2} +\frac{\partial^2T}{\partial z^2}\right)+\Phi $$

由于速度散度为零，此时耗散函数简化为

$$ \Phi=2\mu\boldsymbol{E}:\boldsymbol{E} $$

其分量表达式为

$$ \begin{aligned} \Phi=\mu\Bigg[&2\left(\frac{\partial u}{\partial x}\right)^2 +2\left(\frac{\partial v}{\partial y}\right)^2 +2\left(\frac{\partial w}{\partial z}\right)^2\\ &+\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)^2 +\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)^2 +\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)^2\Bigg] \end{aligned} $$

不可压缩温度方程的指标形式为

$$ \rho c_v\left(\frac{\partial T}{\partial t} +u_j\frac{\partial T}{\partial x_j}\right) =\rho q+k\frac{\partial^2T}{\partial x_j\partial x_j}+\Phi $$

其中

$$ \Phi=2\mu\varepsilon_{ij}\varepsilon_{ij},\qquad \varepsilon_{ij} =\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right) $$

4、无黏绝热流动

对于无黏流动，$\Phi=0$。若同时没有体积热源与热传导，则内能方程化为

$$ \rho\frac{\mathrm{D}e}{\mathrm{D}t} =-p\nabla\cdot\boldsymbol{V} $$

由连续性方程

$$ \frac{\mathrm{D}\rho}{\mathrm{D}t} =-\rho\nabla\cdot\boldsymbol{V} $$

可得比体积的变化率

$$ \frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{1}{\rho}\right) =-\frac{1}{\rho^2}\frac{\mathrm{D}\rho}{\mathrm{D}t} =\frac{1}{\rho}\nabla\cdot\boldsymbol{V} $$

因此，内能方程可改写为

$$ \frac{\mathrm{D}e}{\mathrm{D}t} +p\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{1}{\rho}\right)=0 $$

记单位质量熵为 $s$。对成分固定、处于局部热力学平衡的简单可压缩流体，热力学关系为

$$ T\,\mathrm{d}s=\mathrm{d}e+p\,\mathrm{d}\left(\frac{1}{\rho}\right) $$

沿流体质点取随体导数，得到

$$ T\frac{\mathrm{D}s}{\mathrm{D}t} =\frac{\mathrm{D}e}{\mathrm{D}t} +p\frac{\mathrm{D}}{\mathrm{D}t}\left(\frac{1}{\rho}\right)=0 $$

故光滑的无黏绝热流动满足

$$ \frac{\mathrm{D}s}{\mathrm{D}t}=0 $$

即每个流体质点的熵沿运动保持不变；若初始熵在全流场相同，则流场保持均匀熵。

若流动还不可压缩，并采用 $\mathrm{d}e=c_v\,\mathrm{d}T$，则

$$ \frac{\mathrm{D}T}{\mathrm{D}t}=0 $$

这表示温度沿质点轨迹不变。只有在初始温度均匀等条件下，才能进一步得到全流场温度为同一常数。