粘性流体动量方程（NS 方程）推导

接着推导粘性流体的动量守恒方程，也叫 NS 方程。

1. 简介

动量守恒定律也可以描述为：流进控制体的动量和流出控制体的动量的差值与作用在控制体的体积力（对于粘性流体，通常用体积力）之和等于控制体的动量变化率，这里动量流量为标量，规定流入控制体为正，流出控制体为负，动量变化率在数值上等于合力。

2. 动量守恒方程推导

考虑如下控制体，粘性流体和理想流体的区别是表面力除了受压力外，还受切向力的作用。我们做些约定，用 τij 表示表面应力，其中第一个下标 i 代表应力所在平面的外法线方向，第二个下标 j 代表应力的方向，例如 τxy 表示作用在与 x 轴垂直的平面上 y 方向的切应力；τxx 表示作用在与 x 轴垂直的平面上 x 方向的切应力，也就是正应力，参考材料力学，通常也把这个应力标作 σxx。根据切应力互等定律，互换下标的每一对应力是相等的，即 τij=τji，相关的证明过程可以翻阅材料力学教程。另外，我们将应力的增量放在坐标轴的正向上。

控制体的力和应力状态如下图所示，从视图上看，立方体各对面我们命名为前 - 后（法向为 y）；下 - 上（法向为 z）；左 - 右（方向为 x）；

接下来我们根据以上的受力状态分别计算三个方向上的合力。

2.1 X 方向动量方程

所有力和应力中在 x 方向的分量包括下图的红色标识，所有分量之和为：

\[\begin{array}{l} \left(\sigma_{x x}+\frac{\partial \sigma_{x x}}{\partial x} d x\right) d y d z-\sigma_{x x} d y d z+\left(\tau_{y z}+\frac{\partial \tau_{y x}}{\partial y} d y\right) d x d z-\tau_{y x} d x d z+\left(\tau_{z x}+\frac{\partial \tau_{z x}}{\partial z} d z\right) d x d y-\tau_{z x} d x d y \\ +f_{x} d x d y d z=\left(f_{x}+\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right) d x d y d z \end{array}\]

接着，计算控制体 x 方向上的动量流量差值。动量流量定义为：

对于控制体，质量通量来自于六个面，动量均指向 x 方向，亦即速度在 x 的分量方向，如下图。故动量流量差值为：

\[\begin{array}{l} \rho v_{x} v_{x} d y d z-\left[\rho v_{x} v_{x}+\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x} d x\right] d y d z+\rho v_{y} v_{x} d x d z-\left[\rho v_{y} v_{x}+\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y} d y\right] d x d z \\ +\rho v_{z} v_{x} d x d y-\left[\rho v_{z} v_{x}+\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z} d z\right] d x d y \\ =-\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x} d x y d z-\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z} d x d y d z \end{array}\]

X 方向的动量变化率为：

\[\frac{\partial \left( \rho v_x \right)}{\partial t} d x d y d z\]

于是根据动量守恒定律，可以得到 x 方向上的动量方程：

\[\begin{array}{l} \frac{\partial\left(\rho v_{x}\right)}{\partial t} d x d y d z=-\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x} d x d y d z-\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z} d x d y d z+ \\ \left(f_{x}+\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right) d x d y d z \end{array}\]

整理得：

\[\frac{\partial\left(\rho v_{x}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z}=f_{x}+\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right)\]

2.2 Y 方向动量方程

所有力和应力中在 y 方向的分量包括下图的红色标识，所有分量之和为：

\[\begin{array}{l} \left(\sigma_{y y}+\frac{\partial \sigma_{y y}}{\partial y} d y\right) d x d z-\sigma_{y y} d x d z+\left(\tau_{x y}+\frac{\partial \tau_{x y}}{\partial x} d x\right) d y d z-\tau_{x y} d y d z+\left(\tau_{z y}+\frac{\partial \tau_{z y}}{\partial z} d z\right) d x d y-\tau_{z y} d x d y \\ +f_{y} d x d y d z=\left(f_{y}+\frac{\partial \sigma_{y y}}{\partial y}+\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \tau_{z y}}{\partial z}\right) d x d y d z \end{array}\]

