3 有限体积法：推导方程

3 有限体积法：推导方程

基本原理和目标

（注意：这一节看不懂没关系，在后面的推导中会慢慢用到）

  • 质量、动量和能量的守恒

    1. 流体的质量守恒
    2. 动量改变的速度 = 一个流体粒子上受到的力的总和（牛顿第二定律）
    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\)，也就是代表了该流体微元的总质量。

然后求质量关于时间的一阶导，就是质量的净增加速度了。

穿过流体元素边界，流入流体元素的净质量流动速率为：

​image​

其中，约定流入流体微元的质量流动标正号，流出的标负号。

注意这里是如何运用泰勒级数展开的（以x方向，即图中横着的方向为例）：

  1. 在点x处的质量流动速率密度（即单位面积的质量流率）为\(\rho u\)；

  2. 使用一阶泰勒级数展开：

    \[\begin{align}f(x+\delta x)=f(x)+\frac{df(x)}{dx} \delta x\end{align} \]

  3. 在\(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\)（即流出处）同理，此处不再展示。

  4. 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} \]

其具体推导过程为：

  1. 原实质导数参见方程17：\(\frac{D\phi}{Dt}=\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\nabla\phi\)​

  2. 首先，推广到密度和物理量一起，即考虑\(\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} \]

  3. 下一步，考虑\(\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} \]

  4. 将过程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} \]

  5. 根据连续性方程，质量守恒要求：

    \[\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} \]

  6. 注意到方程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\)速度关于时间的导数。

作用于流体元上的力：

  • 表面力：作用于流体表面，由于流体元与其周围流体或边界的接触产生。

    1. 压力力（pressure forces）

      • 由于流体中的压力梯度作用于流体元表面产生的力；
      • 垂直于流体元表面，指向外界（不过这里好像也不一定，要看具体流动是出微元还是进）；
      • 导致流体加减速，是推动流体流动的重要因素。

    2. 黏性力（viscous forces）

      • 由于流体的黏性（流体分子之间的内摩擦），在流体元表面产生的剪切应力；
      • 在流体层之间起阻力作用，特别是在层流和接触边界处；
      • 造成速度梯度（摩擦-相对速度），使得流体从高速区域向低速区域传递动量，最终影响整个流体的流动特性；
      • 补充：摩擦力是分子间的吸引力（如范德华力）、原子间的电磁排斥力、原子间电子交换或共享形成弱键合（如氢键）等等电磁力在宏观下的体现。

  • 体力（body forces，不是HP）：外界力场直接作用于流体的每个体积单元内部。

    1. 重力，影响流体压力分布、密度梯度等。
    2. 离心力（centrifugal force）：在旋转系统（如地球自转引起的大气、海洋环流）中，由于旋转参考系下的惯性效应产生的力，指向远离旋转轴的方向。
    3. 科里奥利力（Coriolis force）：旋转参考系下的虚拟力，在地球表面上主要由地球自转引起。在北半球，流体运动向右偏；在南半球，流体运动向左偏。
    4. 电磁力（electromagnetic force）：带电流体元在电场和磁场中所受的洛伦兹力，在工程力学中常用于控制电解液、等离子体的运动。

一般地，在动量方程中，我们将表面力（压力和黏性力）视为主要项并将它们分开考虑，而将体力视为源项。

通过压力和9个粘性应力分量，我们定义了流体元素的应力状态。

​image​

该图是为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方向移动，如下图所示。

​image​

注意，由于流体微元左边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} \]

仍然等同于速度场的空间变化。也就是说，体积应变率由速度梯度张量的对角线给出，它告诉我们，在流动过程中，一个流体粒子/流体微元会伸展多少。​​

​2024_11_11_20_55_11​

如图所示，一个二维的流体微元，由于速度梯度，产生剪切变形和角变形。

道理同上面一维形变，这里不再赘述。

注意，角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} \]

  • 动量守恒：

    1. 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} \]

    2. 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} \]

    3. 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} \]

    4. 能量

      \[\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} \]

    5. 状态方程

      \[\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} \]

‍

posted @ 2024-11-05 19:30  灰鲤  阅读(41)  评论(0)  编辑  收藏  举报