CFD中的能量方程
质量守恒、动量守恒、能量守恒分别为三大守恒定律。由这些守恒定律可以分别推导出连续性方程、动量方程以及能量方程。对于可压缩瞬态问题，连续性方程可以用来求解密度\(\rho\)，动量方程可以用来求解速度\(\bfU\)，能量方程可以用来获得温度相关量。不可压缩问题的处理要稍微复杂些。能量守恒，顾名思义，定义为绝热系统的总能量是一个常数，总能量不能凭空的产生和消失。能量的单位为\(\mathrm{J}=\mathrm{N}\cdot \mathrm{m}=\frac{\mathrm{kg}\cdot \mathrm{m}^2}{\mathrm{s}^2}\)。在CFD中，可以用下述变量表示能量相关量：
比内能\(e\)：单位质量的内能；
比机械能\(K\)：单位质量的动能，可以定义为\(K=\frac{1}{2}|\mathbf{U}|^2\)；
比总能\(E\)：单位质量的总能量。\(E=K+e\)；
比焓\(h\)：单位质量的焓。\(h=c_p T\)；
上述变量单位均为\(\mathrm{m}^2/\mathrm{s}^2\)。能量守恒表示其中比能\(E\)为守恒变量。在无源项和沉降项的时候，比能其值不随时间变化。只会相互转换。在推导能量方程的过程中，还需要考虑做功的影响。热力学第一定律表明物体内能的增加，等于物体吸热的能量变化和对物体所作的功的总和。考虑一拉格朗日粒子，其每单位体积比能的时间变化率，即为每单位质量的比能的变化率再乘以密度：
(1)\[
\rho\frac{\rD E}{\rD t}=\rho \frac{\mathrm{D}(K+e)}{\mathrm{D}t} \left[\frac{\mathrm{kg}}{\mathrm{s}^3\cdot\mathrm{m}}\right]
\]
方程(1)够成能量方程，也即比能方程的左边。
现考虑热通量导致的能量变化。定义热通量为\(\bfq\)，其单位为\(\mathbf{q}\left[\frac{\mathrm{w}}{\mathrm{m}^2}=\frac{\mathrm{kg}}{\mathrm{s}^3}\right]\)。从物理意义上理解，热通量与连续性方程中的通量\(\rho\bfU\)类似。\(\rho\bfU\)表示速度导致的每单位时间每单位面积的质量变化。\(\bfq\)则表示温度导致的每单位时间每单位面积的能量变化。参考连续性方程，在连续性方程右侧，速度导致的每单位时间、每单位体积的质量变化可以用\(-\nabla\cdot(\rho\bfU)\)来表示。同样的，能量方程中热通量导致的每单位时间、每单位体积的能量变化可以用\(-\nabla\cdot(\bfq)\)来表示。
注意：
热通量（heat flux）为一个矢量。同时，虽然其被称之为热通量，但其表示每单位时间、每单位面积的能量。或许叫能量通量更好。
热通量\(\bfq\)通常并不是CFD关注的主要变量。同时注意到，热通量直接与温度梯度相关，有下述公式：
(2)\[
\bfq=-k\nabla T
\]
其中\(k\)表示热导率。因此，能量方程右侧的热通量相关项\(-\nabla\cdot(\bfq)\)可以写成温度相关的函数\(\nabla\cdot(k\nabla T)\)。
除了热通量会引起能量的变化之外。力对物体做功，必然引起能量的变化。现在回想动量方程的相关项。在动量方程中，单位面积上的应力张量\(\boldsymbol{\sigma}\)导致的动量变化可以表示为\(\nabla\cdot\boldsymbol{\sigma}\)。在能量方程中，单位面积上的应力张量相对应的功率定义为单位面积上的应力张量与速度的乘积\(\boldsymbol{\sigma}\cdot\bfU\)。类似的，能量方程中单位面积上的应力张量导致的能量变化可以表示为\(\nabla\cdot(\boldsymbol{\sigma}\cdot\bfU)\)。现将\(\boldsymbol{\sigma}\)进行分解为压力贡献以及粘性贡献：\(\boldsymbol{\sigma}=\boldsymbol{\tau}-p\bfI\)。有：
(3)\[
\nabla\cdot(\boldsymbol{\sigma}\cdot\bfU)=\nabla\cdot((\boldsymbol{\tau}-p\bfI)\cdot\bfU)
=\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)-\nabla\cdot(p\bfU)
\]
恒等式：
\[
\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)=(\nabla\cdot\boldsymbol{\tau})\cdot\bfU+\boldsymbol{\tau}:\nabla\bfU
\]
至此，能量方程可以表示为：
(4)\[
\rho\frac{\rD E}{\rD t}-\nabla\cdot(k\nabla T)+\nabla\cdot(p\bfU)=\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)
\]
也即：
(5)\[
\frac{\p \rho E}{\p t}+\nabla\cdot(\rho \bfU E)-\nabla\cdot(k\nabla T)+\nabla\cdot(p\bfU)=\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)
\]
其可以继续分解为：
(6)\[
\frac{\p \rho e}{\p t}+\nabla\cdot(\rho \bfU e)
+\frac{\p \rho K}{\p t}+\nabla\cdot(\rho \bfU K)
-\nabla\cdot(k\nabla T)+\nabla\cdot(p\bfU)
=\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)
\]
同时我们应用焓等式：
(7)\[
e=h-p/\rho
\]
将方程(7)代入到(6)中有：
(8)\[
\frac{\p \rho h}{\p t}+\nabla\cdot(\rho \bfU h)
+\frac{\p \rho K}{\p t}+\nabla\cdot(\rho \bfU K)
-\nabla\cdot(k\nabla T)
-\frac{\p p}{\p t}
=\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)
\]
在开源CFD代码OpenFOAM中，热通量项通常并不与温度变量联系起来，而是与比内能\(e\)建立关系：
(9)\[
\bfq=-\alpha_{eff}\nabla e
\]
其中\(\alpha_{eff}\)表示层流热扩散系数以及湍流热扩散系数之和。并且，对于亚音速流动，粘性应力做功项\(\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)\)通常被忽略。
注意：
超音速rhoCentralFoam求解器考虑了\(\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)\)的贡献。
这样，OpenFOAM中植入的能量方程通常具有两种形式，一种是比能方程：
(10)\[
\frac{\p \rho e}{\p t}+\nabla\cdot(\rho \bfU e)
+\frac{\p \rho K}{\p t}+\nabla\cdot(\rho \bfU K)
-\nabla\cdot(\alpha_{eff}\nabla e)+\nabla\cdot(p\bfU)
=0
\]
一种是比焓方程：
(11)\[
\frac{\p \rho h}{\p t}+\nabla\cdot(\rho \bfU h)
+\frac{\p \rho K}{\p t}+\nabla\cdot(\rho \bfU K)
-\nabla\cdot(\alpha_{eff}\nabla e)
-\frac{\p p}{\p t}
=0
\]
另外，如果方程不可压缩或弱可压缩，一般求解比内能方程。比内能方程可以结合动量方程推导而来。首先考虑动量方程
(12)\[
\rho \frac{\mathrm{D} \mathbf{U}}{\mathrm{D} t}=\nabla \cdot \tau-\nabla p
\]
将动量方程中的每个速度分量方程乘以速度分量并加和有：
(13)\[
\rho \frac{\mathrm{D}\left[\frac{1}{2}\left(u_{x}^{2}+u_{y}^{2}+u_{z}^{2}\right)\right]}{\mathrm{D} t}
=\rho\frac{\rD K}{\rD t}
=(\nabla \cdot \tau) \cdot \mathbf{U}-\mathbf{U} \cdot \nabla p
\]
方程(13)与(6)相减有：
(14)\[
\frac{\p \rho e}{\p t}+\nabla\cdot(\rho \bfU e)
+p\nabla\cdot\bfU
-\nabla\cdot(k\nabla T)
=\boldsymbol{\tau}:\nabla\bfU
\]
对于不可压缩或弱可压缩，其中\(\boldsymbol{\tau}:\nabla\bfU\)通常被忽略，因此比内能方程可以写为：
(15)\[
\frac{\p \rho e}{\p t}+\nabla\cdot(\rho \bfU e)
+p\nabla\cdot\bfU
-\nabla\cdot(k\nabla T)
=0
\]
总结：
如果方程不可压缩或弱可压缩，一般求解比内能方程(15)；
如果方程可压缩，一般求解比能方程(10)或者比焓方程(11)；
如果是超音速流动，需要考虑\(\nabla\cdot(\boldsymbol{\tau}\cdot\bfU)\)的贡献；