3.1 扩散项的离散格式（OpenFOAM理论笔记系列）

第三章 离散格式

3.1 扩散项的离散格式

3.1.1 面法向梯度的计算

有限体积法中，物理量都存储在控制体体心上，而由于我们在进行有限体积离散时使用了高斯定理，因此在最终组建代数方程时，我们往往都用的是控制面上的值。在这里，就牵扯倒一个如何根据体心值去计算面上值的过程，这就是离散格式要解决的问题。对于一个标量$\varphi$,按照1.3.2节的方法，我们有：
$${\int }_{V_P}\nabla \cdot (\Gamma \nabla \varphi )dV=\sum _f\left[\overrightarrow{S_f}\cdot {\left(\Gamma \nabla \varphi \right)}_f\right]=\sum _f\left[\overrightarrow{S_f}\cdot {\Gamma }_f{\left(\nabla \varphi \right)}_f\right]\text{\hspace{1em}(1.29)}$$
显然，扩散项的离散最终落实倒了面扩散系数$\Gamma$和面上梯度${\left(\nabla \varphi \right)}_f$的离散上。对于面扩散系数的离散格式，我们往往采用中心格式，即使用当前控制体体心P与相邻控制体体心N的值进行线性插值：

• 图3.1 面上扩散系数的离散
  
  $${\Gamma }_f={\Gamma }_P+\frac{\overline{Pf}}{\overline{PN}}\left({\Gamma }_N-{\Gamma }_P\right)\text{\hspace{1em}(3.1)}$$

相对复杂一点的是面上梯度的计算，因为我们并不知道面上梯度${\left(\nabla \varphi \right)}_f$的方向。在这里，我们首先将面上的梯度${\left(\nabla \varphi \right)}_f$分解为两个方向，即面的法向分量${\left(\nabla \varphi \right)}_f^n=(\vec{n}\cdot {\left(\nabla \varphi \right)}_f)\vec{n}$和切向分量${\left(\nabla \varphi \right)}_f^t=(\vec{t}\cdot {\left(\nabla \varphi \right)}_f)\vec{t}$，此时扩散项中的 $\overrightarrow{S_f}\cdot {\left(\nabla \varphi \right)}_f$可以写为（明确面矢量的方向是垂直于控制面的）：
$$\begin{array}{l}\overrightarrow{S_f}\cdot {\left(\nabla \varphi \right)}_f=\overrightarrow{S_f}\cdot (\vec{n}\cdot {\left(\nabla \varphi \right)}_f)\vec{n}+\overrightarrow{S_f}\cdot (\vec{t}\cdot {\left(\nabla \varphi \right)}_f)\vec{t}\\ =(\vec{n}\cdot {\left(\nabla \varphi \right)}_f)\overrightarrow{S_f}\cdot \vec{n}\\ =(\vec{n}\cdot {\left(\nabla \varphi \right)}_f)|\overrightarrow{S_f}|\end{array}\text{\hspace{1em}(3.2)}$$
式(3.2)中$\vec{n}\cdot {\left(\nabla \varphi \right)}_f$即我们常说的面法向梯度，其是面上梯度$\left(\nabla \varphi \right)$在网格控制面法向上的分量。对于如图3.1所示的正交网格系统(即控制体中心连线垂直于网格控制面)，面法向梯度显然可以按照下式进行计算：
$$\vec{n}\cdot {\left(\nabla \varphi \right)}_f={\left[{\left(\nabla \varphi \right)}_f\right]}_{PN}=\frac{{\varphi }_N-{\varphi }_P}{\overline{PN}}=\frac{{\varphi }_N-{\varphi }_P}{|\vec{d}|}\text{\hspace{1em}(3.3)}$$
式{3.3}中$\vec{d}$指向量$\overrightarrow{PN}$。式3.3即我们计算面法向梯度的基本公式，但是这个式子在网格非正交时是需要修正的。在非正交情况下，我们可以将面法向矢量$\vec{n}$分解为两个矢量$\vec{\Delta }$和$\vec{k}$，其中$\vec{\Delta }$方向与$\vec{d}$平行，即：
$$\vec{n}\cdot {\left(\nabla \varphi \right)}_f=\vec{\Delta }\cdot {\left(\nabla \varphi \right)}_f+\vec{k}\cdot {\left(\nabla \varphi \right)}_f\text{\hspace{1em}(3.4)}$$
其中$\vec{\Delta }\cdot {\left(\nabla \varphi \right)}_f$称为正交项，显然其为面上梯度沿着$\vec{d}$方向的分量。注意到在式(3.3)中，$\vec{n}$是与$\vec{d}$同方向的单位矢量，且$\vec{\Delta }$与式(3.3)中的$\vec{n}$同方向，因此$\vec{\Delta }\cdot {\left(\nabla \varphi \right)}_f=|\vec{\Delta }|\frac{{\varphi }_N-{\varphi }_P}{|\vec{d}|}$。因此式(3.4)可转化为：
$$\vec{n}\cdot {\left(\nabla \varphi \right)}_f=|\vec{\Delta }|\frac{{\varphi }_N-{\varphi }_P}{|\vec{d}|}+\vec{k}\cdot {\left(\nabla \varphi \right)}_f\text{\hspace{1em}(3.5)}$$
显然，想要求得 $\vec{n}\cdot {\left(\nabla \varphi \right)}_f$，我们需要知道 $|\vec{\Delta }|$以及$\vec{k}\cdot {\left(\nabla \varphi \right)}_f$的值。考虑到 $\vec{n}=\vec{\Delta }+\vec{k}$且$\vec{n}$为面的单位法向矢量，在网格结构确定的情况下是已知的，因此我们只需要知道 $\vec{\Delta }$即可求得 $|\vec{\Delta }|$和 $\vec{k}$的值。那么 $\vec{k}\cdot {\left(\nabla \varphi \right)}_f$的值怎么求呢？

