流体力学守恒方程

贡献者： _Eden_

预备知识　流体运动的描述方法，物质导数（实质导数），重积分、面积分、体积分（简明微积分）

　　人们一般用密度 $ρ$ ，速度 $v$ 等物理量描述流体中每个位置流体微元的性质，并列出所谓的流体力学方程来描述每个流体微元随时间的演化。它们之间应当满足一定的守恒方程，例如连续性方程，动量守恒方程，能量守恒方程等。在这一节中，我们将对流体力学的守恒方程进行简要的推导和介绍。

1. 连续性方程（流守恒方程）

　　让我们对物质体（Material Volumn）进行分析，取一个随流体一起运动的被闭曲面 $S$ 包围的体积 $V$ ，对其中的密度进行积分得到这一个物质体的质量。物质体在随流体一起运动的过程中，其质量应当是不变的，那么我们有

\begin{matrix} (1) & \frac{d}{d t} \int_{V (t)} ρ d V = 0 . \end{matrix}

进一步将它拆成两部分，一部分是每个位置密度的偏导数随时间的变化，另一部分是闭合曲面

S

的变化带来的质量变化。我们得到

\begin{matrix} (2) & \int \frac{\partial ρ}{\partial t} d V + \int_{S} ρ u \cdot \hat{n} d S = 0 \Rightarrow \int [\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u)] d V = 0 . \end{matrix}

如上，我们得到了连续性方程的微分形式

\begin{matrix} (3) & \frac{\partial ρ}{\partial t} + \nabla \cdot (ρ u) = 0 . \end{matrix}

利用

\nabla \cdot (ρ u) = ρ \nabla \cdot u + u \cdot \nabla ρ

，上式可以改写为

\begin{matrix} (4) & \frac{d ρ}{d t} + ρ \nabla \cdot u = 0, \end{matrix}

我们也可以用另一种方式来考察这一结果。直接从

\frac{d}{d t} (ρ δ V) = 0

出发（这表明一个流体微元随时间演化的过程中质量是不变的），得到

\begin{matrix} (5) & ρ \frac{d δ V}{d t} + δ V \frac{d ρ}{d t} = 0 \Rightarrow \frac{d δ V}{d t} = - δ V \frac{1}{ρ} \frac{d ρ}{d t} . \end{matrix}

利用式 4 ，上式得到了以下结果：

\begin{matrix} (6) & \frac{1}{δ V} \frac{d δ V}{d t} = \nabla \cdot u, \end{matrix}

这个结果是不难理解的。取

δ V

为一个无限小立方体，计算六个面流体运动的流进流出，可以轻易得到这个表达式，然后事实上不论

δ V

的形状如何，上式都是成立的，因为我们总能把

δ V

分割成许许多多个小立方体，而每个小立方体都满足上式，进而利用

\nabla \cdot u

在这一小区域的连续性，就可以得到这个结果。

　　类似上面的推导，我们先对物质体列出积分方程：

\begin{matrix} (7) & \frac{d}{d t} \int_{V (t)} ρ u d V = \int ρ g d V + \int_{S} f d S . \end{matrix}

其中

g

表示作用于每个流体微元上的体力，例如重力就是体力的一种。对上式左侧进行化简，可以得到

\begin{matrix} (8) & \int \frac{\partial (ρ u)}{\partial t} d V + \int_{S} ρ u (u \cdot \hat{n}) d S = \int ρ g d V + \int_{S} \hat{n} \cdot \overset{\leftrightarrow}{T} d S . \end{matrix}

其中

T

是流体的应力张量，

T_{i j}

表示在法线为

x_{i}

方向的单位面元上，面外对面内的面力的

x_{j}

分量。可以证明，为了保证角动量守恒，应力张量是二阶对称的张量。最终我们有矢量表达式

\begin{matrix} (9) & ρ \frac{\partial u}{\partial t} + ρ (u \cdot \nabla) u = ρ g + \nabla \cdot \overset{\leftrightarrow}{T}, \end{matrix}

