CFD理论 | 聊一聊有限体积法！

weixin

一、求解对象
首先，提一下我们的主角Navier-Stokes方程：

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

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

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

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

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

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

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

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

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

三、有限体积法工作原理
首先，我们需要沿着网格P 对Navier-Stokes方程进行积分：

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

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

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

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

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

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

假设源项SU 被添加到Navier-Stokes方程（S 表示标量）：

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

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

这就是线性源项：

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

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

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

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

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

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

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

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

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

如果是非线性源项呢？

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

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

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

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

2. 对流项与扩散项

对流项与扩散项具有散度计算 (▽·)，因此他们更难处理。

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

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

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

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

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

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

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

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

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

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

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

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

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

如下所示，对流项会带来对角线项和非对角项，这是相邻网格连通性的体现：\mathcal{M}=\left(M11  M12  M13  ...  0  M21  M22  0  ...  0 0 M23  M33  ...  0  ┆ ┆ ┆ ┆ ┆ 0 0 0  ...  Mmm\right)。

四、小结
Navier-Stokes方程中的每一项都能够一一积分：

每一项对矩阵的贡献不同。我们需要做的就是将这些贡献加起来组成完整的矩阵形式MU=B 这些贡献会矩阵带来对角线项和非对角线项，如：

参考：
[1] Jasak H . Error Analysis and Estimation for the Finite Volume Method With Applications to Fluid Flows. 1996.Fluid Mechanics101

来源：BB学长微信公众号（ID：fj_chenzhibin），作者：CFD控。