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

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

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

$${\left(\frac{\partial \varphi }{\partial x}\right)}_e=\frac{{\varphi }_E-{\varphi }_C}{x_E-x_C}=\frac{{\varphi }_E-{\varphi }_C}{\delta x_e}$$
单元形心$C$处的变量梯度也可得出
$${\left(\frac{\partial \varphi }{\partial x}\right)}_C=\frac{{\varphi }_e-{\varphi }_w}{x_e-x_w}=\frac{\frac{{\varphi }_E+{\varphi }_C}{2}-\frac{{\varphi }_C+{\varphi }_W}{2}}{\Delta x_C}=\frac{{\varphi }_E-{\varphi }_W}{2\Delta x_C}$$
这个常叫做“中心差分”，对于2维情况，与此类似，可得

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

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

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

方法1：紧致框架

对于上图(a)和(b)所示的2维和3维情况，一种最简单的方法是直接由面两侧单元形心值平均来获得面上变量的值，即
$${\varphi }_f=g_C{\varphi }_C+(1-g_C){\varphi }_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$为两点之间的距离值，若该面恰好位于两单元中心的中间位置，则有
$${\varphi }_f=\frac{{\varphi }_C+{\varphi }_F}{2}$$
这个方法非常简便易行，然而遗憾的是，只有当线段$CF$和面${\vec{S}}_f$的交点与面${\vec{S}}_f$的中心重合时（即上图（a）（b）情况），该方法才能保证2阶精度，而大多数情况下，直接由上式算得的梯度精度是没法保证的。

然而，对于通常的非正交网格和非结构网格来说，是没有办法来保证这个条件的，比如上图中的（c）（d）情况，网格的偏斜（非正交，skewness（non-conjunctionality））使得线段$CF$和面${\vec{S}}_f$交点为$f^{\prime }$，与面形心$f$是不重合的。此时，需要把插值得到的${\varphi }_{f^{\prime }}$修正以得到${\varphi }_f$，修正方程为
$${\varphi }_f={\varphi }_{f^{\prime }}+correction={\varphi }_{f^{\prime }}+(\nabla \varphi {)}_{f^{\prime }}\cdot ({\vec{r}}_f-{\vec{r}}_{f^{\prime }})$$
即，用梯度$(\nabla \varphi {)}_{f^{\prime }}$来修正，修正也可以展开为如下形式：
ϕ f = g C { ϕ C + ( ∇ ϕ ) C ⋅ ( r ⃗ f − r ⃗ C ) } + ( 1 − g C ) { ϕ F + ( ∇ ϕ ) F ⋅ ( r ⃗ f − r ⃗ F ) } = ϕ f ′ + g C ( ∇ ϕ ) C ⋅ ( r ⃗ f − r ⃗ C ) + ( 1 − g C ) ( ∇ ϕ ) F ⋅ ( r ⃗ f − r ⃗ 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) ϕf​=gC​{ϕC​+(∇ϕ)C​⋅(r f​−r C​)}+(1−gC​){ϕF​+(∇ϕ)F​⋅(r f​−r F​)}=ϕf′​+gC​(∇ϕ)C​⋅(r f​−r C​)+(1−gC​)(∇ϕ)F​⋅(r f​−r F​)
即， { ϕ C + ( ∇ ϕ ) C ⋅ ( r ⃗ f − r ⃗ C ) } \{\phi_C+(\nabla\phi)_C\cdot(\vec r_f-\vec r_C)\} {ϕC​+(∇ϕ)C​⋅(r f​−r C​)}是把$C$处的值修正到$f$处，而 { ϕ F + ( ∇ ϕ ) F ⋅ ( r ⃗ f − r ⃗ F ) } \{ \phi_F+(\nabla\phi)_F\cdot(\vec r_f-\vec r_F) \} {ϕF​+(∇ϕ)F​⋅(r f​−r F​)}是把$F$处的值修正到$f$处，然后再做加权插值处理，就得到了${\varphi }_f$。

不难发现，${\varphi }_f$的修正计算要用到$(\nabla \varphi {)}_C$和$(\nabla \varphi {)}_F$，而$(\nabla \varphi {)}_C$和$(\nabla \varphi {)}_F$则是用${\varphi }_f$算得的，也就是说，这里${\varphi }_f$的计算是不能一蹴而就的，而是一个迭代计算的过程，但是过多的迭代反而会造成解的振荡，一般做两次迭代就好了。

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

选择1

