1.2.1流动与传热的通用控制方程（OpenFOAM理论笔记系列）

1.2系统的控制方程

1.2.1流动与传热的通用控制方程

通用的流体流动与传热方程如下：

连续性方程：
∂ρ∂t+∇⋅(ρU⃗)=0(1.4){{\partial \rho}\over{\partial t}}+\nabla\cdot(\rho \vec U)=0 \tag {1.4} ∂t∂ρ+∇⋅(ρU)=0(1.4)

动量方程：
∂(ρU⃗)∂t+∇⋅(ρU⃗U⃗)=∇⋅τ−∇p+ρg⃗(1.5){{\partial (\rho\vec U)}\over{\partial t}}+\nabla\cdot(\rho \vec U\vec U)=\nabla\cdot\tau-\nabla p+\rho \vec g \tag {1.5} ∂t∂(ρU)+∇⋅(ρUU)=∇⋅τ−∇p+ρg(1.5)

能量方程：
∂(ρe)∂t+∇⋅(ρeU⃗)=−∇⋅q⃗−∇⋅(pU⃗)+∇⋅(τ⋅U⃗)+ρg⃗⋅U⃗+ρQ(1.6){{\partial (\rho e)}\over{\partial t}}+\nabla\cdot(\rho e\vec U)=-\nabla\cdot\vec q-\nabla\cdot(p\vec U)+\nabla\cdot(\tau\cdot\vec U)+\rho\vec g\cdot\vec U+\rho Q \tag {1.6} ∂t∂(ρe)+∇⋅(ρeU)=−∇⋅q−∇⋅(pU)+∇⋅(τ⋅U)+ρg⋅U+ρQ(1.6)

状态方程：
p=p(ρ,T),Ｅi=Ei(rho,T)(1.7)p=p(\rho,T)\space,\space Ｅ_i=E_i(rho,T) \tag {1.7} p=p(ρ,T) , Ｅi=Ei(rho,T)(1.7)

接下来我们对上述方程依次进行讨论：