　　或者我们有一种更简单的推导方法。注意到 $d (ρ u δ V) / d t = ρ δ V d u / d t + u d (ρ δ V) / d t = ρ δ V d u / d t$ 。所以从式 7 出发立刻可以得到

\begin{matrix} (10) & ρ \frac{d u}{d t} = ρ g + \nabla \cdot \overset{\leftrightarrow}{T}, \end{matrix}

这和式 9 是等价的。

　　应力张量 $T_{i j}$ 一般表示为 $- p + τ_{i j}$ ， $p$ 表示压强， $τ_{i j}$ 则是由流体的粘性带来的项， $τ_{i j}$ 是对称张量。

　　最后让我们讨论能量守恒方程。流体微元的能量一般分为内能 $e$ 和动能 $\frac{1}{2} ρ u^{2}$ ，前者和流体的热力学状态（压强、温度或密度）有关，而后者和流体微元的运动有关。下面我们来具体分析能量守恒方程的形式。

\begin{matrix} (11) & \frac{d}{d t} \int_{V (t)} ρ (e + \frac{1}{2} | u |^{2}) d V = \int ρ g \cdot u d V + \int_{S} f \cdot u d S - \int_{S} q \cdot n d A . \end{matrix}

e

为单位质量的内能。面力

f

在速度方向上作用一段距离导致做功，

q

表示面元上的热流方向。

　　那么经过一系列的化简我们可以得到

\begin{matrix} (12) & \begin{aligned} ρ \frac{d}{d t} (e + \frac{1}{2} | u |^{2}) & = ρ g \cdot u + \frac{\partial}{\partial x_{i}} (T_{i j} u_{j}) - \nabla \cdot q \\ = ρ g_{i} u_{i} + (- p \frac{\partial u_{j}}{\partial x_{j}} + τ_{i j} \frac{\partial u_{j}}{\partial x_{i}}) + u_{j} (- \frac{\partial p}{\partial x_{j}} + \frac{\partial τ_{i j}}{\partial x_{i}}) - \frac{\partial q_{i}}{\partial x_{i}}, \end{aligned} \end{matrix}

以上就是我们得到的能量守恒方程。下面我们对动量守恒方程点乘

u

，来考察能量守恒方程中的动能守恒

\begin{matrix} (13) & ρ \frac{d}{d t} (\frac{1}{2} | u |^{2}) = ρ g_{i} u_{i} + u_{j} \frac{\partial T_{i j}}{\partial x_{i}} = ρ g_{i} u_{i} + u_{j} (- \frac{\partial p}{\partial x_{j}} + \frac{\partial τ_{i j}}{\partial x_{i}}) . \end{matrix}

将式 12 和式 13 相减，我们得到了内能守恒：

\begin{matrix} (14) & \begin{aligned} ρ \frac{d e}{d t} & = T_{i j} \frac{\partial u_{j}}{\partial x_{i}} - \frac{\partial q_{i}}{\partial x_{i}} \\ = - p \frac{\partial u_{j}}{\partial x_{j}} + τ_{i j} \frac{\partial u_{j}}{\partial x_{i}} - \frac{\partial q_{i}}{\partial x_{i}} = - \frac{p}{ρ} \frac{d ρ}{d t} - \nabla \cdot q \\ \Rightarrow \frac{d e}{d t} & = p \frac{d ν}{d t} + τ_{i j} S_{i j} - \nabla \cdot q, (ν = 1 / ρ) . \end{aligned} \end{matrix}

　　如果将应力张量的表达式式 2 代入，并利用热传导定律将 $q$ 表达为 $- k \nabla T$ ，上述方程可以改写为

\begin{matrix} (15) & ρ \frac{d e}{d t} = - p \frac{\partial u_{m}}{\partial x_{m}} + 2 μ {(S_{i j} - \frac{1}{3} \frac{\partial u_{m}}{\partial x_{m}} δ_{i j})}^{2} + μ_{ν} {(\frac{\partial u_{m}}{\partial x_{m}})}^{2} + \frac{\partial}{\partial x_{i}} (k \frac{\partial T}{\partial x_{i}}) . \end{matrix}