接着，计算控制体 y 方向上的动量流量差值。动量流量定义为：

对于控制体，质量通量来自于六个面，动量均指向 y 方向，亦即速度在 y 的分量方向，如下图。故动量流量差值为：

\[\begin{array}{l} \rho v_{x} v_{y} d y d z-\left[\rho v_{x} v_{y}+\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x} d x\right] d y d z+\rho v_{y} v_{y} d x d z-\left[\rho v_{y} v_{y}+\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y} d y\right] d x d z \\ +\rho v_{z} v_{y} d x d y-\left[\rho v_{z} v_{y}+\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z} d z\right] d x d y \\ =-\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x} d x y d z-\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z} d x d y d z \end{array}\]

y 方向的动量变化率为：

\[\frac{\partial \left( \rho v_y \right)}{\partial t} d x d y d z\]

于是根据动量守恒定律，可以得到 y 方向上的动量方程：

\[\begin{array}{l} \frac{\partial\left(\rho v_{y}\right)}{\partial t} d x d y d z=-\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x} d x d y d z-\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z} d x d y d z+ \\ \left(f_{y}+\frac{\partial \sigma_{y y}}{\partial y}+\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \tau_{z y}}{\partial z}\right) d x d y d z \end{array}\]

整理得：

\[\frac{\partial\left(\rho v_{y}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z}=f_{y}+\left(\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}+\frac{\partial \tau_{z y}}{\partial z}\right)\]

2.3 Z 方向动量方程

所有力和应力中在 z 方向的分量包括下图的红色标识，所有分量之和为：

\[\begin{array}{l} \left(\sigma_{z z}+\frac{\partial \sigma_{z z}}{\partial z} d z\right) d x d y-\sigma_{z z} d x d y+\left(\tau_{x z}+\frac{\partial \tau_{x z}}{\partial x} d x\right) d y d z-\tau_{x z} d y d z+\left(\tau_{y z}+\frac{\partial \tau_{y z}}{\partial y} d y\right) d x d z-\tau_{y z} d x d z \\ +f_{z} d x d y d z=\left(f_{z}+\frac{\partial \sigma_{z z}}{\partial z}+\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}\right) d x d y d z \end{array}\]

接着，计算控制体 z 方向上的动量流量差值。动量流量定义为：

对于控制体，质量通量来自于六个面，动量均指向 z 方向，亦即速度在 z 的分量方向，如下图。故动量流量差值为：

\[\begin{array}{l} \rho v_{x} v_{z} d y d z-\left[\rho v_{x} v_{z}+\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x} d x\right] d y d z+\rho v_{y} v_{z} d x d z-\left[\rho v_{y} v_{z}+\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y} d y\right] d x d z \\ +\rho v_{z} v_{z} d x d y-\left[\rho v_{z} v_{z}+\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z} d z\right] d x d y \\ =-\frac{\partial\left(\rho v_{x} v_{z}\right)}{\partial x} d x d y d z-\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z} d x d y d z \end{array}\]

z 方向的动量变化率为：

\[\frac{\partial \left( \rho v_z \right)}{\partial t} d x d y d z\]

于是根据动量守恒定律，可以得到 y 方向上的动量方程：

\[\begin{array}{l} \frac{\partial\left(\rho v_{z}\right)}{\partial t} d x d y d z=-\frac{\partial\left(\rho v_{x} v_{z}\right)}{\partial x} d x d y d z-\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y} d x d y d z-\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z} d x d y d z+ \\ \left(f_{z}+\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \sigma_{z z}}{\partial z}\right) d x d y d z \end{array}\]

整理得：

\[\frac{\partial\left(\rho v_{z}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{z}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z}=f_{z}+\left(\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \sigma_{z z}}{\partial z}\right)\]

综上，获得 x、y、z 三个方向的动量方程组：

