CFD中的严格温度方程

本文仅发布在博客园(https://www.cnblogs.com/josephliu2300/p/19124423 )和github仓库。未经许可禁止转载。如果你在其他网站见到本文，那必定是违规转载，我严厉谴责这种生儿子没屁眼的行为。
mynw.cn 在发布之后6分钟就违规转载了我的文章，悬挂大量广告，而且该网站不提供任何沟通、申诉渠道。该网站在文末标注的出处也不是我发布的博客园，而是它本身的网站，这是十分无耻的行为，对此表示严厉谴责。
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}) }\]

来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！