导读：有限体积法的工作原理。

一、求解对象

首先提一下我们的主角Navier-Stokes方程：

为了能够通过SIMPLE算法求解Navier-Stokes方程，我们需要将该方程改写为下列矩阵形式的方程：

其中为了系数矩阵，是每个网格未知的速度矢量，为右侧项。

一旦我们将Navier-Stokes方程写成线性化矩阵形式，我们就能够通过各种线性方程求解器来求解迭代过程中的速度场。

因此本文主讲解重点是将Navier-Stokes方程转换为矩阵形式，然后通过有限体积法求解。

二、有限体积法

为了方便演示，我们以稳定、不可压缩的Navier-Stokes方程为例：

采用三维多面体网格进行离散化。

在二阶有限体积法中，流动变量（、、)沿着网格线性变化，这些流动变量，包括压力、温度、速度等会被存储网格中心（）。

我们同样需要考虑相邻网格，通常每个网格都有个相邻网格，流动变量存储在这些网格中心（）。

这就是有限体积最原始的样子。

三、有限体积法工作原理

首先，我们需要沿着网格对Navier-Stokes方程进行积分：

通过积分的加法运算，意味着我们可以把上式括号里不同项分开，再对他们逐个积分:

在方程的左边，我们有对流项。方程右边有压力梯度、扩散项，还有一些源项，这里的源项是重力项。

3.1 源项

我们先来看一下源项，由于重力加速度是恒定，因此可以提到积分号外部，得到：

为网格的体积。回顾一下，前面提到的矩阵形式：

因此源项中，重力乘以网格体积完全不依赖速度，可以添加到方程右边（向量），而任何与速度成正比的东西都会被加到左边的系数矩阵中，所以常数源项，比如重力，很容易处理。

如果源项线性依赖速度要怎么处理？

假设源项被添加到Navier-Stokes方程(表示标量)：

同样需要将该源项在网格内进行积分：

（首先是标量，可以放到积分外；由于这是二阶的有限体积分法，速度在网格内的变化也是线性的，因此可以通过w网格中心的速度以及其他位置到网格中心点的距离获得速度值，从而对速度的体积积分进行展开；最后在将括号里面的加法进行展开，第二项的积分恒等于0，这是网格中心的定义）

这就是线性源项:

重新回到矩阵形式，如果我们把源项挪到左边可以添加到矩阵中，或者添加到矩阵中。

因此对于源项我们有两种不同的处理方法：将它作为隐式项添加到左边，或者让将它作为显式项加到右边。

这样处理会有什么区别呢？

如果是隐式处理，将矩阵 相加，则会得到下列矩阵：

这是一个对角矩阵，因此我们现在关注的是网格的中心，矩阵内非对角项则是与相邻的网格相关。

如果是显式处理，则会变为：

在CFD求解中，求解器通常会动态选择使用显式处理或隐式处理，如果S项是负的，而我们想要对角线项的贡献是正的，则应该选择隐式处理，可以增强对角线优势，使矩阵更加稳定。

如果S项是负的，这就减少了对角线的优势，所以应该选择显式处理。

因此源项的处理取决于我们是否要改善求解的对角线优势。

如果是非线性源项呢？

假如我们面对的是下列源项：

我们通常会尝试线性化源项：

对源项的积分，由于上一次迭代的速度是已知的，因此可以得到：

这是一种混合的隐式+显式处理方法。同理，我们可以采用同样的方式处理速度的立方项。

3.2 对流项与扩散项

对流项与扩散项具有散度计算（），因此他们更难处理。

由于篇幅有限，这里主要讲解对流项。

我们一般是用散度定理处理对流与扩散项。

散度定理指的是对任意向量的散度积分等价于向量的曲面积分：

因此我们可以通过散度定理修改对流项及扩散项，算出积分：

速度，单位法向量和表面的点积就是流出表面的体积流量：

括号外面的速度（)就是我们需要求解的量。

接下来我们需要分解曲面积分：

速度沿着面的变化是线性的，因此我们可以取面中心的速度作为速度面积分的近似【1】。

这样一来我们就消除了积分，只需要做求和。

但是这里其实是未知的，我们只知道网格中心的变量值。因此我们需要进行插值。

在CFD中提供多种插值方式来求解面速度，包括：迎风、二阶/线性迎风、中心差分、QUICK。

小结一下，对流项我们可以写为：

因此我们可以将对流项加入矩阵形式当中：

如下所示，对流项会带来对角线项和非对角项，这是相邻网格连通性的体现：\mathcal{M}=\left(

\right)。

4 小结

（1）Navier-Stokes方程中的每一项都能够一一积分：

→有限体积法基础文档每一项对矩阵的贡献不同。我们需要做的就是将这些贡献加起来组成完整的矩阵形式（2）这些贡献会矩阵带来对角线项和非对角线项