点$f^{\prime }$就选择在线段$CF$和面${\vec{S}}_f$的交点处，以$\vec{n}$代表面积单位矢量，即，$\vec{n}={\vec{S}}_f/||{\vec{S}}_f||$，以$\vec{e}$代表沿着$CF$的单位矢量，即$\vec{e}=\overrightarrow{CF}/||\overrightarrow{CF}||$，则$f^{\prime }$的位置可由几何关系算出，为
$${\vec{r}}_{f^{\prime }}=\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^{\prime }$的向量$\overrightarrow{C{f}^{\prime }}$，再与${\vec{r}}_C$相加便是点$f^{\prime }$的位置矢量。
得到$f^{\prime }$的位置后，可以直接算出$g_C$的值
$$g_C=\frac{||{\vec{r}}_F-{\vec{r}}_{f^{\prime }}||}{||{\vec{r}}_F-{\vec{r}}_C||}=\frac{d_{F{f}^{\prime }}}{d_{FC}}$$
计算流程如下：

1. 计算${\varphi }_{f^{\prime }}$使用${\varphi }_{f^{\prime }}=g_C{\varphi }_C+(1-g_C){\varphi }_F$
2. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_{f^{\prime }}{\vec{S}}_f$
   接下来用下面步骤来修正梯度场
3. 更新${\varphi }_f$使用${\varphi }_f={\varphi }_{f^{\prime }}+g_C(\nabla \varphi {)}_C\cdot ({\vec{r}}_f-{\vec{r}}_C)+(1-g_C)(\nabla \varphi {)}_F\cdot ({\vec{r}}_f-{\vec{r}}_F)$
4. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$
5. 返回步骤3再迭代一次。

选择2

点$f^{\prime }$就选择在线段$CF$的中点，相当于$g_C=0.5$，这就使得计算简单多了，其计算流程如下：

1. 计算${\varphi }_{f^{\prime }}$使用${\varphi }_{f^{\prime }}=\frac{{\varphi }_C+{\varphi }_F}{2}$
2. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_{f^{\prime }}{\vec{S}}_f$
   接下来用下面步骤来修正梯度场
3. 更新${\varphi }_f$使用${\varphi }_f={\varphi }_{f^{\prime }}+\frac{(\nabla \varphi {)}_C+(\nabla \varphi {)}_F}{2}\cdot \left({\vec{r}}_f-\frac{{\vec{r}}_C+{\vec{r}}_F}{2}\right)$
4. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$
5. 返回步骤3再迭代一次。

选择3