\[\begin{array}{l} \frac{\partial\left(\rho v_{x}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z}=f_{x}+\left(\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right) \\ \frac{\partial\left(\rho v_{y}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z}=f_{y}+\left(\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y}+\frac{\partial \tau_{z y}}{\partial z}\right) \\ \frac{\partial\left(\rho v_{z}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{z}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z}=f_{z}+\left(\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \sigma_{z z}}{\partial z}\right) \end{array}\]

显然，上述方程数量小于未知数的数量，方程组不封闭，无解。需要引入其他的方程才能使方程组封闭，在上述方程组的基础上引入了以下本构方程，也就是将应力和速度建立关系。

\[\begin{aligned} \sigma_{x x} &=-p+2 \mu \frac{\partial v_{x}}{\partial x}-\frac{2}{3} \mu\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right) \\ \sigma_{y y} &=-p+2 \mu \frac{\partial v_{y}}{\partial y}-\frac{2}{3} \mu\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right) \\ \sigma_{z z} &=-p+2 \mu \frac{\partial v_{z}}{\partial z}-\frac{2}{3} \mu\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right) \\ \tau_{x y} &=\tau_{y x}=\mu\left(\frac{\partial v_{x}}{\partial y}+\frac{\partial v_{y}}{\partial x}\right) \\ \tau_{x z} &=\tau_{z x}=\mu\left(\frac{\partial v_{x}}{\partial z}+\frac{\partial v_{z}}{\partial x}\right) \\ \tau_{y z} &=\tau_{y z}=\mu\left(\frac{\partial v_{z}}{\partial y}+\frac{\partial v_{y}}{\partial z}\right) \end{aligned}\]

于是联立方程组（10）和（11），可以得到完整形式的 NS 方程组：

\[\begin{array}{l} \frac{\partial\left(\rho v_{x}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{x}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{x}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{x}\right)}{\partial z}=f_{x}-\frac{\partial p}{\partial x}+ \\ \frac{\partial}{\partial x}\left[\mu\left(2 \frac{\partial v_{x}}{\partial x}-\frac{2}{3}\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right)\right)\right]+\frac{\partial}{\partial y}\left[\mu\left(\frac{\partial v_{x}}{\partial y}+\frac{\partial v_{y}}{\partial x}\right)\right]+\frac{\partial}{\partial z}\left[\mu\left(\frac{\partial v_{x}}{\partial z}+\frac{\partial v_{z}}{\partial x}\right)\right] \\ \frac{\partial\left(\rho v_{y}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{y}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{y}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{y}\right)}{\partial z}=f_{y}-\frac{\partial p}{\partial y}+ \\ \frac{\partial}{\partial y}\left[\mu\left(2 \frac{\partial v_{y}}{\partial y}-\frac{2}{3}\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right)\right)\right]+\frac{\partial}{\partial x}\left[\mu\left(\frac{\partial v_{x}}{\partial y}+\frac{\partial v_{y}}{\partial x}\right)\right]+\frac{\partial}{\partial z}\left[\mu\left(\frac{\partial v_{y}}{\partial z}+\frac{\partial v_{z}}{\partial y}\right)\right] \\ \frac{\partial\left(\rho v_{z}\right)}{\partial t}+\frac{\partial\left(\rho v_{x} v_{z}\right)}{\partial x}+\frac{\partial\left(\rho v_{y} v_{z}\right)}{\partial y}+\frac{\partial\left(\rho v_{z} v_{z}\right)}{\partial z}=f_{z}-\frac{\partial p}{\partial z}+ \\ \frac{\partial}{\partial z}\left[\mu\left(2 \frac{\partial v_{z}}{\partial z}-\frac{2}{3}\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{y}}{\partial y}+\frac{\partial v_{z}}{\partial z}\right)\right)\right]+\frac{\partial}{\partial x}\left[\mu\left(\frac{\partial v_{x}}{\partial z}+\frac{\partial v_{z}}{\partial x}\right)\right]+\frac{\partial}{\partial y}\left[\mu\left(\frac{\partial v_{y}}{\partial z}+\frac{\partial v_{z}}{\partial y}\right)\right] \end{array}\]

作者：鱼花生