3 有限体积法：推导方程
3 有限体积法：推导方程
基本原理和目标
（注意：这一节看不懂没关系，在后面的推导中会慢慢用到）
-
质量、动量和能量的守恒
- 流体的质量守恒
- 动量改变的速度 = 一个流体粒子上受到的力的总和（牛顿第二定律）
- 能量改变的速度 = 一个流体粒子吸收的热量，和作用在其上的功的总和（热力学第一定律）
-
推导出控制流动的偏微分方程
-
通过牛顿模型描述粘性应力，推导出纳维-斯托克斯方程
-
控制流体行为的方程和传输方程是相似的
-
把传输方程的积分形式应用于有限时间间隔和有限控制体积
-
把物理现象分为椭圆型、抛物线型、双曲线型
-
将流体看作一个连续体（continuum）
-
考虑一个小流体元素，其边长为\(\delta x\)，\(\delta y\)，\(\delta z\)
三维物理量守恒
质量守恒
一个流体元素中质量的增加速度 = 流入流体元素的净流动速度
质量增加速度：
\[\begin{align}\frac{\partial}{\partial t}(\rho\delta x\delta y\delta z)=\frac{\partial\rho}{\partial t}(\delta x\delta y\delta z)\end{align}
\]
注意：这里\(\rho\delta x\delta y\delta z\)可理解为\(\rho \delta V\)，也就是代表了该流体微元的总质量。
然后求质量关于时间的一阶导，就是质量的净增加速度了。
穿过流体元素边界，流入流体元素的净质量流动速率为：
其中，约定流入流体微元的质量流动标正号，流出的标负号。
注意这里是如何运用泰勒级数展开的（以x方向，即图中横着的方向为例）：
-
在点x处的质量流动速率密度（即单位面积的质量流率）为\(\rho u\)；
-
使用一阶泰勒级数展开：
\[\begin{align}f(x+\delta x)=f(x)+\frac{df(x)}{dx} \delta x\end{align} \] -
在\(x-\frac{1}{2} \delta x\)（即流入处）：
\[\begin{align} f(x-\frac{1}{2}\delta x)&=f(x)-\frac{\partial f(x)}{\partial x}\cdot \frac{1}{2}\delta{x}\\ &=\rho u-\frac{\partial (\rho u)}{\partial x}\cdot \frac{1}{2}\delta{x} \end{align} \]在\(x+\frac{1}{2} \delta x\)（即流出处）同理，此处不再展示。
-
x方向流入的真正质量流动速率为（即乘上个截面面积）：
\[\begin{align}\left(\rho u-\frac{\partial(\rho u)}{\partial x}\frac12\delta x\right)\delta y\delta z\end{align} \]
把六个方向（x，y，z方向的流入和流出）加一下，可得：
\[\begin{align}
\begin{gathered}
\left(\rho u-\frac{\partial(\rho u)}{\partial x}\frac12\delta x\right)\delta y\delta z-\left(\rho u+\frac{\partial(\rho u)}{\partial x}\frac12\delta x\right)\delta y\delta z \\
+\left(\rho v-\frac{\partial(\rho v)}{\partial y}\frac{1}{2}\delta y\right)\delta x\delta z-\left(\rho v+\frac{\partial(\rho v)}{\partial y}\frac{1}{2}\delta y\right)\delta x\delta z \\
+\left(\rho w-\frac{\partial(\rho w)}{\partial z}\frac{1}{2}\delta z\right)\delta y\delta x-\left(\rho w+\frac{\partial(\rho w)}{\partial z}\frac{1}{2}\delta z\right)\delta y\delta x
\end{gathered}\end{align}
\]
化简可得最终穿过流体元素边界，流入流体元素的净质量流动速率：
\[\begin{align}\left(\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}\right)\delta x\delta y\delta z\end{align}
\]
那么，我们可以得到质量守恒方程：
\[\begin{align}
\frac{\partial\rho}{\partial t}(\delta x\delta y\delta z)+\left(\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}\right)\delta x\delta y\delta z=0\end{align}
\]
把\(\delta x \delta y \delta z\)消掉，可得：
\[\begin{align}&\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u)}{\partial x}+\frac{\partial(\rho v)}{\partial y}+\frac{\partial(\rho w)}{\partial z}=0
\\&\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\boldsymbol{u})=0\end{align}
\]
注意，这里的Nabla算子：
\[\begin{align}
\nabla=\frac{\partial}{\partial x}\vec{i}+\frac{\partial}{\partial y}\vec{j}+\frac{\partial}{\partial z}\vec{k}
\end{align}
\]
上述方程可表示在可压缩流体中的某一点，非稳态、三维流动下的质量守恒或连续性方程。
而对于不可压缩流动，密度是守恒的，此时：
\[\begin{align}
&\nabla\cdot\boldsymbol{u}=0\\
\quad&\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0
\end{align}
\]
推广：物理量守恒
现设有一物理量在每单位质量的值为\(\phi\)。
沿着一个流体粒子的轨迹，对函数\(\phi\)的总（total）或实质（substantive）导数可表示为：
\[\begin{align}
\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial t}+\frac{\partial\phi}{\partial x}\frac{dx}{dt}+\frac{\partial\phi}{\partial y}\frac{dy}{dt}+\frac{\partial\phi}{\partial z}\frac{dz}{dt}\end{align}
\]
注意，流体粒子会随着流动而移动，故：
\[\begin{align}
\frac{dx}{dt}=u,~\frac{dy}{dt}=v,~\frac{dz}{dt}=w
\end{align}
\]
由此可得，\(\phi\)的实质导数可写作：
\[\begin{align}
\frac{D\phi}{Dt}&=\frac{\partial\phi}{\partial t}+u\frac{\partial\phi}{\partial x}+v\frac{\partial\phi}{\partial y}+w\frac{\partial\phi}{\partial z}\\
&=\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi
\end{align}
\]
该方程表示，该物理量随着时间的变化+流体粒子在空间中的运动与物理量在空间中变化的叠加效应=单位质量物质的\(\phi\)属性沿着一个流体微元运动轨迹的总变化率。
那么，对于一个流体粒子而言，其物理量\(\phi\)每单位体积的变化速度为：
\[\begin{align}
\rho\frac{D\phi}{Dt}=\rho\left(\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi\right)
\end{align}
\]
下一步，我们向任意守恒的物理量（arbitrary conserved property）推广（generalize）出其守恒方程：
（注意：别直接用这个方程，这个是来自PPT的，请参见下面的方程27）
\[\begin{align}
\frac{\partial(\rho\phi)}{\partial t}+\nabla\cdot(\rho\phi\boldsymbol{u})=\rho\left[\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi\right]+\phi\left[\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\boldsymbol{u})\right]=\rho\frac{D\phi}{Dt}
\end{align}
\]
其具体推导过程为：
-
原实质导数参见方程17：\(\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi\)
-
首先，推广到密度和物理量一起，即考虑\(\rho \phi\)随时间的变化率，此时使用乘积法则：
\[\begin{align} d(uv)=udv+vdu \end{align} \]利用该法则，有：
\[\begin{equation} \frac{\partial(\rho\phi)}{\partial t}=\phi\frac{\partial\rho}{\partial t}+\rho\frac{\partial\phi}{\partial t} \end{equation} \] -
下一步，考虑\(\rho \phi\)在空间上的变化，加入流体运动速度\(\boldsymbol{u}\)（这里加粗代表x，y，z三个方向的速度矢量），即对\(\rho \phi \boldsymbol{u}\)取散度，三个物理量都可能随着位置变化：
\[\begin{equation} \nabla\cdot(\rho\phi\vec{u})=\phi[\nabla\cdot(\rho\vec{u})]+(\rho\vec{u})\cdot\nabla\phi \end{equation} \] -
将过程2和过程3的结果，即时间导数和空间上的散度项结合起来：
\[\begin{align} \frac{\partial(\rho\phi)}{\partial t}+\nabla\cdot(\rho\phi\vec{u})=\left(\phi\frac{\partial\rho}{\partial t}+\rho\frac{\partial\phi}{\partial t}\right)+(\phi(\nabla\cdot(\rho\vec{u}))+\rho\vec{u}\cdot\nabla\phi) \end{align} \] -
根据连续性方程，质量守恒要求：
\[\begin{equation} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\vec{u})=0 \end{equation} \]注意到方程23中：
\[\begin{equation} \phi \frac{\partial\rho}{\partial t}+\phi \cdot \nabla\cdot(\rho\vec{u})=0 \end{equation} \]故结合后的方程可写作：
\[\begin{equation} \frac{\partial(\rho\phi)}{\partial t}+\nabla\cdot(\rho\phi\vec{u})=\rho\frac{\partial\phi}{\partial t}+\rho\vec{u}\cdot\nabla\phi \end{equation} \] -
注意到方程26中右侧项恰等于\(\rho\)乘以\(\frac{D\phi}{Dt}\)。故可将推广后的守恒方程写作：
\[\begin{equation} \frac{\partial(\rho\phi)}{\partial t}+\nabla\cdot(\rho\phi\vec{u})=\rho\left(\frac{\partial\phi}{\partial t}+\vec{u}\cdot\nabla\phi\right)=\rho\frac{D\phi}{Dt} \end{equation} \]其物理意义为：
\[物理量\phi随时间的变化率+净流出速率=总变化率 \]
动量守恒
牛顿第二定律表明，一个流体粒子动量的改变=作用在流体粒子上力的总和。
对于一个流体粒子，在x, y, z方向，每单位体积动量增加的速度为：
\[\begin{equation}
\rho\frac{Du}{Dt},~\rho\frac{Dv}{Dt},~\rho\frac{Dw}{Dt}
\end{equation}
\]
注意这里：动量\(P=mv\)（这里v不代表特定于y方向的速度，而是泛指速度），每单位体积动量即为\(\rho v\)，则每单位体积动量增加的速度即为密度\(\times\)速度关于时间的导数。
作用于流体元上的力：
-
表面力：作用于流体表面，由于流体元与其周围流体或边界的接触产生。
-
压力力（pressure forces）
- 由于流体中的压力梯度作用于流体元表面产生的力；
- 垂直于流体元表面，指向外界（不过这里好像也不一定，要看具体流动是出微元还是进）；
- 导致流体加减速，是推动流体流动的重要因素。
-
黏性力（viscous forces）
- 由于流体的黏性（流体分子之间的内摩擦），在流体元表面产生的剪切应力；
- 在流体层之间起阻力作用，特别是在层流和接触边界处；
- 造成速度梯度（摩擦-相对速度），使得流体从高速区域向低速区域传递动量，最终影响整个流体的流动特性；
- 补充：摩擦力是分子间的吸引力（如范德华力）、原子间的电磁排斥力、原子间电子交换或共享形成弱键合（如氢键）等等电磁力在宏观下的体现。
-
-
体力（body forces，不是HP）：外界力场直接作用于流体的每个体积单元内部。
- 重力，影响流体压力分布、密度梯度等。
- 离心力（centrifugal force）：在旋转系统（如地球自转引起的大气、海洋环流）中，由于旋转参考系下的惯性效应产生的力，指向远离旋转轴的方向。
- 科里奥利力（Coriolis force）：旋转参考系下的虚拟力，在地球表面上主要由地球自转引起。在北半球，流体运动向右偏；在南半球，流体运动向左偏。
- 电磁力（electromagnetic force）：带电流体元在电场和磁场中所受的洛伦兹力，在工程力学中常用于控制电解液、等离子体的运动。
一般地，在动量方程中，我们将表面力（压力和黏性力）视为主要项并将它们分开考虑，而将体力视为源项。
通过压力和9个粘性应力分量，我们定义了流体元素的应力状态。
该图是为x方向动量方程准备的，故只画出了6个表面沿x方向走的力。
注意，图中所示的物理量变化速率仍然运用了一阶泰勒级数展开，具体推导过程见质量守恒部分，方程2-5：
\[\begin{align}f(x+\delta x)=f(x)+\frac{df(x)}{dx} \delta x\end{align}
\]
考虑x方向，由于压力P和应力分量\(\tau_{xx}\)，\(\tau_{yx}\)和\(\tau_{zx}\)造成的总力。
如上图所示，考虑一对面（E，W）（东/右侧壁面，西/左侧壁面，被x轴穿过的一对，上北下南左西右东），把图示力加一下，得：
\[\begin{equation}
\begin{aligned}
&\left[\left(p-\frac{\partial p}{\partial x}\frac{1}{2} \delta x\right)-\left(\tau_{xx}-\frac{\partial\tau_{xx}}{\partial x}\frac{1}{2} \delta x\right)\right]\delta y\delta z \\
&+\left[-\left(p+\frac{\partial p}{\partial x}\frac{1}{2}\delta x\right)+\left(\tau_{xx}+\frac{\partial\tau_{xx}}{\partial x}\frac{1}{2}\delta x\right)\right]\delta y\delta z\\&=\left(-\frac{\partial p}{\partial x}+\frac{\partial\tau_{xx}}{\partial x}\right)\delta x\delta y\delta z
\end{aligned}
\end{equation}
\]
注意，这里考虑了4个力：两个表面上的压力和黏性力，而下面两个方程是不考虑压力的，因为压力垂直于流体元表面 。
类似地，对于另一对面（N，S）（上面和底下面，被y轴穿过），有：
\[\begin{equation}
-\left(\tau_{yx}-\frac{\partial\tau_{yx}}{\partial y}\frac{1}{2}\delta y\right)\delta x\delta z+\left(\tau_{yx}+\frac{\partial\tau_{yx}}{\partial y}\frac{1}{2}\delta y\right)\delta x\delta z=\frac{\partial\tau_{yx}}{\partial y}\delta x\delta y\delta z
\end{equation}
\]
而对于（T，B）（前面和后面，被z轴穿过），有：
\[\begin{equation}
-\left(\tau_{zx}-\frac{\partial\tau_{zx}}{\partial z}\frac{1}{2}\delta z \right)\delta x\delta y+\left(\tau_{zx}+\frac{\partial\tau_{zx}}{\partial z}\frac{1}{2}\delta z \right)\delta x\delta y=\frac{\partial\tau_{zx}}{\partial z}\delta x\delta y\delta z
\end{equation}
\]
把上述方程30、31、32相加，可得由上述表面力导致的，该流体微元每单位体积受到的总力，也就是该流体微元在x方向的每单位体积动量变化速度：
\[\begin{equation}
\rho \frac{Du}{Dt}=\frac{\partial(-p+\tau_{xx})}{\partial x}+\frac{\partial\tau_{yx}}{\partial y}+\frac{\partial\tau_{zx}}{\partial z}
\end{equation}
\]
类似地，y方向的动量方程为：
\[\begin{equation}
\rho \frac{Dv}{Dt}=\frac{\partial\tau_{xy}}{\partial x}+\frac{\partial\left(-p+\tau_{yy}\right)}{\partial y}+\frac{\partial\tau_{zy}}{\partial z}
\end{equation}
\]
z方向的动量方程：
\[\begin{equation}
\rho\frac{Dw}{Dt}=\frac{\partial\tau_{xz}}{\partial x}+\frac{\partial\tau_{yz}}{\partial y}+\frac{\partial(-p+\tau_{zz})}{\partial z}
\end{equation}
\]
能量守恒
基于热力学第一定律：
\[\begin{equation}
\begin{aligned}
流体微元能量增加速度=\\净热量增加速度+作用于流体微元上的净功增加速度
\end{aligned}
\end{equation}
\]
定义一个流体微元/流体粒子每单位体积的能量变化速率为：
\[\rho \frac{DE}{Dt}
\]
做功
接下来，以x方向为例：
\[\begin{equation}
\begin{aligned}
&作用于流体粒子上的表面力所做功的速率=\\
&力\times 力作用方向上的速度
\end{aligned}
\end{equation}
\]
这样，我们把表示x方向合力的方程30乘以x方向速度分量\(u\)，得：
\[\begin{equation}
\begin{aligned}
&\left[\left(pu-\frac{\partial(pu)}{\partial x}\frac{1}{2} \delta x\right)-\left(\tau_{xx}u-\frac{\partial(\tau_{xx}u)}{\partial x}\frac{1}{2} \delta x\right)-\left(pu+\frac{\partial(pu)}{\partial x}\frac{1}{2} \delta x\right)\right] \\
&\left.+\left(\tau_{xx}u+\frac{\partial(\tau_{xx}u)}{\partial x}\frac{1}{2}\delta x\right)\right]\delta y\delta z \\
&+\left[-\left(\tau_{yx}u-\frac{\partial\left(\tau_{yx}u\right)}{\partial y}\frac{1}{2}\delta y\right)+\left(\tau_{yx}u+\frac{\partial\left(\tau_{yx}u\right)}{\partial y}\frac{1}{2}\delta y\right)\right]\delta x\delta z \\
&+\left[-\left(\tau_{zx}u-\frac{\partial(\tau_{zx}u)}{\partial z}\frac{1}{2} \delta z\right)+\left(\tau_{zx}u+\frac{\partial(\tau_{zx}u)}{\partial z}\frac{1}{2} \delta z\right)\right]\delta x\delta y
\end{aligned}
\end{equation}
\]
加减乘除一下，可得x方向表面力的净功率为：
\[\begin{equation}
\left[\frac{\partial\left(u(-p+\tau_{xx})\right)}{\partial x}+\frac{\partial\left(u\tau_{xy}\right)}{\partial y}+\frac{\partial(u\tau_{xz})}{\partial z}\right]\delta x\delta x\delta y
\end{equation}
\]
类似地，有y和z方向表面力对x方向流动的净作用功率：
\[\begin{equation}
\left[\frac{\partial\left(v\tau_{xy}\right)}{\partial x}+\frac{\partial\left(v\left(-p+\tau_{yy}\right)\right)}{\partial y}+\frac{\partial\left(v\tau_{yz}\right)}{\partial z}\right]\delta x\delta y\delta z
\end{equation}
\]
\[\begin{equation}
\left[\frac{\partial(v\tau_{xz})}{\partial x}+\frac{\partial(v\tau_{xy})}{\partial y}+\frac{\partial\left(v(-p+\tau_{zz})\right)}{\partial z}\right]\delta x\delta y\delta z
\end{equation}
\]
将上述方程39、40、41相加，可得表面力对流体粒子作用的总功率为：
\[\begin{equation}\begin{aligned}
&-\nabla\cdot(p\boldsymbol{u}) \\
&+\frac{\partial(u\tau_{xx})}{\partial x}+\frac{\partial\big(u\tau_{xy}\big)}{\partial y}+\frac{\partial(u\tau_{xz})}{\partial z}\\&+\frac{\partial\big(v\tau_{xy}\big)}{\partial x}+\frac{\partial\big(v\tau_{yy}\big)}{\partial y}+\frac{\partial\big(v\tau_{yz}\big)}{\partial z}\\&+\frac{\partial(w\tau_{xz})}{\partial x}+\frac{\partial\big(w\tau_{yz}\big)}{\partial y}
+\frac{\partial(w\tau_{zz})}{\partial z}
\end{aligned}\end{equation}
\]
热
由于热流穿过流体微元边界，流体微元每单位体积的净热增加速度为：
\[\begin{equation}
-\frac{\partial q_x}{\partial x}-\frac{\partial q_y}{\partial y}-\frac{\partial q_z}{\partial z}=-\nabla\cdot \vec{q}
\end{equation}
\]
此时，我们引入傅里叶定律：
\[\begin{equation}
\vec{q}=-k\nabla T
\end{equation}
\]
这里：
- \(\vec{q}\)是热流密度向量，表示热量通过单位面积的速率；
- k是材料的热导率；
- \(\nabla T\)是温度梯度。
这样，我们能重写方程44的推导结果为：
\[\begin{equation}
-\nabla \cdot \vec{q}=\nabla \cdot (k\nabla T)
\end{equation}
\]
动能
为便于进一步推导内能i或温度T的方程，我们建立动能变化的方程。
该方程实际上就是把动量方程再乘上速度\(\vec{u}\)。
\[\begin{equation}
\begin{aligned}
&\rho\frac{DK}{Dt} \\
&=- \boldsymbol{u}\cdot\nabla p+u\left(\frac{\partial\tau_{xx}}{\partial x}+\frac{\partial\tau_{xy}}{\partial y}+\frac{\partial\tau_{xz}}{\partial z}\right)+v\left(\frac{\partial\tau_{yx}}{\partial x}+\frac{\partial\tau_{yy}}{\partial y}+\frac{\partial\tau_{yz}}{\partial z}\right) \\
&+w\left(\frac{\partial\tau_{zx}}{\partial x}+\frac{\partial\tau_{zy}}{\partial y}+\frac{\partial\tau_{zz}}{\partial z}\right)+\boldsymbol{u}\cdot S_E
\end{aligned}
\end{equation}
\]
注意：为得到此方程，需要把方程33-35加起来（这里只考虑了压力和黏性力），四则运算一下，再补一个源项\(S_E\)，最后再乘以\(\boldsymbol{u}\)。
总能量方程与内能
总能量方程应为：
\[\begin{equation}
\rho\frac{De}{Dt}=-\nabla\cdot(p\boldsymbol{u})+\nabla\cdot(\boldsymbol{u}\cdot\boldsymbol{\tau})+\nabla\cdot(k\nabla T)+S_i++\boldsymbol{u}\cdot\boldsymbol{S}_E
\end{equation}
\]
其中：
- \(-\nabla\cdot(p\boldsymbol{u})\)是压力对流体做功的速率，\(\nabla\cdot(\boldsymbol{u}\cdot\boldsymbol{\tau})\)是黏性力对流体做功的贡献，二者共同组成描述表面力对流体微元作用的总功率，即方程42；
- \(\nabla \cdot (k\nabla T)\)是热传导的贡献，见方程43-45；
- 新加入的元素：\(S_i\)是单位体积的内部热源项，\(\boldsymbol{u}\cdot\boldsymbol{S}_E\)表示体力对流体微元的作用。
（PS：这个体力的项我没搞明白，它要不要乘以\(\rho\)？还有它到底该不该存在，为什么就突然引入动能方程了？虽然最后是要被抵消的）
而总动能方程简写为：
\[\begin{equation}
\rho\frac{DK}{Dt}=-\boldsymbol{u}\cdot\nabla p+\boldsymbol{u}\cdot\nabla\cdot\boldsymbol{\tau}+\boldsymbol{u}\cdot\boldsymbol{S}_E
\end{equation}
\]
总能量-总动能=总内能i，对每一小项分别计算：
-
压力项：
总能量方程中的为\(-\nabla\cdot (\rho \boldsymbol{u})\)，动能方程中的为\(-\boldsymbol{u}\cdot \nabla p\)，根据链式法则可得：\[\begin{equation} -\nabla\cdot(p\boldsymbol{u})-(-\boldsymbol{u}\cdot\nabla p)=-p\nabla\cdot\boldsymbol{u} \end{equation} \] -
黏性力项：
总能量方程中的为\(\nabla \cdot (\boldsymbol{u}\cdot \boldsymbol{\tau})\)，而动能方程中的为\(\boldsymbol{u} \cdot \nabla \cdot \boldsymbol{\tau}\)，二者相减的结果为：\[\begin{equation} \begin{aligned} (\nabla\boldsymbol{u}):\boldsymbol{\tau}&=\tau_{xx}\frac{\partial u}{\partial x}+\tau_{xy}\frac{\partial u}{\partial y}+\tau_{xz}\frac{\partial u}{\partial z}\\&+\tau_{yx}\frac{\partial v}{\partial x}+\tau_{yy}\frac{\partial v}{\partial y}+\tau_{yz}\frac{\partial v}{\partial z}\\&+\tau_{zx}\frac{\partial w}{\partial x}+\tau_{zy}\frac{\partial w}{\partial y}+\tau_{zz}\frac{\partial w}{\partial z} \end{aligned} \end{equation} \] -
热传导项和内部热源项保留。
将上述方程49和50合并进来，得到内能方程：
\[\begin{equation}
\rho\frac{Di}{Dt}=-p\nabla\cdot\boldsymbol{u}+(\nabla\boldsymbol{u}):\boldsymbol{\tau}+\nabla\cdot(k\nabla T)+S_i
\end{equation}
\]
展开之后的形式为：
\[\begin{equation}
\begin{aligned}
&\rho\frac{Di}{Dt} \\
&=- p\nabla\cdot\boldsymbol{u}+\nabla\cdot(k\nabla T)\\&+\tau_{xx}\frac{\partial u}{\partial x}+\tau_{yx}\frac{\partial u}{\partial y}+\tau_{zx}\frac{\partial u}{\partial z}\\&+\tau_{xy}\frac{\partial v}{\partial x}+\tau_{yy}\frac{\partial v}{\partial y}+\tau_{zy}\frac{\partial v}{\partial z} \\
&+\tau_{xz}\frac{\partial w}{\partial x}+\tau_{yz}\frac{\partial w}{\partial y}+\tau_{zz}\frac{\partial w}{\partial z}+S_{i}
\end{aligned}
\end{equation}
\]
进一步地，对于不可压缩流体，有：
- 密度恒定；
- 连续性方程，\(\nabla \cdot \boldsymbol{u}=0\);
- 有内能与温度的关系：\(i=cT\)，其中c为比热容。
那么此时由于连续性，不存在压力梯度，压力项消失，内能方程简化为：
\[\begin{equation}
\rho\frac{Di}{Dt}=(\nabla\boldsymbol{u}):\boldsymbol{\tau}+\nabla\cdot(k\nabla T)+S_i
\end{equation}
\]
\[\begin{equation}
\begin{aligned}
&\rho c\frac{DT}{Dt} \\
&=\nabla\cdot(k\nabla T)\\&+\tau_{xx}\frac{\partial u}{\partial x}+\tau_{yx}\frac{\partial u}{\partial y}+\tau_{zx}\frac{\partial u}{\partial z}\\&+\tau_{xy}\frac{\partial v}{\partial x}+\tau_{yy}\frac{\partial v}{\partial y}+\tau_{zy}\frac{\partial v}{\partial z} \\
&+\tau_{xz}\frac{\partial w}{\partial x}+\tau_{yz}\frac{\partial w}{\partial y}+\tau_{zz}\frac{\partial w}{\partial z}+S_{i}
\end{aligned}
\end{equation}
\]
N-S方程
一个流体元素的线性形变有9个分量：3个线性伸长分量（linear elongating components），与6个剪切线性形变分量（shearing linear doformation components）。
-
线性伸长分量：
\[\begin{equation} e_{xx}=\frac{\partial u}{\partial x},e_{yy}=\frac{\partial v}{\partial y},e_{zz}=\frac{\partial w}{\partial z} \end{equation} \] -
剪切线性形变分量：
\[\begin{equation} \begin{gathered} e_{xy}=e_{yx}=\frac{1}{2}\Big(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\Big) \\ e_{xz}=e_{zx}=\frac{1}{2}\Big(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\Big) \\ e_{yz}=e_{zy}=\frac{1}{2}\Big(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\Big) \end{gathered} \end{equation} \]
线性形变：考虑一个流体微元AB，正沿x方向移动，如下图所示。
注意，由于流体微元左边A和右边B处的速度不同（速度场的空间变化），故同样时间，运动距离不同，形成形变。
x方向的应变速率：
\[\begin{equation}
\begin{aligned}
&=\frac{d}{dt}\biggl(\frac{A^{\prime}B^{\prime}-AB}{AB}\biggr)=\frac{1}{\delta x}\biggl(\frac{d\delta x}{dt}\biggr) \\
&=\frac{1}{dt}\Bigg(\frac{\delta x+\left(u+\frac{\partial u}{\partial x}\delta x\right)dt -udt-\delta x}{\delta x}\Bigg)=\frac{\partial u}{\partial x}
\end{aligned}
\end{equation}
\]
等同于速度场的空间变化。
据此，我们将形变速率推广到三维，考虑一个边长为\(\delta x\)，\(\delta y\)和\(\delta z\)的立方体流体微元。其体积应变率是：
\[\begin{equation}
\begin{aligned}
\frac{1}{\delta V}\frac{d}{dt}(\delta V)& =\frac1{\delta x\delta y\delta z}\frac d{dt}(\delta x\delta y\delta z) \\
&=\frac1{\delta x}\frac d{dt}(\delta x)+\frac1{\delta y}\frac d{dt}(\delta y)+\frac1{\delta z}\frac d{dt}(\delta z) \\
\lim_{\delta V\text{一一}\to0}\left(\frac{1}{\delta V}\frac{d}{dt}(\delta V)\right)& =\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \\
&=Div(u)=\nabla\cdot \boldsymbol{u}
\end{aligned}
\end{equation}
\]
仍然等同于速度场的空间变化。也就是说，体积应变率由速度梯度张量的对角线给出，它告诉我们，在流动过程中，一个流体粒子/流体微元会伸展多少。
如图所示，一个二维的流体微元，由于速度梯度，产生剪切变形和角变形。
道理同上面一维形变，这里不再赘述。
注意，角O改变了\(90-(d\alpha + d\beta)\)度，因此其角变形（angular deformation，或单位变形 unit deformation）为\(d\alpha + d\beta\)度。
求角变形：
-
\(d\alpha\)：由于角度较小，故认为\(sin(d\alpha)=d\alpha\)，又有\(sin(d\alpha)=\frac{B B'}{O B'}\)。
\[\begin{equation} d\alpha =\frac{BB^{\prime}}{OB^{\prime}}=\frac{\left[\left(u+\frac{\partial u}{\partial y}\delta y \right)dt -udt \right]}{\delta y}=\frac{\partial u}{\partial y}dt \end{equation} \] -
\(d\beta\)：类似上面，有
\[\begin{equation} d\beta =\frac{AA'}{OA'}=\frac{\left[\left(v+\frac{\partial v}{\partial x}\delta x\right)dt -vdt \right]}{\delta x}=\frac{\partial v}{\partial x}dt \end{equation} \]
接下来求剪切应变速率：
\[\begin{equation}
\begin{aligned}&=\frac{(d\alpha+d\beta)}{dt}=\frac{\left(\frac{\partial u}{\partial y}dt +\frac{\partial v}{\partial x}dt\right)}{dt}\\&=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\quad\end{aligned}
\end{equation}
\]
虽然该结果仍然体现了速度场的空间分布，但注意：这里变成了x方向速度u沿着y方向上的空间变化，加上y方向速度v沿着x方向的空间变化，也就是不同方向上的速度分量梯度。如上图点O和点A的速度不同，可以很好地体现速度场在空间不同位置之间的变化。
根据公式57，我们已经知道体积应变速率为\(\nabla \cdot \boldsymbol{u}\)。
进一步地，粘性应力和形变速率是成正比的。对于牛顿流体，粘性应力可以用流体的动态粘性系数\(\mu\)，和应变速率描述。
注意到有3个正应力分量和6个剪切应力分量，类似前述。
-
正应力分量（normal stress components）：
\[\begin{equation} \begin{aligned}\tau_{xx}&=2\mu\frac{\partial u}{\partial x}+\lambda\nabla\cdot\mathbf{u}\\\tau_{yy}&=2\mu\frac{\partial v}{\partial y}+\lambda\nabla\cdot\mathbf{u}\\\tau_{zz}&=2\mu\frac{\partial w}{\partial z}+\lambda\nabla\cdot\mathbf{u}\end{aligned} \end{equation} \]注意，正应力包含两部分：一部分是流体在该方向上的变形速率，用\(\mu\)加权；另一部分是流体的体积变化速率（即散度），用第二粘性系数\(\lambda\)加权。
-
剪切应力分量（shear stress components）：
\[\begin{equation} \begin{align*} \tau_{xy} = \tau_{yx} &= \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\ \tau_{xz} = \tau_{zx} &= \mu \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\ \tau_{yz} = \tau_{zy} &= \mu \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \end{align*}\end{equation} \]它们表示不同方向上的速度分量梯度所产生的应力。其中，\(\tau_{yx}\)表示y方向的速度梯度在x方向上的变化导致的应力。
下一步就是正式进入到N-S方程了。
我们以x方向的流体输运作为例子。
\[\begin{equation}
\begin{align*}
\rho\frac{Du}{Dt}=-\frac{\partial p}{\partial x}+\frac{\partial}{\partial x}\Big[2\mu\frac{\partial u}{\partial x}+\lambda\nabla\cdot\boldsymbol{u}\Big]+\frac{\partial}{\partial y}\Big[\mu\Big(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\Big)\Big]\\+\frac{\partial}{\partial z}\Big[\mu\Big(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\Big)\Big]+S_{Mx}
\end{align*}
\end{equation}
\]
看到等式右边，压力梯度按从高到低的方向推动流体，加上x方向受到的1个正应力分量和2个剪切应力分量。等式左边是流体密度\(\times\)流体质点的加速度，总的来说是一个惯性项。
把方程63中的3个应力分量展开，可得：
\[\begin{equation}
\begin{aligned}
&\frac{\partial}{\partial x}\bigg[2\mu\frac{\partial u}{\partial x}+\lambda\nabla\cdot\boldsymbol{u}\bigg]+\frac{\partial}{\partial y}\bigg[\mu\bigg(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\bigg)\bigg]+\frac{\partial}{\partial z}\bigg[\mu\bigg(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\bigg)\bigg] \\
&=\frac{\partial}{\partial x}\Big(\mu\frac{\partial u}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(\mu\frac{\partial u}{\partial y}\Big)+\frac{\partial}{\partial z}\Big(\mu\frac{\partial u}{\partial z}\Big)+\Big[\frac{\partial}{\partial x}\Big(\mu\frac{\partial u}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(\mu\frac{\partial v}{\partial y}\Big)+\frac{\partial}{\partial z}\Big(\mu\frac{\partial w}{\partial z}\Big)\Big] \\
&+\frac\partial{\partial x}(\lambda\nabla\cdot\boldsymbol{u})=\nabla\cdot(\mu\nabla u)+S_{Mx}
\end{aligned}
\end{equation}
\]
具体是如何变换的呢？
- \(\frac\partial{\partial x}\left(\mu\frac{\partial u}{\partial x}\right)+\frac\partial{\partial y}\left(\mu\frac{\partial u}{\partial y}\right)+\frac\partial{\partial z}\left(\mu\frac{\partial u}{\partial z}\right)\)是梯度项，表示黏性对速度场的扩散作用，变成\(\nabla\cdot(\mu\nabla u)\)；
- \(\Big[\frac{\partial}{\partial x}\Big(\mu\frac{\partial u}{\partial x}\Big)+\frac{\partial}{\partial y}\Big(\mu\frac{\partial v}{\partial y}\Big)+\frac{\partial}{\partial z}\Big(\mu\frac{\partial w}{\partial z}\Big)\Big] +\frac\partial{\partial x}(\lambda\nabla\cdot\boldsymbol{u})\)包括散度项和体积黏度项，它们在黏性应力项中的贡献较少，因此隐藏在\(S_{Mx}\)动量源项中。
这样，我们可以得到一个更简洁漂亮的N-S动量方程形式：
\[\begin{equation}
\rho\frac{Du}{Dt}=-\frac{\partial p}{\partial x}+\nabla\cdot(\mu\nabla u)+S_{Mx}
\end{equation}
\]
类似地，对于y方向和z方向流动，有：
\[\begin{equation}\begin{aligned}
\rho\frac{Dv}{Dt}&=- \frac{\partial p}{\partial y}+\nabla\cdot(\mu\nabla v)+S_{My}\\\rho\frac{Dw}{Dt}&=- \frac{\partial p}{\partial z}+\nabla\cdot(\mu\nabla w)+S_{Mz}\end{aligned}
\end{equation}
\]
总结：有限体积法的基础
描述时间依赖的三维流体流动和传热的系统方程，守恒形式或散度形式为：
-
质量守恒：
\[\begin{equation} \frac{\partial\rho}{\partial x}+\nabla\cdot(\rho\boldsymbol{u})=0 \end{equation} \] -
动量守恒：
-
x方向
\[\begin{equation} \frac{\partial\rho u}{\partial t}+\nabla\cdot(\rho u\boldsymbol{u})=-\frac{\partial p}{\partial x}+\nabla\cdot(\mu\nabla u)+S_{Mx} \end{equation} \] -
y方向
\[\begin{equation} \frac{\partial\rho v}{\partial t}+\nabla\cdot(\rho v\boldsymbol{u})=-\frac{\partial p}{\partial y}+\nabla\cdot(\mu\nabla v)+S_{My} \end{equation} \] -
z方向
\[\begin{equation} \frac{\partial\rho w}{\partial t}+\nabla\cdot(\rho w\boldsymbol{u})=-\frac{\partial p}{\partial z}+\nabla\cdot(\mu\nabla w)+S_{Mz} \end{equation} \] -
能量
\[\begin{equation} \frac{\partial\rho i}{\partial t}+\nabla\cdot(\rho i\boldsymbol{u})=-p\nabla\boldsymbol{u}+\nabla\cdot(k\nabla\boldsymbol{u})+\Phi+S_i \end{equation} \] -
状态方程
\[\begin{equation} p=p(\rho,T)\quad and\quad i=i(\rho,T) \end{equation} \]
-
5个流体偏微分方程，2个额外的代数方程，再加上牛顿黏性模型；此时共7个未知数需求解，因此该系统在数学上是封闭的。
广义传输方程的微分形式和积分形式
引入一个通用变量\(\phi\)，可将所有流体流动方程（包括温度、污染物等标量的方程）写成如下形式：
\[\begin{equation}
\frac{\partial\rho\phi}{\partial t}+\nabla\cdot(\rho\phi\boldsymbol{u})=\nabla\cdot(\Gamma\nabla\phi)+S_\phi
\end{equation}
\]
各项的含义：
- 左1 \(\frac{\partial\rho\phi}{\partial t}\)：变化率项
- 左2 \(\nabla\cdot(\rho\phi\boldsymbol{u})\)：对流项
- 右1 \(\nabla\cdot(\Gamma\nabla\phi)\)：扩散项
- 右2 \(S_{\phi}\)：源项
通过把\(\phi\)设置为1，u，v，w，I，并选择合适的扩散系数\(\Gamma\)和源项，我们能获得一组特定方程。
对方程73，在三维控制体积内积分：
\[\begin{equation}
\int_{CV}\frac{\partial\rho\phi}{\partial t}dV+\int_{CV}\nabla\cdot(\rho\phi\boldsymbol{u})dV=\int_{CV}\nabla\cdot(\Gamma\nabla\phi)dV+\int_{CV}S_{\varphi}dV
\end{equation}
\]
使用高斯散度定理：
\[\begin{equation}
\begin{aligned}\int\limits_{CV}\nabla\boldsymbol{a}dV&=\int\limits_{A}\boldsymbol{n}\cdot\boldsymbol{a}dA\end{aligned}
\end{equation}
\]
注意，这里\(\boldsymbol{n}\cdot\boldsymbol{a}\)的物理含义为：正交于面积元dA的单位向量\(\boldsymbol{n}\)方向上的，向量\(\boldsymbol{a}\)的通量。
这样，方程74可写作：
\[\begin{equation}
\begin{aligned}
&\frac{\partial}{\partial t}\Bigg(\int\limits_{CV}\rho\phi dV\Bigg)+\int\limits_{A}\boldsymbol{n}\cdot(\rho\phi\boldsymbol{u})dA\\&=\int\limits_{A}\boldsymbol{n}\cdot(\Gamma\nabla\phi)dA+\int\limits_{CV}S_{\phi}dV\end{aligned}
\end{equation}
\]
各项的意义：
- 左1：控制体积内，流体微元物理性质\(\phi\)总量的变化率；
- 左2：由于对流引起的，流体微元物理性质\(\phi\)的净减少率；
- 右1：由于扩散引起的，流体微元物理性质\(\phi\)的净增加率；
- 右2：由于源项引起的，流体微元物理性质\(\phi\)的增加率。
考虑如下2种特殊情况：
-
对于稳态流动，\(\frac{\partial}{\partial t}=0\)，故式76的左1项消失，可写为：
\[\begin{equation} \int_A\boldsymbol{n}\cdot(\rho\phi\boldsymbol{u})dV=\int_A\boldsymbol{n}\cdot( \Gamma\nabla\phi)dV+\int_{CV}S_\phi dV \end{equation} \] -
对于依赖于时间的问题，还需要关于时间积分一次。最常用的形式为：
\[\begin{equation} \begin{gathered} \int\limits_{\Delta t}\frac{\partial}{\partial t}\Bigg(\int\limits_{CV}\rho\phi dV\Bigg)dt+\int\limits_{\Delta t}\int\limits_{A}\boldsymbol{n}\cdot d(\rho\phi\boldsymbol{u})A dt \\ =\int\limits_{\Delta t}\int\limits_{A}\boldsymbol{n}\cdot(\Gamma\nabla\phi)dA dt+\int\limits_{\Delta t}\int\limits_{CV}S_{\phi}dV dt \end{gathered} \end{equation} \]