点$f^{\prime }$的位置选择要保证距离$f{f}^{\prime }$是最短的，即，$f{f}^{\prime }$是垂直于$CF$的，这使得第1步迭代计算变得更加准确。$f^{\prime }$的位置计算很简单，直接把向量$\overrightarrow{Cf}$投影到$\overrightarrow{CF}$上就妥了，即
$${\vec{r}}_{f^{\prime }}={\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^{\prime }}$使用${\vec{r}}_{f^{\prime }}={\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$使用$g_C=\frac{||{\vec{r}}_F-{\vec{r}}_{f^{\prime }}||}{||{\vec{r}}_F-{\vec{r}}_C||}$
3. 计算${\varphi }_{f^{\prime }}$使用${\varphi }_{f^{\prime }}=g_C{\varphi }_C+(1-g_C){\varphi }_F$
4. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_{f^{\prime }}{\vec{S}}_f$
   接下来用下面步骤来修正梯度场
5. 计算$\nabla {\varphi }_{f^{\prime }}$使用$\nabla {\varphi }_{f^{\prime }}=g_C\nabla {\varphi }_C+(1-g_C)\nabla {\varphi }_F$
6. 更新${\varphi }_f$使用${\varphi }_f={\varphi }_{f^{\prime }}+\nabla {\varphi }_{f^{\prime }}\cdot ({\vec{r}}_f-{\vec{r}}_{f^{\prime }})$
7. 计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$
8. 返回步骤5再迭代一次。

---

例1

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

解法a. Green-Gauss方法求解$\nabla {\varphi }_C$，不含修正
计算面形心（2维情况下就是面中心）坐标值
$$\begin{gathered} 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 \end{gathered}$$
即
$$f_1(8.5,11)$$
同样方法算得单元C的其余面的形心坐标
$$f_2(10,6.5),f_3(14.5,7),f_4(17.25,11.5),f_5(14.75,15.5),f_6(10.5,15.5)$$
计算面积矢量${\vec{S}}_{f_1}$
$$\begin{gathered} {\vec{S}}_{f_1}=\Delta y\mathbf{i}-\Delta x\mathbf{j}=(y_{n_2}-y_{n_1})\mathbf{i}-(x_{n_2}-x_{n_1})\mathbf{j} \\ =(8-14)\mathbf{i}-(8-9)\mathbf{j}=-6\mathbf{i}+\mathbf{j} \end{gathered}$$
同样方法算得其余面的面积矢量
$${\vec{S}}_{f_2}=-3\mathbf{i}-4\mathbf{j}{\vec{S}}_{f_3}=4\mathbf{i}-5\mathbf{j}{\vec{S}}_{f_4}=5\mathbf{i}-0.5\mathbf{j}{\vec{S}}_{f_5}=3\mathbf{i}+5.5\mathbf{j}{\vec{S}}_{f_6}=-3\mathbf{i}+3\mathbf{j}$$
计算插值系数$g_{C_1}$
$$\begin{gathered} g_{C_1}=\frac{F_1{f}_1}{F_1{f}_1+C{f}_1} \\ {F}_1{f}_1=\sqrt{(4.5-8.5{)}^2+(9.5-11{)}^2}=4.272 \\ C{f}_1=\sqrt{(13-8.5{)}^2+(11-11{)}^2}=4.5 \end{gathered}$$
算得
$$g_{C_1}=0.487$$
同样，算得其它面的插值系数
$$g_{C_2}=0.427,g_{C_3}=0.502,g_{C_4}=0.538,g_{C_5}=0.492,g_{C_6}=0.455$$
计算面形心的变量值${\varphi }_f$
$${\varphi }_{f_1}=g_{C_1}{\varphi }_C+(1-g_{C_1}){\varphi }_{F_1}=0.487\ast 167+(1-0.487)\ast 56.75=110.442$$
同样算得其它面形心的变量值
$${\varphi }_{f_2}=91.364,{\varphi }_{f_3}=123.674,{\varphi }_{f_4}=206.27,{\varphi }_{f_5}=263.012,{\varphi }_{f_6}=158.28$$
计算单元形心C处的梯度值
$$\begin{gathered} \nabla {\varphi }_C=\frac{1}{V_C}\sum _{f=1}^6{\varphi }_f{\vec{S}}_f \\ =\frac{1}{76}\left\{\left[\begin{array}{c}110.442(-6\mathbf{i}+\mathbf{j})+91.364(3\mathbf{i}-4\mathbf{j})+123.674(4\mathbf{i}-5\mathbf{j})\\ +206.27(5\mathbf{i}-0.5\mathbf{j})+263.012(3\mathbf{i}+5.5\mathbf{j})+158.28(-3\mathbf{i}+3\mathbf{j})\end{array}\right]\right\} \\ =11.889\mathbf{i}+12.433\mathbf{j} \end{gathered}$$

解法b. Green-Gauss方法，含skewness correction，$f^{\prime }$选择为线段CF的中点
计算CF的中点$f^{\prime }$处的变量值，使用公式${\varphi }_{f^{\prime }}=({\varphi }_C+{\varphi }_F)/2$，得
$${\varphi }_{f_1}=111.875,{\varphi }_{f_2}=101{\varphi }_{f_3}=123.5,{\varphi }_{f_4}=209.5,{\varphi }_{f_5}=261.5,{\varphi }_{f_6}=159$$
计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_{f^{\prime }}{\vec{S}}_f$，得
$$\begin{gathered} \nabla {\varphi }_C=\frac{1}{V_C}\sum _{f=1}^6{\varphi }_f{\vec{S}}_f \\ =\frac{1}{76}\left\{\left[\begin{array}{c}111.875(-6\mathbf{i}+\mathbf{j})+101(3\mathbf{i}-4\mathbf{j})+123.5(4\mathbf{i}-5\mathbf{j})\\ +209.5(5\mathbf{i}-0.5\mathbf{j})+261.5(3\mathbf{i}+5.5\mathbf{j})+159(-3\mathbf{i}+3\mathbf{j})\end{array}\right]\right\} \\ =11.510\mathbf{i}+11.854\mathbf{j} \end{gathered}$$
接下来修正梯度值
计算${\vec{d}}_f={\vec{r}}_f-\frac{{\vec{r}}_C+{\vec{r}}_F}{2}$，得
$$\begin{gathered} {\vec{d}}_{f_1}=-0.25\mathbf{i}+0.75\mathbf{j},{\vec{d}}_{f_2}=-0.5\mathbf{i}-0.5\mathbf{j},{\vec{d}}_{f_3}=-0.5\mathbf{i}-0.25\mathbf{j}, \\ {\vec{d}}_{f_4}=-0.25\mathbf{i}+\mathbf{j},{\vec{d}}_{f_5}=-0.25\mathbf{i},{\vec{d}}_{f_6}=0.5\mathbf{i}+\mathbf{j} \end{gathered}$$
更新${\varphi }_f$使用${\varphi }_f={\varphi }_{f^{\prime }}+\frac{(\nabla \varphi {)}_C+(\nabla \varphi {)}_F}{2}\cdot \left({\vec{r}}_f-\frac{{\vec{r}}_C+{\vec{r}}_F}{2}\right)$，得
$${\varphi }_{f_1}=115.631,{\varphi }_{f_2}=91.909{\varphi }_{f_3}=115.766,{\varphi }_{f_4}=224.113{\varphi }_{f_5}=265.564,{\varphi }_{f_6}=176.554$$
计算$\nabla {\varphi }_C$使用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$，得
$$\nabla {\varphi }_C=11.594\mathbf{i}+13.781\mathbf{j}$$
再迭代一次修正过程，得
$$\nabla {\varphi }_C=11.652\mathbf{i}+15.793\mathbf{j}$$
如果继续迭代修正，你会发现，这个值压根不会收敛，而且会越来越大直至崩溃，也就是说，修正上一次两次就OK了，既不会浪费时间，也不会让值太离谱了。

---

方法2：扩展框架

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

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

角点处的值由围绕角点的单元形心处的值加权平均算出，加权系数取为距离的倒数（距离越远，影响越小），即
$${\varphi }_n=\frac{\sum _{k=1}^{NB(n)}\frac{{\varphi }_{F_k}}{||{\vec{r}}_n-{\vec{r}}_{F_k}||}}{\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维问题为例，${\varphi }_f$为
$${\varphi }_f=\frac{{\varphi }_{n1}+{\varphi }_{n2}}{2}$$
紧接着，可算得单元形心的梯度$\nabla {\varphi }_C$为
$$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$$
对于3维情况，面形心的值${\varphi }_f$由角点值的距离加权平均计算，即
$${\varphi }_f=\frac{\sum _{k=1}^{nb(f)}\frac{{\varphi }_{n_k}}{||{\vec{r}}_f-{\vec{r}}_{n_k}||}}{\sum _{k=1}^{nb(f)}\frac{1}{||{\vec{r}}_f-{\vec{r}}_{n_k}||}}$$
而单元形心梯度值的计算方法照旧不变。

---

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

用方法2：扩展框架计算单元中心梯度值
先计算角点与其周围单元形心的距离值$d_{n{F}_k}=||{\vec{r}}_n-{\vec{r}}_{F_k}||$，得
$d_{n_{1}C}=5,d_{n_1{F}_6}=4.472,d_{n_1{F}_1}=6.364$
$d_{n_{2}C}=5.831,d_{n_2{F}_1}=3.808,d_{n_2{F}_2}=5$
$d_{n_{3}C}=6.083,d_{n_3{F}_2}=4.472,d_{n_3{F}_3}=5.220$
$d_{n_{4}C}=4.472,d_{n_4{F}_3}=5.5,d_{n_4{F}_4}=5.099$
$d_{n_{5}C}=5.408,d_{n_5{F}_4}=6.021,d_{n_5{F}_5}=6.185$
$d_{n_{6}C}=6.083,d_{n_6{F}_5}=5,d_{n_6{F}_6}=5.099$
用距离倒数作为权系数，由角点周围单元形心值获取角点值，即
$${\varphi }_n=\frac{\sum _{k=1}^{NB(n)}\frac{{\varphi }_{F_k}}{||{\vec{r}}_n-{\vec{r}}_{F_k}||}}{\sum _{k=1}^{NB(n)}\frac{1}{||{\vec{r}}_n-{\vec{r}}_{F_k}||}}$$
得
$${\varphi }_{n_1}=131.008,{\varphi }_{n_2}=79.708{\varphi }_{n_3}=87.317,{\varphi }_{n_4}=168.416{\varphi }_{n_5}=254.144,{\varphi }_{n_6}=228.840$$
计算面形心处的变量值，用面角点值平均来获取，例如${\varphi }_{f1}=({\varphi }_{n1}+{\varphi }_{n2})/2$，得
$${\varphi }_{f_1}=105.358,{\varphi }_{f_2}=83.512{\varphi }_{f_3}=127.866,{\varphi }_{f_4}=211.280{\varphi }_{f_5}=241.492,{\varphi }_{f_6}=179.924$$
计算单元中心梯度值，用$\nabla {\varphi }_C=\frac{1}{V_C}\sum _{f-nb(C)}{\varphi }_f{\vec{S}}_f$，得
$$\nabla {\varphi }_C=11.446\mathbf{i}+11.767\mathbf{j}$$

---

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

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

考虑上图，单元$C$有第1层邻近单元和第2层邻近单元，那么，如果单元形心的梯度$\nabla {\varphi }_C$是精确的话，有
$${\varphi }_F={\varphi }_C+(\nabla \varphi {)}_C\cdot ({\vec{r}}_F-{\vec{r}}_C)$$
在最小二乘法中，设法让上式算得的单元邻近单元值的加权加和值最小，即找到如下函数的最小值
G C = ∑ k = 1 N B ( C ) { w k [ ϕ F k − ( ϕ C + ( ∇ ϕ ) C ⋅ r ⃗ C F k ) ] 2 } = ∑ k = 1 N B ( C ) { w k [ Δ ϕ k − ( Δ x k ( ∂ ϕ ∂ x ) C + Δ y k ( ∂ ϕ ∂ y ) C + Δ z k ( ∂ ϕ ∂ z ) C ) ] 2 } 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}{\partial x}\right)_C + \Delta y_k\left(\frac{\partial\phi}{\partial y}\right)_C + \Delta z_k\left(\frac{\partial\phi}{\partial z}\right)_C \right) \right]^2 \right\} GC​=k=1∑NB(C)​{wk​[ϕFk​​−(ϕC​+(∇ϕ)C​⋅r CFk​​)]2}=k=1∑NB(C)​{wk​[Δϕk​−(Δxk​(∂x∂ϕ​)C​+Δyk​(∂y∂ϕ​)C​+Δzk​(∂z∂ϕ​)C​)]2}
其中$w_k$为加权系数
$$\begin{gathered} \Delta {\varphi }_k={\varphi }_{F_k}-{\varphi }_C \\ \Delta x_k={\vec{r}}_{C{F}_k}\cdot \vec{i} \\ \Delta y_k={\vec{r}}_{C{F}_k}\cdot \vec{j} \\ \Delta z_k={\vec{r}}_{C{F}_k}\cdot \vec{k} \end{gathered}$$
$G_C$函数的最小值应满足如下条件
$$\frac{\partial G_C}{\partial \left(\frac{\partial \varphi }{\partial x}\right)}=\frac{\partial G_C}{\partial \left(\frac{\partial \varphi }{\partial y}\right)}=\frac{\partial G_C}{\partial \left(\frac{\partial \varphi }{\partial z}\right)}=0$$
得到了如下含三个未知量的三个方程
$$\begin{gathered} \sum _{k=1}^{NB(C)}\left\{2w_k\Delta x_k\left[-\Delta {\varphi }_k+\Delta x_k{\left(\frac{\partial \varphi }{\partial x}\right)}_C+\Delta y_k{\left(\frac{\partial \varphi }{\partial y}\right)}_C+\Delta z_k{\left(\frac{\partial \varphi }{\partial z}\right)}_C\right]\right\}=0 \\ \sum _{k=1}^{NB(C)}\left\{2w_k\Delta y_k\left[-\Delta {\varphi }_k+\Delta x_k{\left(\frac{\partial \varphi }{\partial x}\right)}_C+\Delta y_k{\left(\frac{\partial \varphi }{\partial y}\right)}_C+\Delta z_k{\left(\frac{\partial \varphi }{\partial z}\right)}_C\right]\right\}=0 \\ \sum _{k=1}^{NB(C)}\left\{2w_k\Delta z_k\left[-\Delta {\varphi }_k+\Delta x_k{\left(\frac{\partial \varphi }{\partial x}\right)}_C+\Delta y_k{\left(\frac{\partial \varphi }{\partial y}\right)}_C+\Delta z_k{\left(\frac{\partial \varphi }{\partial z}\right)}_C\right]\right\}=0 \end{gathered}$$
可以写成矩阵形式
$$\left[\begin{array}{ccc}\sum _{k=1}^{NB(C)}{w}_k\Delta x_k\Delta x_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta x_k\Delta y_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta x_k\Delta z_k\\ \sum _{k=1}^{NB(C)}{w}_k\Delta y_k\Delta x_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta y_k\Delta y_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta y_k\Delta z_k\\ \sum _{k=1}^{NB(C)}{w}_k\Delta z_k\Delta x_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta z_k\Delta y_k& {\sum }_{k=1}^{NB(C)}{w}_k\Delta z_k\Delta z_k\end{array}\right]\left[\begin{array}{c}{\left(\frac{\partial \varphi }{\partial x}\right)}_C\\ \\ {\left(\frac{\partial \varphi }{\partial y}\right)}_C\\ \\ {\left(\frac{\partial \varphi }{\partial z}\right)}_C\end{array}\right]=\left[\begin{array}{c}\sum _{k=1}^{NB(C)}{w}_k\Delta x_k\Delta {\varphi }_k\\ \sum _{k=1}^{NB(C)}{w}_k\Delta y_k\Delta {\varphi }_k\\ \sum _{k=1}^{NB(C)}{w}_k\Delta z_k\Delta {\varphi }_k\end{array}\right]$$
求解上述方程组，便可得到单元形心的梯度值$\nabla {\varphi }_C$，只要左端项矩阵不奇异，即存在解。$w_k$的选择是至关重要的，比如把它设成1就不太合适了，常见的设置方法是距离的倒数，即
$$w_k=\frac{1}{||{\vec{r}}_{F_k}-{\vec{r}}_C||}=\frac{1}{\sqrt{\Delta x_{F_k}^2+\Delta y_{F_k}^2+\Delta z_{F_k}^2}}$$
也可以用距离的$n$次幂的倒数来做加权系数，即
$$w_k=\frac{1}{||{\vec{r}}_{F_k}-{\vec{r}}_C|{|}^n}$$
幂次$n$可以是1,2,3等。

