FVM in CFD 学习笔记_第9章_梯度计算

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 9 Gradient Computation

FVM in CFD 学习笔记_第9章_梯度计算

在CFD的FVM的离散过程当中，在单元形心和面形心处变量 $\phi$的梯度是经常要用到的物理量，那么如何由单元形心处的变量 $\phi$来获取单元形心和面形心处的变量 $\phi$的梯度 $\nabla \phi$呢？本节便讲讲FVM in CFD中梯度的计算方法。

1 笛卡尔网格上梯度的计算

因为笛卡尔网格横平竖直，有很好的正交特性，因此梯度的计算变得十分简单快捷了，对于1维问题，且均匀网格，则面 $e$上的变量梯度可直接得出：

$\left(\frac{\partial\phi}{\partial x}\right)_e=\frac{\phi_E-\phi_C}{x_E-x_C}=\frac{\phi_E-\phi_C}{\delta x_e}$
单元形心 $C$处的变量梯度也可得出
$\left(\frac{\partial\phi}{\partial x}\right)_C=\frac{\phi_e-\phi_w}{x_e-x_w}=\frac{\displaystyle \frac{\phi_E+\phi_C}{2}-\frac{\phi_C+\phi_W}{2}}{\Delta x_C}=\frac{\phi_E-\phi_W}{2 \Delta x_C}$
这个常叫作“中心差分”，对于2维状况，与此相似，可得

$\left(\frac{\partial\phi}{\partial x}\right)_C=\frac{\phi_E-\phi_W}{x_E-x_W},\quad \left(\frac{\partial\phi}{\partial y}\right)_C=\frac{\phi_N-\phi_S}{x_N-x_S}$
当处理非正交网格或非结构网格时，状况就复杂得多了，这就须要用到更加通用方法，即Green-Gauss梯度法和最小二乘梯度法等。

2 Green-Gauss Gradient（格林-高斯梯度）

这是计算梯度方法中应用最广的一个，即，单元形心处变量的梯度能够由面形心处的变量值与面积矢量复合后相加和除以单元体积来获取，用到的数学定理是Green-Gauss或梯度定理，即
$\int_V\nabla\phi dV=\oint_{\partial V}\phi d\vec S \Rightarrow \overline{\nabla\phi} V =\oint_{\partial V}\phi d\vec S \Rightarrow \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_f \vec S_f$
即
$\nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_f \vec S_f$
其中 $V_C$为单元体积， $f$表明单元的某个面，而 $\vec S_f$为该面的面积矢量，面形心的变量值 $\phi_f$是未知的，仍旧须要计算出来，否则上面这个公式是用不了的。 $\phi_f$的计算方法有两种，一是用紧致框架（compact stencil）由面左右单元（所属单元和邻居单元，即面两侧单元）形心值来计算，二是用扩展框架（extended stencil）先用面的角点周围单元上的值来获得面的角点值，而后再由角点值来平均获得面形心的值。

方法1：紧致框架

对于上图(a)和(b)所示的2维和3维状况，一种最简单的方法是直接由面两侧单元形心值平均来得到面上变量的值，即
$\phi_f=g_C\phi_C+(1-g_C)\phi_F$
其中 $g_C$为几何权重系数，等于
$g_C=\frac{||\vec r_F-\vec r_f||}{||\vec r_F-\vec r_C||}=\frac{d_{Ff}}{d_{FC}}$
其中 $\vec r$为距离矢量，而 $d$为两点之间的距离值，若该面刚好位于两单元中心的中间位置，则有
$\phi_f=\frac{\phi_C+\phi_F}{2}$
这个方法很是简便易行，然而遗憾的是，只有当线段 $CF$和面 $\vec S_f$的交点与面 $\vec S_f$的中心重合时（即上图（a）（b）状况），该方法才能保证2阶精度，而大多数状况下，直接由上式算得的梯度精度是无法保证的。

