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}+\frac{\partial}{\partial y}+\frac{\partial}{\partial z} \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} \]
标签：方程,partial,推导,align,phi,体积,rho,delta,frac From： https://www.cnblogs.com/yukibvvd/p/18527577/3-limited-volume-method-derivative-equation-18