4 梯度插值到面

同样地，由面两侧的单元值插值获取，然后再做修正，修正的原则是让面梯度沿着$CF$方向的值与由$C$和$F$处$\varphi$值定义的局部梯度相等，如此，面梯度$\nabla {\varphi }_f$为
$$\nabla {\varphi }_f=\overline{\nabla {\varphi }_f}+\left[\frac{{\varphi }_F-{\varphi }_C}{d_{CF}}-(\overline{\nabla {\varphi }_f}\cdot {\vec{e}}_{CF})\right]{\vec{e}}_{CF}$$
其中
$$\overline{\nabla {\varphi }_f}=g_C\nabla {\varphi }_C+g_F\nabla {\varphi }_F$$
$${\vec{e}}_{CF}=\frac{{\vec{d}}_{CF}}{d_{CF}}$$
$${\vec{d}}_{CF}={\vec{r}}_F-{\vec{r}}_C$$

补充一些个人领悟：这里从单元梯度计算面梯度的过程，其思想是这样子的，在${\vec{e}}_{CF}$方向上的梯度由$\frac{{\varphi }_F-{\varphi }_C}{d_{CF}}{\vec{e}}_{CF}$给出；而在垂直于${\vec{e}}_{CF}$方向上的梯度则由原本两侧owner和neighbour单元（所属和邻近单元）的梯度值插值到面上的梯度值$\overline{\nabla {\varphi }_f}=g_C\nabla {\varphi }_C+g_F\nabla {\varphi }_F$在垂直于${\vec{e}}_{CF}$方向上的分量来给出，即$\overline{\nabla {\varphi }_f}-(\overline{\nabla {\varphi }_f}\cdot {\vec{e}}_{CF}){\vec{e}}_{CF}$；如此一来，把${\vec{e}}_{CF}$方向和垂直于${\vec{e}}_{CF}$方向这两个方向上梯度值复合在一起，便是面上的梯度值了，即
$$\nabla {\varphi }_f=\overline{\nabla {\varphi }_f}+\left[\frac{{\varphi }_F-{\varphi }_C}{d_{CF}}-(\overline{\nabla {\varphi }_f}\cdot {\vec{e}}_{CF})\right]{\vec{e}}_{CF}$$

