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$