3.2.1 对流项离散格式的基本介绍 | 3.2.2 中心离散格式与一阶迎风格式（OpenFOAM理论笔记系列）

3.2 对流项的离散格式

3.2.1 对流项离散格式的基本介绍

在开始本节的讨论开始，笔者首先要说明，本节所介绍的内容相对于整个对流项离散格式的开发的历史可以说是极其简略的。可以这么说，整个计算流体动力学发展的历史有一半以上是由对流项离散格式写成的。本节的内容以Versteeg和Malalasekera的著作《An Introduction to Computational Fluid Dynamic: The Finite Volume Method Second edition》中的对流项离散格式介绍内容为框架，以对流项离散格式发展的历史过程为顺序，重点介绍几种重要的离散格式和他们的特点，分析OpenFOAM中对流项离散格式的植入过程，并最终为我们实际的模拟实践提供指导。
在第一章中，我们曾经使用高斯定理对对流项进行过有限体积离散：
$$\begin{gathered} {\int }_{V_P}\nabla \cdot (\rho \varphi \vec{U})dV=\sum _f\left[\overrightarrow{S_f}\cdot {\left(\rho \varphi \vec{U}\right)}_f\right]=\sum _f\left[\overrightarrow{S_f}\cdot (\rho \overrightarrow{U_f}){\varphi }_f\right] \\ =\sum _f(\mathbf{F}{\varphi }_f)\text{\hspace{1em}(1.28)} \end{gathered}$$
所谓对流项离散格式，主要解决的就是式(1.28)中面上值${\varphi }_f$如何计算的问题。根据我们上一节的讨论，使用与面上扩散系数相同的中心离散格式似乎是一种理所当然的选择，在下一节中我们会讨论中心离散格式在对流项离散格式中的表现，读者将会看到，中心离散格式在对流项离散过程中带来的一些问题以及这些问题如何推动对流项离散格式进一步发展。另一方面，纯对流的物理问题在自然界中很少存在，扩散在各种物理过程中广泛存在，包括N-S方程在内输运方程都属于对流-扩散方程，因此我们在讨论对流项离散格式的时候往往不单独讨论对流方程，而是考虑对流项离散格式对于对流-扩散方程计算效果的影响。另一方面，由于扩散项的离散格式基本没有什么选项可供选择(面上扩散系数采用中心离散格式，面法向梯度采用体心插值)，因此对流项离散格式基本上也决定了整个对流-扩散方程计算的质量。

为了方便叙述，在本节中我们以一维的对流扩散方程为例:
$$\frac{\mathrm{d}}{\mathrm{d}x}(\rho u\varphi )=\frac{\mathrm{d}}{\mathrm{d}x}\left(\Gamma \frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)\text{\hspace{1em}(3.12)}$$
在如下的网格系统中，上述方程可用之前提到的标准有限体积离散过程离散为：

• 图3.3 考虑相邻两个网格的一维网格系统¹
  