5 代码讲解

5.1 uFVM

在uFVM中，cfdComputeGradientGauss0和cfdComputeGradientNodal函数是用来计算单元梯度的，Gauss梯度函数输入的是单元变量值phi的数组，返回的是单元梯度值phiGrad的数组，其中面形心变量是用几何权系数从两侧单元形心值插值获取的，并没有做非正交偏斜修正操作，因为这个函数里的0表明是不带修正的。cfdComputeGradientGauss0和cfdComputeGradientNodal函数分别对应本章理论部分的紧致框架和扩展框架，cfdComputeGradientGauss0函数及个人理解如下

function phiGrad = cfdComputeGradientGauss0(phi, theMesh)
%===================================================
% written by the CFD Group @ AUB, Fall 2006
%===================================================
%
if(nargin<2)    % 如果函数输入参数小于2个，则提取网格信息theMesh
    theMesh = cfdGetMesh;
end

theSize = size(phi);    % 输入的单元场phi的维数 [长度, 分量数目]
theNumberOfComponents = theSize(2); % 变量phi的分量数目（标量1、矢量3、张量9）
% 如果分量数目多于3个，无法处理，直接报错，
% 因为标量的梯度为矢量，而矢量的梯度为张量，那么张量的梯度是啥呢？？没有呢~
if(theNumberOfComponents > 3)
    echo('********* ERROR **********');
    exit;
