CFD中通常求解能量方程以模拟温度变化，这需要在比能/比焓和温度之间的转化。但如果问题包含多个相，且多相之间具有相同的温度（热平衡假设），则没有明确的能量方程可以求解。此时只能求解温度方程。作为基础，本文推导了严格的单相温度方程。
从焓方程开始
在OpenFOAM的fluid模块中，比焓\(h\)的能量方程如此植入：
\[\begin{align}
\frac{\partial (\rho h)}{\partial t} + \nabla\cdot(\rho U h) + \frac{\partial(\rho K)}{\partial t} + \nabla\cdot(\rho U K) - \frac{\partial p}{\partial t} + \nabla\cdot q = S_h
\end{align}
\]
比焓的全导数
\[\begin{align}
\mathrm{d}h = (\frac{\partial h}{\partial T})_p \mathrm{d}T + (\frac{\partial h}{\partial p})_T \mathrm{d}p
\end{align}
\]
第一个偏导数是定压比热\(C_p\)，第二个偏导数相对抽象，它代表等温压缩造成的能量变化，记为\(Ж\)。
将比焓的全微分带入比焓方程，逐项拆解。
瞬态项和对流项变为：
\[\begin{align}\frac{\partial (\rho h)}{\partial t} + \nabla\cdot(\rho U h) &= \underbrace{h \frac{\partial \rho}{\partial t} + h \nabla\cdot(\rho U)}_{=0} + \rho\frac{\partial h}{\partial t} + \rho U\cdot\nabla h \\
&= \rho C_p \frac{\partial T}{\partial t} + \rho Ж\frac{\partial p}{\partial t} + \rho U\cdot C_p\nabla T + Ж\rho U\cdot\nabla p \\
&= \rho C_p \frac{\partial T}{\partial t} + \underbrace{\nabla\cdot (\rho U C_p T)}_\text{fvm::div} - \underbrace{T \nabla\cdot(\rho U C_p)}_\text{fvm::SuSp} + Ж(\rho\frac{\partial p}{\partial t}+\rho U\cdot\nabla p)
\end{align}
\]
扩散项变为：
\[\begin{align}
\nabla\cdot q=-\nabla\cdot(\kappa\nabla T)
\end{align}
\]
瞬态项、扩散项采用隐式离散，其余没有标注的项均采用显式离散。
从状态方程计算\(Ж\)
\[\begin{align}
(\frac{\partial H}{\partial p})_T &= (\frac{\partial S}{\partial p})_T + V \\
&=V - T (\frac{\partial V}{\partial T})_p \\
h &= \frac{H}{M_w} \\
Ж &= \frac{1}{M_w}(V - T (\frac{\partial V}{\partial T})_p) \\
&= \frac{1}{\rho} + \tilde{R}(z + T(\frac{\partial z}{\partial T})_p)
\end{align}
\]
式中\(H\)的单位是J/mol，\(h\)的单位是J/kg。
显然对于理想气体，\(Ж =0\)。对于非理想气体，需要从状态方程得到\((\frac{\partial V}{\partial T})_p\)的表达式。
RK方程
对于RK方程，
\[(\frac{\partial V}{\partial T})_p = \frac{\frac{R V^2}{p} + (\frac{a}{2 p \sqrt{T}} + \frac{R b}{p})V + \frac{ab}{2 p T \sqrt{T}}}{ 3V^2 - 2\frac{RTV}{p} + (\frac{a}{p\sqrt{T}}-b^2-\frac{bRT}{p}) }
\]
本文章来自 https://github.com/nnSemenov/Mikeno-13/blob/master/doc/rigorousTemperatureEquation.typ ，作者是本人。