事实上，我们还有另一种求解${\left(\nabla \varphi \right)}_f$的方法，由第一章式(1.3)可得：
$${\int }_{V_P}\nabla \varphi dV={\left(\nabla \varphi \right)}_P{V}_p\text{\hspace{1em}(3.6)}$$
同时，类比式(1.20)与式(1.2５)可得：
$${\int }_{V_P}\nabla \varphi dV=\sum _f\left({\varphi }_f\overrightarrow{S_f}\right)\text{\hspace{1em}(3.7)}$$
联立式(3.6)和式(3.7)可得：
$${\left(\nabla \varphi \right)}_P=\frac{1}{V_P}\sum _f\left({\varphi }_f\overrightarrow{S_f}\right)\text{\hspace{1em}(3.8)}$$
而面上梯度${\left(\nabla \varphi \right)}_f$可以仿照式(3.1)进行插值：
$${\left(\nabla \varphi \right)}_f={\left(\nabla \varphi \right)}_P+\frac{\overline{Pf}}{\overline{PN}}\left({\left(\nabla \varphi \right)}_N-{\left(\nabla \varphi \right)}_P\right)\text{\hspace{1em}(3.9)}$$
那么既然我们可以使用式(3.9)来计算${\left(\nabla \varphi \right)}_f$，那为什么之前还要使用式(3.3)来计算面法向梯度呢？这主要出于节省计算资源和提高计算精度的考虑。首先，用式(3.9)乘以单位面法向量显然要比直接使用式(3.3)计算面法向梯度需要更多的计算量。更重要的是，如果使用泰勒展开的方法展开式(3.3)和式(3.9)，将会发现式(3.9)的计算误差是式(3.3)的4倍¹，即使式(3.9)多出来的计算资源消耗还是可以克服的，但这个计算误差的放大则是无法接受的，因此我们不选用式(3.9)来计算面法向梯度。

然而在网格非正交的情况下，我们无法简单使用(3.3)来计算面法向梯度，我们必须使用式(3.9)来计算式(3.4)中的非正交部。最终，在非正交情况下，式(3.4)最终转化为：

$$\begin{array}{l}\vec{n}\cdot {\left(\nabla \varphi \right)}_f=\vec{\Delta }\cdot {\left(\nabla \varphi \right)}_f+\vec{k}\cdot {\left(\nabla \varphi \right)}_f=|\vec{\Delta }|\frac{{\varphi }_N-{\varphi }_P}{|\vec{d}|}+\left(\vec{n}-\vec{\Delta }\right)\cdot {\left(\nabla \varphi \right)}_f^{\ast }\\ \text{其中:}\\ {\left(\nabla \varphi \right)}_f^{\ast }={\left(\nabla \varphi \right)}_P+\frac{\overline{Pf}}{\overline{PN}}\left({\left(\nabla \varphi \right)}_N-{\left(\nabla \varphi \right)}_P\right)\\ {\left(\nabla \varphi \right)}_P=\frac{1}{V_P}{\left[{\sum }_f\left({\varphi }_f\overrightarrow{S_f}\right)\right]}_P\\ {\left(\nabla \varphi \right)}_N=\frac{1}{V_N}{\left[{\sum }_f\left({\varphi }_f\overrightarrow{S_f}\right)\right]}_N\end{array}\text{\hspace{1em}(3.10)}$$
显然，我们只需要知道$\vec{\Delta }$就可以根据式(3.10)计算出网格非正交情况下的面法向梯度了。根据不同的$\vec{\Delta }$设计方式，存在多种非正交修正方法：这里我们介绍其中的一种最小值修正方法(Minimum correction approach)