end
%-----------------------------------------------------
% INTERIOR FACES contribution to gradient
%-----------------------------------------------------
iFaces = 1:theMesh.numberOfInteriorFaces;   % 内部面标识列表
iBFaces = theMesh.numberOfInteriorFaces+1:theMesh.numberOfFaces; % 边界面标识列表
iElements = 1:theMesh.numberOfElements;     % 单元标识列表
iBElements = theMesh.numberOfElements+1:theMesh.numberOfElements+...
    theMesh.numberOfBFaces; % 边界单元标识列表

iOwners = [theMesh.faces(iFaces).iOwner]';  % 内部面owner单元标识列表
iNeighbours = [theMesh.faces(iFaces).iNeighbour]';  % 内部面neighbour单元标识列表

Sf = [theMesh.faces(iFaces).Sf]';           % 内部面面积矢量列表
gf = [theMesh.faces(iFaces).gf]';           % 内部面的插值几何权重系数列表

%-----------------------------------------------------
% Initialize phiGrad Array
%-----------------------------------------------------

% phiGrad存储在单元场上，总共有单元数目+边界单元数目个量，每个量有分量*3个值
phiGrad = zeros(theMesh.numberOfElements+theMesh.numberOfBElements, ...
    3,theNumberOfComponents); 

% 每个分量依次计算梯度值
for iComponent=1:theNumberOfComponents
    % 插值由面两侧单元形心处变量的值计算面形心处变量的值
    phi_f = gf.*phi(iNeighbours,iComponent) + (1-gf).*phi(iOwners,iComponent);
    % 对内部面做循环，将 面形心处变量值*面面积矢量 累加到面的own单元形心梯度值上
    % 累减到面的neighbour单元形心梯度值上去
    for iFace=iFaces
        phiGrad(iOwners(iFace),:,iComponent) = ...
                phiGrad(iOwners(iFace),:,iComponent) + phi_f(iFace)*Sf(iFace,:);
        phiGrad(iNeighbours(iFace),:,iComponent) = ...
            phiGrad(iNeighbours(iFace),:,iComponent) - phi_f(iFace)*Sf(iFace,:);
    end