对于连续性方程，ρ\rhoρ为流体的密度，U⃗\vec UU为流体的速度。在不可压缩的情况下，密度不随时空坐标变化，其可化简为：
∇⋅U⃗=0(1.8)\nabla \cdot \vec U =0 \tag{1.8} ∇⋅U=0(1.8)
对于动量方程，τ\tauτ为粘性剪切应力，g⃗\vec gg是重力加速度。当流体不可压缩时，密度可以提到各导数运算符外。为了能够求解动量方程，我们必须知道粘性剪切应力的表达形式，常见的流体为牛顿流体，其粘性剪切应力与形变遵循线性关系。考虑不可压缩流体，其形变率张量可以写为:
S=12(∇U⃗+∇U⃗T)(1.9)\textbf{S}={1\over2}(\nabla\vec U+{\nabla\vec U}^T) \tag{1.9} S=21(∇U+∇UT)(1.9)
粘性剪切应力与剪切率成线性关系，即：
τ=2μS=μ[∇U⃗+(∇U⃗)T](1.10)\tau=2\mu S=\mu\left[\nabla\vec U+(\nabla\vec U)^T\right] \tag{1.10} τ=2μS=μ[∇U+(∇U)T](1.10)
μ\muμ即动力粘度。对于可压缩的牛顿流体，τ\tauτ的表达形式会更加复杂并且要把压力也归入粘性应力的作用中,我们用σ\sigmaσ表示完整的牛顿粘性应力则有：
σ=−(p+23μ∇⋅U⃗)I+μ[∇U⃗+(∇U⃗)T](1.11)\sigma=-(p+{2\over3}\mu\nabla\cdot\vec U)\textbf{I}+\mu\left[\nabla\vec U+(\nabla\vec U)^T\right] \tag{1.11} σ=−(p+32μ∇⋅U)I+μ[∇U+(∇U)T](1.11)
在本文中，我们主要关注不可压缩流体的流动，因此之后我们均采用式(1.10)的表达方式。将式(1.10)带入式(1.5)：
∂(ρU⃗)∂t+∇⋅(ρU⃗U⃗)=∇⋅μ[∇U⃗+(∇U⃗)T]−∇p+ρg⃗(1.12){{\partial (\rho\vec U)}\over{\partial t}}+\nabla\cdot(\rho \vec U\vec U)=\nabla\cdot{\mu\left[\nabla\vec U+(\nabla\vec U)^T\right]}-\nabla p+\rho \vec g \tag{1.12} ∂t∂(ρU)+∇⋅(ρUU)=∇⋅μ[∇U+(∇U)T]−∇p+ρg(1.12)
需要注意的是：
∇U⃗=[∂u∂x∂v∂x∂w∂x∂u∂y∂v∂y∂w∂y∂u∂z∂v∂z∂w∂z](1.13)\nabla \vec{U}=\left[\begin{array}{ccc} \frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} & \frac{\partial w}{\partial x} \\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} & \frac{\partial w}{\partial y} \\ \frac{\partial u}{\partial z} & \frac{\partial v}{\partial z} & \frac{\partial w}{\partial z} \end{array}\right] \tag{1.13} ∇U=⎣⎡∂x∂u∂y∂u∂z∂u∂x∂v∂y∂v∂z∂v∂x∂w∂y∂w∂z∂w⎦⎤(1.13)
∇⋅(∇U⃗T)=∇⋅[∂u∂x∂u∂y∂u∂z∂v∂x∂v∂y∂v∂z∂w∂x∂w∂y∂w∂z]=[∂∂x(∂u∂x)+∂∂y(∂v∂x)+∂∂z(∂w∂x)∂∂x(∂u∂y)+∂∂y(∂v∂y)+∂∂z(∂w∂y)∂∂x(∂u∂z)+∂∂y(∂v∂z)+∂∂z(∂w∂z)]=0(1.14)\nabla \cdot\left(\nabla \vec{U}^{\mathrm{T}}\right)=\nabla \cdot\left[\begin{array}{ccc} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} & \frac{\partial u}{\partial z} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} & \frac{\partial v}{\partial z} \\ \frac{\partial w}{\partial x} & \frac{\partial w}{\partial y} & \frac{\partial w}{\partial z} \end{array}\right]=\left[\begin{array}{c} \frac{\partial}{\partial x}\left(\frac{\partial u}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{\partial v}{\partial x}\right)+\frac{\partial}{\partial z}\left(\frac{\partial w}{\partial x}\right) \\ \frac{\partial}{\partial x}\left(\frac{\partial u}{\partial y}\right)+\frac{\partial}{\partial y}\left(\frac{\partial v}{\partial y}\right)+\frac{\partial}{\partial z}\left(\frac{\partial w}{\partial y}\right) \\ \frac{\partial}{\partial x}\left(\frac{\partial u}{\partial z}\right)+\frac{\partial}{\partial y}\left(\frac{\partial v}{\partial z}\right)+\frac{\partial}{\partial z}\left(\frac{\partial w}{\partial z}\right) \end{array}\right]=0 \tag{1.14} ∇⋅(∇UT)=∇⋅⎣⎢⎡∂x∂u∂x∂v∂x∂w∂y∂u∂y∂v∂y∂w∂z∂u∂z∂v∂z∂w⎦⎥⎤=⎣⎢⎡∂x∂(∂x∂u)+∂y∂(∂x∂v)+∂z∂(∂x∂w)∂x∂(∂y∂u)+∂y∂(∂y∂v)+∂z∂(∂y∂w)∂x∂(∂z∂u)+∂y∂(∂z∂v)+∂z∂(∂z∂w)⎦⎥⎤=0(1.14)
其中u,v,wu,v,wu,v,w是速度在x,y,zx,y,zx,y,z方向上的分量。在式(1.14)中第二个等号后的向量的x分量中，提出∂∂x\frac{\partial}{\partial x}∂x∂后带入不可压缩流动的连续性方程(1.8)进行化简可得0值，其余分量也进行类似处理。
将式1.14带入式1.12，并将密度项提出约去即可得到最常见的不可压缩牛顿流体的动量方程：
∂U⃗∂t+∇⋅(U⃗U⃗)=∇⋅ν∇U⃗−∇pρ+g⃗(1.15){{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U-\nabla {p\over\rho}+ \vec g \tag{1.15} ∂t∂U+∇⋅(UU)=∇⋅ν∇U−∇ρp+g(1.15)
式中ν\nuν即为运动粘度，ν=μ/ρ\nu={\mu/\rho}ν=μ/ρ。

对于能量方程，其等号左边为比能e的全微分。比能由单位质量的机械能Ek=12∣U⃗∣2E_k={1\over 2}|\vec U|^2Ek=21∣U∣2与内能EiE_iEi组成。方程右侧−∇⋅q⃗-\nabla\cdot\vec q−∇⋅q表示热传导对控制体的加热作用,q⃗\vec qq一般可用傅里叶导热定律计算：q⃗=−k∇T\vec q=-k\nabla Tq=−k∇T，其中k为热导率，T为温度。−∇⋅(pU⃗)-\nabla\cdot(p\vec U)−∇⋅(pU)表示压力的功率。∇⋅(τ⋅U⃗)\nabla\cdot(\tau\cdot\vec U)∇⋅(τ⋅U)表示粘性应力的功率，对于牛顿流体可以带入式(1.10)和式(1.11)进行计算。ρg⃗⋅U⃗\rho\vec g\cdot\vec Uρg⋅U代表了重力的功率。ρQ\rho QρQ表示了控制体的内热源带来的热量。

由于动量方程兼有机械能守恒方程的作用，在实际计算中往往将完整的能量方程中的机械能部分抽出，仅组建内能的守恒方程。对于不可压缩流动，为了方便理解，将动量方程(1.5)左侧写成全微分的形式，提出密度项，并左右各自点乘速度矢量可得：
ρdU⃗dt⋅U⃗=(∇⋅τ)⋅U⃗−∇p⋅U⃗+ρg⃗⋅U⃗(1.16)\rho{{d \vec U}\over{d t}}\cdot\vec U=(\nabla\cdot\tau)\cdot\vec U-\nabla p\cdot\vec U+\rho \vec g\cdot\vec U \tag{1.16} ρdtdU⋅U=(∇⋅τ)⋅U−∇p⋅U+ρg⋅U(1.16)
对于方程左侧，为了便于理解，我们在x分量上进行分析：
(ρdU⃗dt⋅U⃗)x=(ρdudt)u=ρd12u2dt(1.17)\left(\rho{{d \vec U}\over{d t}}\cdot\vec U\right)_x=\left(\rho{ d u\over{dt}}\right)u= \rho{{d{1\over2}u^2}\over{dt}} \tag{1.17} (ρdtdU⋅U)x=(ρdtdu)u=ρdtd21u2(1.17)
注意，式(1.17)在不可压缩情况下并不成立：
[d(ρU⃗)dt⋅U⃗]x=[d(ρu)dt]u=(ρdudt+udρdt)u=ρd12u2dt+u2dρdt=d12ρu2dt+12u2dρdt(1.18)\left[{{d (\rho\vec U)}\over{d t}}\cdot\vec U\right]_x=\left[{d (\rho u)\over{dt}}\right]u=\left( \rho{{du}\over{dt}}+u{{d\rho}\over{dt}}\right)u\\= \rho{{d{1\over2}u^2}\over{dt}}+u^2{{d\rho}\over{dt}}={{d {1\over2}\rho u^2}\over{dt}}+{1\over2}u^2{{d\rho}\over{dt}} \tag{1.18} [dtd(ρU)⋅U]x=[dtd(ρu)]u=(ρdtdu+udtdρ)u=ρdtd21u2+u2dtdρ=dtd21ρu2+21u2dtdρ(1.18)
仿照式(1.17)写出其他三个方向的分量并带入式(1.16)可得：
ρdEkdt=(∇⋅τ)⋅U⃗−∇p⋅U⃗+ρg⃗⋅U⃗(1.19)\rho{{d E_k}\over{d t}}=(\nabla\cdot\tau)\cdot\vec U-\nabla p\cdot\vec U+\rho \vec g\cdot\vec U \tag{1.19} ρdtdEk=(∇⋅τ)⋅U−∇p⋅U+ρg⋅U(1.19)
式(1.6)减去式(1.19)可得：
ρ∂Ei∂t+ρ∇⋅(EiU⃗)=−∇⋅q⃗−p∇⋅U⃗+∇⋅(τ⋅U⃗)+ρQ(1.20)\rho{{\partial E_i}\over{\partial t}}+\rho\nabla\cdot( E_i\vec U)=-\nabla\cdot\vec q-p\nabla\cdot\vec U+\nabla\cdot(\tau\cdot\vec U)+\rho Q \tag {1.20} ρ∂t∂Ei+ρ∇⋅(EiU)=−∇⋅q−p∇⋅U+∇⋅(τ⋅U)+ρQ(1.20)
考虑到不可压缩流体的连续性方程(1.8)，忽略粘性耗散以及假设系统无内热源，式(1.20)可转化为：
ρ∂Ei∂t+ρ∇⋅(EiU⃗)=−∇⋅q⃗(1.20)\rho{{\partial E_i}\over{\partial t}}+\rho\nabla\cdot( E_i\vec U)=-\nabla\cdot\vec q \tag {1.20} ρ∂t∂Ei+ρ∇⋅(EiU)=−∇⋅q(1.20)
引入傅里叶导热定律q⃗=−k∇T\vec q=-k\nabla Tq=−k∇T并考虑到Ei=CvTE_i=C_vTEi=CvT，式(1.20)可以转化为更加常用的温度方程：
∂T∂t+∇⋅(TU⃗)=∇⋅(α∇T)(1.20){{\partial T}\over{\partial t}}+\nabla\cdot( T\vec U)=\nabla\cdot(\alpha\nabla T )\tag {1.20} ∂t∂T+∇⋅(TU)=∇⋅(α∇T)(1.20)
上式中α\alphaα为热扩散率α=k/Cvρ\alpha=k/C_v\rhoα=k/Cvρ，式1.20仅在热容在整个场内分布均匀的时候成立，如果热容在场内的空间倒数不为０则不能整合为热扩散率的形式。

状态方程用来封闭方程，最常见的状态方程为理想气体的状态方程，推导温度方程(1.20)时用到的Ei=CvTE_i=C_vTEi=CvT也属于一种状态方程。

本文中主要讨论不可压缩牛顿流体的流动问题，综上所述，其控制方程为：

连续性方程
∇⋅U⃗=0(1.8)\nabla \cdot \vec U =0 \tag{1.8} ∇⋅U=0(1.8)

动量方程
∂U⃗∂t+∇⋅(U⃗U⃗)=∇⋅ν∇U⃗−∇pρ+g⃗(1.15){{\partial \vec U}\over{\partial t}}+\nabla\cdot(\vec U\vec U)=\nabla\cdot\nu\nabla\vec U-\nabla {p\over\rho}+ \vec g \tag{1.15} ∂t∂U+∇⋅(UU)=∇⋅ν∇U−∇ρp+g(1.15)

有关流动与传热的通用控制方程更详细的内容可以参考相关文献¹²³⁴。

系列说明：
接触有限体积法有一段时间了，也看了一些资料，但是有时候总觉得看过一遍之后什么也记不住。老话说得好，眼过千遍不如手过一遍，久而久之我就有了写一些比较像样子的笔记的想法。初学的时候曾经写过一本叫“OpenFOAM编程笔记：单相不可压缩流动”的册子，但当时基础太差（现在基础也不好），错误太多，倒不如推倒重来。
本系列将持续更新，欢迎各位挑错交流，挑错交流可以直接留言也可以联系邮箱cloudbird7@foxmail.com。等到预定内容全部写完后我将集结所有内容为一个独立的文件，开放下载。

Versteeg H K , Malalasekera W . An Introduction to Computational Fluid Dynamic: The Finite Volume Method Second edition [M]. Edinburgh Gate：Pearson Education, 2007. 9-23. ↩︎
Jasak, Hrvoje. Error analysis and estimation for the finite volume method with applications to fluid flows.[D]. London：Imperial College London . 1996 . 63-66 ↩︎
李东岳. 无痛苦N-S方程笔记.[EB/OL]. http://www.dyfluid.com/ , 20200519/20200709 ↩︎
李东岳. CFD中的能量方程.[EB/OL]. http://www.dyfluid.com/docs/energy.html, 20200519/20200709 ↩︎