$$(\rho uA\varphi {)}_e-(\rho uA\varphi {)}_w={\left(\Gamma A\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}_e-{\left(\Gamma A\frac{\mathrm{d}\varphi }{\mathrm{d}x}\right)}_w\text{\hspace{1em}(3.13)}$$

速度的连续性方程可以写作：
$$(\rho uA{)}_e-(\rho uA{)}_w=0\text{\hspace{1em}(3.14)}$$
我们用之前对扩散项的处理方法处理式(3.13)右侧的扩散项并采用如下的符号表示式子中的对流通量和扩散通量：
$$\begin{array}{ll}{F}_w=(\rho u{)}_w& F_e=(\rho u{)}_e\\ D_w=\frac{{\Gamma }_w}{\delta x_{WP}}& D_e=\frac{{\Gamma }_e}{\delta x_{PE}}\end{array}\text{\hspace{1em}(3.15)}$$
最终我们可以得到整理后的离散对流-扩散方程以及连续性方程：
$$F_e{\varphi }_e-F_w{\varphi }_w=D_e\left({\varphi }_E-{\varphi }_P\right)-D_m\left({\varphi }_P-{\varphi }_W\right)\text{\hspace{1em}(3.16)}$$
$$F_e-F_w=0\text{\hspace{1em}(3.17)}$$

基于式(3.16)和式(3.17),我们用网格佩克莱数(cell Peclet Number)衡量一个控制体局部对流作用与扩散作用的强弱：
$$Pe=\frac{F}{D}=\frac{\rho u}{\Gamma /\delta x}\text{\hspace{1em}(3.12)}$$
佩克莱数中分子位置为对流通量，分母位置为扩散通量，$l$为特征长度，一般取网格的平均宽度。显然，贝克莱数越大，对流作用越强，贝克莱数越小，扩散作用越强。

为了准确计算对流-扩散方程，离散格式一般需要满足以下三个特性：守恒性(Conservativeness),有界性(Boundedness)，以及输运性(Transportiveness)。

守恒性:
守恒性的概念非常简单,即通过一个控制面的通量需要保持相等,具体来说,在控制体a和控制体b交接的控制面上,流出控制体a的通量与流入控制体b的通量应该保持一致。举例来说，中心离散格式对于扩散通量来说就是守恒的，考虑一维情况，在如图3.3所示的网格系统中，将流经1-4号控制体的扩散通量相加可得：

• 图3.4 中心离散格式的扩散通量守恒性²
  
  $$\begin{array}{l}\left[{\Gamma }_{e_1}\frac{\left({\varphi }_2-{\varphi }_1\right)}{\delta x}-q_A\right]+\left[{\Gamma }_{e_2}\frac{\left({\varphi }_3-{\varphi }_2\right)}{\delta x}-{\Gamma }_{m_2}\frac{\left({\varphi }_2-{\varphi }_1\right)}{\delta x}\right]\\ +\left[{\Gamma }_{e_3}\frac{\left({\varphi }_4-{\varphi }_3\right)}{\delta x}-{\Gamma }_{m_3}\frac{\left({\varphi }_3-{\varphi }_2\right)}{\delta x}\right]+\left[q_B-{\Gamma }_{w_4}\frac{\left({\varphi }_4-{\varphi }_3\right)}{\delta x}\right]\\ =q_B-q_A\end{array}\text{\hspace{1em}(3.13)}$$

显然，中心离散格式在离散通量时没有使流经系统的通量因为离散的原因增加或者减少。

有界性:
有界性是指模拟时，求得的物理量值应当被限制在一定的范围内，例如在多相流的计算中，任何一相的相分数都应该处于0-1之间，涉及能量的计算中，能量也不应该被求解出负值(例如在RANS湍流模型中，如果湍动能求解得到非物理的负值，将导致模拟直接发散)。为使得离散格式具有有界性，我们应该尽可能使得离散后的矩阵方程具有对角占优特征³,对角占优的矩阵系统意味着在缺乏源项的情况下，矩阵系统求解的值受到边界上的值的约束，其结果不会大于边界上的最大值，不会小于边界上的最小值。满足有界性的另一个有利条件是所有离散后的方程具有相同符号的系数(一般为都为正)，这意味着节点之间的变化是正向传递的，即一个节点值的升高会使得周围节点的值也升高，显然这种变化是符合物理规律的。如果离散格式不能保证有界性，那么可能导致模拟直接发散或者得到一个震荡的结果。

输运性：
输运性是指佩克莱数较大时(即对流作用占比较高时)，离散格式在求解面上值时应该更多地考虑来自流动上游的影响。由于扩散运动是没有特定方向的，因此之前我们使用中心离散格式进行扩散项的离散就比较合适，但是使用中心离散格式离散对流项就显得不太合适了，因为它无论流动朝向哪个方向，总是公平的将面上值处理为相邻两个控制体中心的插值。

3.2.2 中心离散格式与一阶迎风格式

在扩散项离散格式的讨论中，我们已经认识了中心离散格式，中心离散格式的思想非常简单，假设物理量呈线性变化，将面上的标量值处理为相邻两体心的线性插值(如图3.3的网格系统)：
$${\varphi }_w={\varphi }_W+\frac{\overline{wW}}{\overline{PW}}\left({\varphi }_P-{\varphi }_W\right)\text{\hspace{1em}(3.14)}$$
最早，人们也尝试使用中心离散格式对对流项进行离散，但数值模拟实践表明中心离散格式并不适合用作对流项离散格式。在这里我们仅给出结论，如果读者对其中的细节感兴趣可以阅读相关参考文献，文献将在本章结束的时候给出。我们按照之前叙述的守恒性，有界性，和输运性进行讨论：

• 守恒性：中心离散格式是守恒的。
• 有界性：当网格佩克莱数小于2时，中心离散格式才是有界准确稳定的。
• 输运性：中心离散格式没有根据流动方向进行插值权重的调整，不具有输运性。

需要指出的是，网格佩克莱数是流体物性，速度以及网格特征尺寸组成的参数，对于给定的物性，如果要降低佩克莱数使得其达到中心格式的稳定条件，就必须使得网格特征长度与速度的乘积尽可能小。这意味着在网格较粗糙的时候我们只能计算雷诺数非常低的流动，要计算高雷诺数流动，就必须大幅细化网格。在实际的模拟中，实用的网格密度条件下中心离散格式往往都是振荡不稳定的。正是基于这一点，中心离散格式才不适合作为对流项的离散格式。

一阶迎风格式是早期人们为克服中心离散格式的缺点发展的一种离散格式。中心离散格式不具有输运性，一阶迎风格式则是将输运性质发挥到极致。在图3.3的网格系统中，如果流动方向为由W至P，则迎风格式认为面上值完全由流动上游W点决定，如果流动方向由P至W,面上值完全由N点确定，即：
$${\varphi }_w=\{\begin{array}{rl}{\varphi }_W& F_w>0\\ {\varphi }_P& F_w<0\end{array}\text{\hspace{1em}(3.15)}$$
在式(3.15)中，$\mathbf{F}$如果大于0，则表示着速度与面矢量的点乘结果大于0，即速度方向从P指向N。

从守恒性，有界性和输运性这三个方面看，一阶迎风格式具有如下特性：

• 守恒性：一阶迎风格式格式是守恒的。
• 有界性：一阶迎风格式是有界无振荡的。
• 输运性：一阶迎风格式具有输运性质。

从上述三个性质看，迎风格式似乎是一种完美的对流项离散格式，但实际上却并不如此。首先，不需要进行泰勒展开，单单从式(3.15)就可以看出一阶迎风格式的精度只有一阶，这意味着如果将一阶迎风格式用于有限体积法，将导致整个方法的精度降为一阶。另外，一阶迎风格式最主要的缺点是当流动与网格线不对齐时(即流动方向不垂直于控制面)，会产生错误的结果。在这种情况下，迎风格式使得被输运的标量分布变得模糊，这种误差具有扩散的特征，一般称为假扩散和数值扩散。假扩散的大小与网格尺寸有关，网格越精细，假扩散的程度越低，结果越准确。

综上所述，中心格式与一阶迎风格式具有几乎完全相悖的性质。中心格式计算准确，但不稳定，容易出现振荡；一阶迎风格式数值耗散大，精度差，但是计算十分稳定；当网格质量非常好的时候，两种格式都可以达到较好的结果。为了调和中心格式与一阶迎风格式这种“相爱相杀”的关系，更多的离散格式被发展了出来。

系列说明：
接触有限体积法有一段时间了，也看了一些资料，但是有时候总觉得看过一遍之后什么也记不住。老话说得好，眼过千遍不如手过一遍，久而久之我就有了写一些比较像样子的笔记的想法。初学的时候曾经写过一本叫“OpenFOAM编程笔记：单相不可压缩流动”的册子，但当时基础太差（现在基础也不好），错误太多，倒不如推倒重来。
本系列将持续更新，欢迎各位挑错交流，挑错交流可以直接留言也可以联系邮箱cloudbird7@foxmail.com。等到预定内容全部写完后我将集结所有内容为一个独立的文件，开放下载。

---

1. Versteeg H K , Malalasekera W . An Introduction to Computational Fluid Dynamic: The Finite Volume Method Second edition [M]. Edinburgh Gate：Pearson Education, 2007. Figure 5.19. ↩︎

2. Versteeg H K , Malalasekera W . An Introduction to Computational Fluid Dynamic: The Finite Volume Method Second edition [M]. Edinburgh Gate：Pearson Education, 2007. Figure 5.7. ↩︎

3. Versteeg H K , Malalasekera W . An Introduction to Computational Fluid Dynamic: The Finite Volume Method Second edition [M]. Edinburgh Gate：Pearson Education, 2007. Figure 143. ↩︎