然而，对于一般的非正交网格和非结构网格来讲，是没有办法来保证这个条件的，好比上图中的（c）（d）状况，网格的偏斜（非正交，skewness（non-conjunctionality））使得线段 $CF$和面 $\vec S_f$交点为 $f'$，与面形心 $f$是不重合的。此时，须要把插值获得的 $\phi_{f'}$修正以获得 $\phi_f$，修正方程为
$\phi_f=\phi_{f'}+correction=\phi_{f'}+(\nabla\phi)_{f'}\cdot(\vec r_f-\vec r_{f'})$
即，用梯度 $(\nabla\phi)_{f'}$来修正，修正也能够展开为以下形式：
$\phi_f=g_C\{\phi_C+(\nabla\phi)_C\cdot(\vec r_f-\vec r_C)\}+(1-g_C)\{ \phi_F+(\nabla\phi)_F\cdot(\vec r_f-\vec r_F) \}\\ \quad \\ =\phi_{f'}+g_C(\nabla\phi)_C\cdot(\vec r_f-\vec r_C)+(1-g_C)(\nabla\phi)_F\cdot(\vec r_f-\vec r_F)$
即， $\{\phi_C+(\nabla\phi)_C\cdot(\vec r_f-\vec r_C)\}$是把 $C$处的值修正到 $f$处，而 $\{ \phi_F+(\nabla\phi)_F\cdot(\vec r_f-\vec r_F) \}$是把 $F$处的值修正到 $f$处，而后再作加权插值处理，就获得了 $\phi_f$。

不难发现， $\phi_f$的修正计算要用到 $(\nabla\phi)_C$和 $(\nabla\phi)_F$，而 $(\nabla\phi)_C$和 $(\nabla\phi)_F$则是用 $\phi_f$算得的，也就是说，这里 $\phi_f$的计算是不能一蹴而就的，而是一个迭代计算的过程，可是过多的迭代反而会形成解的振荡，通常作两次迭代就行了。

另外， $g_C$是与点 $f'$的位置密切相关的，有三种选择方式：

选择1

点 $f'$就选择在线段 $CF$和面 $\vec S_f$的交点处，以 $\vec n$表明面积单位矢量，即， $\vec n=\vec S_f/||\vec S_f||$，以 $\vec e$表明沿着 $CF$的单位矢量，即 $\vec e=\overrightarrow{CF}/||\overrightarrow{CF}||$，则 $f'$的位置可由几何关系算出，为
$\vec r_{f'}=\frac{(\vec r_f - \vec r_C) \cdot \vec n}{\vec e \cdot \vec n}\vec e+\vec r_C=\frac{\vec r_f \cdot \vec n}{\vec e \cdot \vec n}\vec e$
其中， $(\vec r_f - \vec r_C) \cdot \vec n$为点 $C$到面 $\vec S_f$的垂直距离向量，分母 $\vec e \cdot \vec n$为该垂直向量与 $\overrightarrow{CF}$夹角的余弦值 $cos\theta$，因而二者相除便获得了 $C$到 $f'$的向量 $\overrightarrow{Cf'}$，再与 $\vec r_C$相加即是点 $f'$的位置矢量。
获得 $f'$的位置后，能够直接算出 $g_C$的值
$g_C=\frac{||\vec r_F - \vec r_{f'}||}{||\vec r_F - \vec r_{C}||}=\frac{d_{Ff'}}{d_{FC}}$
计算流程以下：

1. 计算 $\phi_{f'}$使用 $\phi_{f'}=g_C\phi_C+(1-g_C)\phi_F$
2. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f$
   接下来用下面步骤来修正梯度场
3. 更新 $\phi_{f}$使用 $\phi_f=\phi_{f'}+g_C(\nabla\phi)_C\cdot(\vec r_f-\vec r_C)+(1-g_C)(\nabla\phi)_F\cdot(\vec r_f-\vec r_F)$
4. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$
5. 返回步骤3再迭代一次。

选择2

点 $f'$就选择在线段 $CF$的中点，至关于 $g_C=0.5$，这就使得计算简单多了，其计算流程以下：

1. 计算 $\phi_{f'}$使用 $\displaystyle \phi_{f'}=\frac{\phi_C+\phi_F}{2}$
2. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f$
   接下来用下面步骤来修正梯度场
3. 更新 $\phi_{f}$使用 $\displaystyle \phi_f=\phi_{f'}+\frac{(\nabla\phi)_C+(\nabla\phi)_F}{2}\cdot\left(\vec r_f-\frac{\vec r_C+\vec r_F}{2}\right)$
4. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$
5. 返回步骤3再迭代一次。

选择3

点 $f'$的位置选择要保证距离 $ff'$是最短的，即， $ff'$是垂直于 $CF$的，这使得第1步迭代计算变得更加准确。 $f'$的位置计算很简单，直接把向量 $\overrightarrow {Cf}$投影到 $\overrightarrow {CF}$上就妥了，即
$\vec r_{f'}=\vec r_C+\frac{\vec r_{Cf} \cdot \vec r_{CF}}{\vec r_{CF} \cdot \vec r_{CF}}(\vec r_F-\vec r_C)$
其计算流程以下：

1. 计算 $\vec r_{f'}$使用 $\displaystyle\vec r_{f'}=\vec r_C+\frac{\vec r_{Cf} \cdot \vec r_{CF}}{\vec r_{CF} \cdot \vec r_{CF}}(\vec r_F-\vec r_C)$
2. 计算 $g_C$使用 $\displaystyle g_C=\frac{||\vec r_F - \vec r_{f'}||}{||\vec r_F - \vec r_C||}$
3. 计算 $\phi_{f'}$使用 $\phi_{f'}=g_C\phi_C+(1-g_C)\phi_F$
4. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f$
   接下来用下面步骤来修正梯度场
5. 计算 $\nabla\phi_{f'}$使用 $\nabla\phi_{f'}=g_C\nabla\phi_C+(1-g_C)\nabla\phi_F$
6. 更新 $\phi_{f}$使用 $\phi_f=\phi_{f'}+\nabla\phi_{f'}\cdot(\vec r_f-\vec r_{f'})$
7. 计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$
8. 返回步骤5再迭代一次。

---

例1

上图所示网格，单元形心 $C$与其邻近单元形心 $F_1$到 $F_6$的坐标为
$C(13, 11), \quad F_1(4.5,9.5), \quad F_2(8,3), \quad F_3(17,3.5), \quad F_4(22,10), \quad F_5(16,20), \quad F_6(7,18)$
角点 $n_1$到 $n_2$的坐标为
$n_1(9,14), \quad n_2(8,8), \quad n_3(12,5), \quad n_4(17,9), \quad n_5(17.5,14), \quad n_6(12,17)$
变量 $\phi$在单元形心的值为
$\phi_C=167, \quad \phi_{F_1}=56.75, \quad \phi_{F_2}=35, \quad \phi_{F_3}=80, \quad \phi_{F_4}=252, \quad \phi_{F_5}=356, \quad\phi_{F_6}=151$
C单元的邻近单元形心处变量 $\phi$的梯度值 $\nabla\phi$为
$\nabla\phi_{F_1}=10.5\bold i+5.5 \bold j ，\quad \nabla\phi_{F_2}=4\bold i+9 \bold j, \quad \nabla\phi_{F_3}=4.5\bold i+18 \bold j, \\ \nabla\phi_{F_4}=11 \bold i+23 \bold j, \quad \nabla\phi_{F_5}=21 \bold i+17 \bold j, \quad \nabla\phi_{F_6}=19\bold i+8 \bold j$
单元C的体积
$V_C=76$
求解单元C形心处的梯度值 $\nabla\phi_C$，使用以下方法：
a. Green-Gauss方法，不含修正
b. Green-Gauss方法，含skewness correction， $f'$选择为线段CF的中点

解法a. Green-Gauss方法求解 $\nabla\phi_C$，不含修正
计算面形心（2维状况下就是面中心）坐标值
$x_{f_1}=(x_{n_1}+x_{n_2})/2=(9+8)/2=8.5 \\ y_{f_1}=(y_{n_1}+y_{n_2})/2=(14+8)/2=11$
即
$f_1(8.5, 11)$
一样方法算得单元C的其他面的形心坐标
$f_2(10, 6.5), \quad f_3(14.5, 7), \quad f_4(17.25, 11.5), \quad f_5(14.75, 15.5), \quad f_6(10.5,15.5)$
计算面积矢量 $\vec S_{f_1}$
$\vec S_{f_1}=\Delta y\bold i - \Delta x\bold j=(y_{n_2}-y_{n_1})\bold i - (x_{n_2}-x_{n_1}) \bold j \\ =(8-14)\bold i - (8-9) \bold j=-6\bold i + \bold j$
一样方法算得其他面的面积矢量
$\vec S_{f_2}=-3\bold i -4 \bold j \quad \vec S_{f_3}=4\bold i -5 \bold j \quad \vec S_{f_4}=5\bold i -0.5 \bold j \quad \vec S_{f_5}=3\bold i +5.5 \bold j \quad \vec S_{f_6}=-3\bold i +3 \bold j$
计算插值系数 $g_{C_1}$
$g_{C_1}=\frac{F_1 f_1}{F_1 f_1+Cf_1} \\ \quad \\ F_1 f_1=\sqrt{(4.5-8.5)^2+(9.5-11)^2}=4.272 \\ \quad \\ Cf_1=\sqrt{(13-8.5)^2+(11-11)^2}=4.5$
算得
$g_{C_1}=0.487$
一样，算得其它面的插值系数
$g_{C_2}=0.427, \quad g_{C_3}=0.502, \quad g_{C_4}=0.538, \quad g_{C_5}=0.492, \quad g_{C_6}=0.455$
计算面形心的变量值 $\phi_f$
$\phi_{f_1}=g_{C_1}\phi_C+(1-g_{C_1})\phi_{F_1}=0.487*167+(1-0.487)*56.75=110.442$
一样算得其它面形心的变量值
$\phi_{f_2}=91.364, \quad \phi_{f_3}=123.674, \quad \phi_{f_4}=206.27, \quad \phi_{f_5}=263.012, \quad \phi_{f_6}=158.28$
计算单元形心C处的梯度值
$\nabla\phi_C = \frac{1}{V_C}\sum_{f=1}^6\phi_{f} \vec S_f \\ =\frac{1}{76}\left\{ \begin{bmatrix} 110.442(-6\bold i + \bold j)+91.364(3\bold i -4 \bold j)+123.674(4\bold i -5 \bold j) \\ +206.27(5\bold i -0.5 \bold j)+263.012(3\bold i +5.5 \bold j)+158.28(-3\bold i +3 \bold j) \end{bmatrix} \right\} \\ = 11.889 \bold i +12.433 \bold j$

解法b. Green-Gauss方法，含skewness correction， $f'$选择为线段CF的中点
计算CF的中点 $f'$处的变量值，使用公式 $\displaystyle \phi_{f'}=(\phi_C+\phi_F)/2$，得
$\phi_{f_1}=111.875, \quad \phi_{f_2}=101 \quad \phi_{f_3}=123.5, \quad \phi_{f_4}=209.5, \quad \phi_{f_5}=261.5, \quad \phi_{f_6}=159$
计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f$，得
$\nabla\phi_C = \frac{1}{V_C}\sum_{f=1}^6\phi_{f} \vec S_f \\ =\frac{1}{76}\left\{ \begin{bmatrix} 111.875(-6\bold i + \bold j)+101(3\bold i -4 \bold j)+123.5(4\bold i -5 \bold j) \\ +209.5(5\bold i -0.5 \bold j)+261.5(3\bold i +5.5 \bold j)+159(-3\bold i +3 \bold j) \end{bmatrix} \right\} \\ = 11.510 \bold i + 11.854 \bold j$
接下来修正梯度值
计算 $\displaystyle \vec d_f = \vec r_f-\frac{\vec r_C+\vec r_F}{2}$，得
$\vec d_{f_1}=-0.25\bold i + 0.75\bold j, \quad \vec d_{f_2}=-0.5\bold i -0.5\bold j, \quad \vec d_{f_3}=-0.5\bold i -0.25\bold j, \\ \vec d_{f_4}=-0.25\bold i + \bold j, \quad \vec d_{f_5}=-0.25\bold i, \quad \vec d_{f_6}=0.5\bold i + \bold j$
更新 $\phi_{f}$使用 $\displaystyle \phi_f=\phi_{f'}+\frac{(\nabla\phi)_C+(\nabla\phi)_F}{2}\cdot\left(\vec r_f-\frac{\vec r_C+\vec r_F}{2}\right)$，得
$\phi_{f_1}=115.631, \quad \phi_{f_2}=91.909 \quad \phi_{f_3}=115.766, \quad \phi_{f_4}=224.113 \quad \phi_{f_5}=265.564, \quad \phi_{f_6}=176.554$
计算 $\nabla\phi_C$使用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$，得
$\nabla\phi_C = 11.594 \bold i + 13.781 \bold j$
再迭代一次修正过程，得
$\nabla\phi_C = 11.652 \bold i + 15.793 \bold j$
若是继续迭代修正，你会发现，这个值压根不会收敛，并且会愈来愈大直至崩溃，也就是说，修正上一次两次就OK了，既不会浪费时间，也不会让值太离谱了。

---

方法2：扩展框架

因为面是由角点构成的，因此用角点处的值的平均来算得面形心的值就理所固然了，那么角点处的值要如何获取呢？用围绕该角点的单元形心处的值来加权平均计算便可。比较拗口哈，看下图：

用 $F_1, F_2, F_3$处的值来算得 $n_1$处的值，用 $F_2, F_3, F_4$处的值来算得 $n_2$处的值，最后再用 $n_1,n_2$处的值来算得 $f$处的值，便可。

角点处的值由围绕角点的单元形心处的值加权平均算出，加权系数取为距离的倒数（距离越远，影响越小），即
$\displaystyle \phi_n=\frac{\displaystyle \sum_{k=1}^{NB(n)}\frac{\phi_{F_k}}{||\vec r_n-\vec r_{F_k}||}}{\displaystyle \sum_{k=1}^{NB(n)} \frac{1}{{||\vec r_n-\vec r_{F_k}||}}}$
其中 $n$表明角点， $F_k$表明邻近单元， $NB(n)$为围绕角点 $n$的单元总数， $||\vec r_n-\vec r_{F_k}||$为角点到邻近单元形心的距离。

角点值获得后，面形心的值也可获得，以2维问题为例， $\phi_f$为
$\phi_f=\frac{\phi_{n1}+\phi_{n2}}{2}$
紧接着，可算得单元形心的梯度 $\nabla\phi_C$为
$\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$
对于3维状况，面形心的值 $\phi_f$由角点值的距离加权平均计算，即
$\displaystyle \phi_f=\frac{\displaystyle \sum_{k=1}^{nb(f)}\frac{\phi_{n_k}}{||\vec r_f-\vec r_{n_k}||}}{\displaystyle \sum_{k=1}^{nb(f)} \frac{1}{{||\vec r_f-\vec r_{n_k}||}}}$
而单元形心梯度值的计算方法照旧不变。

---

例2 与例1数据相同，用方法2：扩展框架计算单元中心梯度值

用方法2：扩展框架计算单元中心梯度值
先计算角点与其周围单元形心的距离值 $d_{nF_k}=||\vec r_n-\vec r_{F_k}||$，得
$d_{n_1C}=5, \quad d_{n_1F_6}=4.472, \quad d_{n_1F_1}=6.364$
$d_{n_2C}=5.831, \quad d_{n_2F_1}=3.808, \quad d_{n_2F_2}=5$
$d_{n_3C}=6.083, \quad d_{n_3F_2}=4.472, \quad d_{n_3F_3}=5.220$
$d_{n_4C}=4.472, \quad d_{n_4F_3}=5.5, \quad d_{n_4F_4}=5.099$
$d_{n_5C}=5.408, \quad d_{n_5F_4}=6.021, \quad d_{n_5F_5}=6.185$
$d_{n_6C}=6.083, \quad d_{n_6F_5}=5, \quad d_{n_6F_6}=5.099$
用距离倒数做为权系数，由角点周围单元形心值获取角点值，即
$\displaystyle \phi_n=\frac{\displaystyle \sum_{k=1}^{NB(n)}\frac{\phi_{F_k}}{||\vec r_n-\vec r_{F_k}||}}{\displaystyle \sum_{k=1}^{NB(n)} \frac{1}{{||\vec r_n-\vec r_{F_k}||}}}$
得
$\phi_{n_1}=131.008, \quad \phi_{n_2}=79.708 \quad \phi_{n_3}=87.317, \quad \phi_{n_4}=168.416 \quad \phi_{n_5}=254.144, \quad \phi_{n_6}=228.840$
计算面形心处的变量值，用面角点值平均来获取，例如 $\phi_{f1}=(\phi_{n1} + \phi_{n2})/2$，得
$\phi_{f_1}=105.358, \quad \phi_{f_2}=83.512 \quad \phi_{f_3}=127.866, \quad \phi_{f_4}=211.280 \quad \phi_{f_5}=241.492, \quad \phi_{f_6}=179.924$
计算单元中心梯度值，用 $\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$，得
$\nabla\phi_C = 11.446 \bold i + 11.767 \bold j$

---

3 Least-Square Gradient（最小二乘梯度）

最小二乘法计算梯度，提供了更高的精度，以及更加灵活的选择，用的框架点也更多，然而其须要计算较多的加权系数，固然计算消耗也比较大。

考虑上图，单元 $C$有第1层邻近单元和第2层邻近单元，那么，若是单元形心的梯度 $\nabla\phi_C$是精确的话，有
$\phi_F=\phi_C+(\nabla\phi)_C\cdot(\vec r_F - \vec r_C)$
在最小二乘法中，设法让上式算得的单元邻近单元值的加权加和值最小，即找到以下函数的最小值
$G_C=\sum_{k=1}^{NB(C)}\{ w_k[\phi_{F_k}-(\phi_C+(\nabla\phi)_C\cdot\vec r_{CF_k})]^2 \} \\ =\sum_{k=1}^{NB(C)}\left\{ w_k\left[\Delta\phi_k - \left( \Delta x_k\left(\frac{\partial\phi}{\p$