end

%-----------------------------------------------------
% BOUNDARY FACES contribution to gradient 边界面对梯度的贡献
%-----------------------------------------------------
iBOwners = [theMesh.faces(iBFaces).iOwner]';    % 边界面的own单元标识列表
phi_b = phi(iBElements,:); % 边界单元的变量值
Sb = [theMesh.faces(iBFaces).Sf]';  % 边界面的面积矢量值
for iComponent=1:theNumberOfComponents  % 每个变量得分量依次计算梯度
    % 每片边界patch依次处理    
    % 将 边界面变量值*边界面面积矢量值 累加到边界面的own单元形心梯度值上去
    for k=1:theMesh.numberOfBFaces  
        phiGrad(iBOwners(k),:,iComponent) = ...
            phiGrad(iBOwners(k),:,iComponent) + phi_b(k)*Sb(k,:);
    end
end

%-----------------------------------------------------
% Get Average Gradient by dividing with element volume
%-----------------------------------------------------
volumes = [theMesh.elements(iElements).volume]';    % 单元体积列表
for iComponent=1:theNumberOfComponents % 每个分量依次处理
    for iElement =1:theMesh.numberOfElements % 做体积平均，完成梯度计算
        phiGrad(iElement,:,iComponent) = phiGrad(iElement,:,iComponent) /...
            volumes(iElement);
    end
end
%-----------------------------------------------------
% Set boundary Gradient equal to associated element
% Gradient
%-----------------------------------------------------
% 对于边界单元来说，直接把它的梯度值用边界面的所属单元形心的梯度值赋给即可，
% 因为边界单元的编码是和边界面的编码顺序一致的，所以直接这么操作是没问题的。
phiGrad(iBElements,:,:) = phiGrad(iBOwners,:,:);
end

节点梯度函数，如下所示，与Gauss梯度函数非常相似，只是它用到了cfdInterpolateFromElementsToNodes和cfdInterpolateFromNodesToFaces来把phi值插值到面上去（先把phi插值由单元形心插值到节点，然后再由节点插值到面上），然后再用Gauss定理去计算单元形心梯度，这个函数如下，由于和Guass梯度函数几乎一样，所以不再注释了

%
%=====================================================
% INTERIOR Gradients
%=====================================================
theNumberOfElements = theMesh.numberOfElements;
theNumberOfBElements = theMesh.numberOfBElements;
theNumberOfInteriorFaces = theMesh.numberOfInteriorFaces;
%
phiNodes = cfdInterpolateFromElementsToNodes(phi);
phi_f = cfdInterpolateFromNodesToFaces(phiNodes);
%
phiGrad = zeros(3,theNumberOfElements+theNumberOfBElements);
%-----------------------------------------------------
% INTERIOR FACES contribution to gradient
%-----------------------------------------------------
fvmFaces = theMesh.faces;
% interpolate phi to faces
%
for iFace=1:theNumberOfInteriorFaces
    %
    theFace = fvmFaces(iFace);
    %
    iElement1 = theFace.iOwner;
    iElement2 = theFace.iNeighbour;
    %
    Sf = theFace.Sf;
    %
    %
    phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_f(iFace)*Sf;
    phiGrad(:,iElement2) = phiGrad(:,iElement2) - phi_f(iFace)*Sf;
end
%=====================================================
% BOUNDARY FACES contribution to gradient
%=====================================================
for iBPatch=1:theNumberOfBElements
    %
    iBFace = theNumberOfInteriorFaces+iBPatch;
    iBElement = theNumberOfElements+iBPatch;
    theFace = fvmFaces(iBFace);
    %
    iElement1 = theFace.iOwner;
    %
    Sb = theFace.Sf;
    phi_b = phi(iBElement);
    %
    phiGrad(:,iElement1) = phiGrad(:,iElement1) + phi_b*Sf;
end
%-----------------------------------------------------
% Get Average Gradient by dividing with element volume
%-----------------------------------------------------
for iElement =1:theNumberOfElements
    theElement = fvmElements(iElement);
    phiGrad(:,iElement) = phiGrad(:,iElement)/theElement.volume;
end
%-----------------------------------------------------
% Set boundary Gradient equal to associated element
% Gradient
%-----------------------------------------------------
for iBPatch = 1:theNumberOfBElements
    iBElement = iBPatch+theNumberOfElements;
    iBFace = iBPatch+theNumberOfInteriorFaces;
    theBFace = fvmFaces(iBFace);
    iOwner = theBFace.iOwner;
    phiGrad(:,iBElement) = phiGrad(:,iOwner);
end