• 图3.2 最小值修正方法的法向量分解
  

如图3.2所示，在最小值修正方法中，$\vec{\Delta }$和$\vec{k}$相互垂直，此时根据向量之间的关系有：
$$\vec{\Delta }=\frac{\vec{d}\cdot \vec{n}}{\vec{d}\cdot \vec{d}}\vec{d}\text{\hspace{1em}(3.11)}$$
在最小值修正方法中，$\vec{k}$的大小取得最小，可以确保相对不精确的非正交项在整个面法向梯度中占比最小。

除了最小值修正方法，还存在着正交修正方法(Orthogonal Correction approach)和过松弛修正方法(Over-relaxed approach),感兴趣的读者可以阅读相关文献²。目前笔者还不确定OF中使用了哪种修正方法，但有可能是过松弛修正法，感兴趣的同学可以自行探索。

3.2.2 OpenFOAM中的扩散项离散格式设置

OpenFoam中设置扩散项的格式为:

Gauss <插值格式> <面法向梯度格式>

举例来说，Gauss linear orthogonal表示对于扩散系数的面上插值，使用线性插值格式(linear)计算，对于面法向梯度，采用orthogonal格式进行离散。

OpenFOAM中有很多的插值格式，但在扩散项的离散中，我们一般使用linear格式进行离散。在涉及多相流计算时，我们也有时使用调和平均(harmonic)来计算面上扩散系数³。

面法向梯度格式主要有三种：orthogonal,corrected和uncorrected。

orthogonal格式是指不对非正交做任何修正，完全按照式(3.3)的方法进行计算，显然，只有在网格正交性非常好的情况下才可以使用这种格式。

但实际中不存在完全正交的网格，因此必须对非正交网格进行修正，即使用corrected格式。同时，corrected格式对隐性离散的正交计算部分以及显性离散的非正交部分应用了亚松弛因子$co{s}^{-1}\alpha$增加对角占优特性，其中$\alpha$是两个相邻控制体体心连线形成的向量与面法向量之间的夹角。

如果想要在正交性非常好的网格中不进行非正交修正并且还使用亚松弛因子，可以使用uncorrected格式，此时其相当于orthogonal格式附加$co{s}^{-1}\alpha$亚松弛因子，也相当于不进行非正交修正的的corrected格式。

当$\alpha$过大时(例如大于70°)，网格正交性非常差，非正交修正会使得计算变得非常不稳定，此时可以使用limited关键词对非正交项进行限制。limited关键词后需要跟一个限制系数$\psi$，根据OpenFOAM用户手册的说明，修正值的取值为：
$$\psi =\{\begin{array}{rl}0& \text{相当于}uncorrected\\ 0.333& \text{非正交修正部分}\le 0.5\times \text{正交修正部分}\\ 0.5& 非正交修正部分\le \text{正交修正部分}\\ 1& \text{相当于}corrected\end{array}$$
$\psi$越小，非正交项占比越小，计算越稳定，但精度低，适合正交性非常不好的情况，$\psi$越大，计算越不稳定，但精度高，适合正交性相对较好的情况。

我们在使用limited关键词时可以配合corrected关键词使用，例如corrected limited 0.5，这是较为常见的用法。同时，我们还可以直接使用limited关键词，例如limited 0.5，其含义与corrected limited 0.5一致，并且limited 0与uncorrected等价，limited 1与corrected等价。

通常，当$\alpha$小于70°时，推荐使用corrected格式，当$\alpha$在70°-85°时，推荐使用limited 0.55，当$\alpha$大于85°时，建议重新设计网格，如果你仍然想要使用当前的网格系统，你可以尝试使用limited 0.333并增加压力速度耦合算法中的非正交修正次数。corrected格式一般使用在网格质量极差的情况下，此时的计算精度非常低，但可以得到非常好的稳定性。另外，除非网格质量极好，否则不推荐使用orthogonal格式。

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

---

1. Jasak, Hrvoje. Error analysis and estimation for the finite volume method with applications to fluid flows.[D]. London：Imperial College London . 1996 . 84 ↩︎

2. Jasak, Hrvoje. Error analysis and estimation for the finite volume method with applications to fluid flows.[D]. London：Imperial College London . 1996 . 85-86 ↩︎

3. 陶文铨. 数值传热学[M]. 西安交通大学出版社, 2001.80-83 ↩︎