面梯度值（面形心处变量梯度值）的计算是由两侧单元梯度值（单元形心处存储的变量梯度值）插值获取的，该操作是通过CFDInterpolateGradientsFromElementsToInteriorFaces函数来实现的。该函数的输入量有插值格式theInterpolationScheme、单元梯度grad，单元phi值，mdot通量值，可以使用迎风或背风格式。如下

function grad_f=cfdInterpolateGradientsFromElementsToInteriorFaces...
    (theInterpolationScheme,grad,phi,mdot_f)
%===================================================
% written by the CFD Group @ AUB, Fall 2006
%===================================================

theMesh= cfdGetMesh;

numberOfInteriorFaces = theMesh.numberOfInteriorFaces;

iOwners = [theMesh.faces(1:numberOfInteriorFaces).iOwner]';
iNeighbours = [theMesh.faces(1:numberOfInteriorFaces).iNeighbour]';
gf = [theMesh.faces(1:numberOfInteriorFaces).gf]';

if(strcmp(theInterpolationScheme,'Average')==1)
    grad_f(:,1) = (1-gf).*grad(iNeighbours,1) + gf.*grad(iOwners,1);
    grad_f(:,2) = (1-gf).*grad(iNeighbours,2) + gf.*grad(iOwners,2);
    grad_f(:,3) = (1-gf).*grad(iNeighbours,3) + gf.*grad(iOwners,3);
    
elseif(strcmp(theInterpolationScheme,'Upwind')==1)
    pos = zeros(size(mdot_f));
    pos((mdot_f>0))=1;
    %
    grad_f(:,1) = pos.*grad(iNeighbours,1) + (1-pos).*grad(iOwners,1);
    grad_f(:,2) = pos.*grad(iNeighbours,2) + (1-pos).*grad(iOwners,2);
    grad_f(:,3) = pos.*grad(iNeighbours,3) + (1-pos).*grad(iOwners,3);
elseif(strcmp(theInterpolationScheme,'Downwind')==1)
    pos = zeros(size(mdot_f));
    pos((mdot_f>0))=1;
    %
    grad_f(:,1) = (1-pos).*grad(iNeighbours,1) + pos.*grad(iOwners,1);
    grad_f(:,2) = (1-pos).*grad(iNeighbours,2) + pos.*grad(iOwners,2);
    grad_f(:,3) = (1-pos).*grad(iNeighbours,3) + pos.*grad(iOwners,3);
elseif(strcmp(theInterpolationScheme,'Average:Corrected')==1)
    grad_f(:,1) = (1-gf).*grad(iNeighbours,1) + gf.*grad(iOwners,1);
    grad_f(:,2) = (1-gf).*grad(iNeighbours,2) + gf.*grad(iOwners,2);
    grad_f(:,3) = (1-gf).*grad(iNeighbours,3) + gf.*grad(iOwners,3);
    
    d_CF = [theMesh.elements(iNeighbours).centroid]' - ...
        [theMesh.elements(iOwners).centroid]';
    dmag=cfdMagnitude(d_CF);
    
    e_CF(:,1)=d_CF(:,1)./dmag;
    e_CF(:,2)=d_CF(:,2)./dmag;
    e_CF(:,3)=d_CF(:,3)./dmag;
    
    local_grad_mag_f = (phi(iNeighbours)-phi(iOwners))./dmag;
    local_grad(:,1) = local_grad_mag_f.*e_CF(:,1);
    local_grad(:,2) = local_grad_mag_f.*e_CF(:,2);
    local_grad(:,3) = local_grad_mag_f.*e_CF(:,3);
    
    local_avg_grad_mag = dot(grad_f',e_CF')';
    local_avg_grad(:,1)=local_avg_grad_mag.*e_CF(:,1);
    local_avg_grad(:,2)=local_avg_grad_mag.*e_CF(:,2);
    local_avg_grad(:,3)=local_avg_grad_mag.*e_CF(:,3);
    grad_f = grad_f - local_avg_grad + local_grad;
else
    disp([theInterpolationScheme;grad;phi;mdot_f]);
    exit;
end

注意，面插值格式theInterpolationScheme若为"Average:Corrected"，则会使用更为准确的修正方式来计算面梯度值，即沿着CF方向的修正，即本章第4部分所讲内容，跟如下式子是完全对应的。
$$\nabla {\varphi }_f=\overline{\nabla {\varphi }_f}+\left[\frac{{\varphi }_F-{\varphi }_C}{d_{CF}}-(\overline{\nabla {\varphi }_f}\cdot {\vec{e}}_{CF})\right]{\vec{e}}_{CF}$$
其中
$$\overline{\nabla {\varphi }_f}=g_C\nabla {\varphi }_C+g_F\nabla {\varphi }_F$$
$${\vec{e}}_{CF}=\frac{{\vec{d}}_{CF}}{d_{CF}}$$
$${\vec{d}}_{CF}={\vec{r}}_F-{\vec{r}}_C$$

5.2 OpenFOAM

to be continued…