学习自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维问题，且均匀网格，则面eee上的变量梯度可直接得出：

(∂ϕ∂x)e=ϕE−ϕCxE−xC=ϕE−ϕCδxe\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} (∂x∂ϕ)e=xE−xCϕE−ϕC=δxeϕE−ϕC
单元形心CCC处的变量梯度也可得出
(∂ϕ∂x)C=ϕe−ϕwxe−xw=ϕE+ϕC2−ϕC+ϕW2ΔxC=ϕE−ϕW2ΔxC\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} (∂x∂ϕ)C=xe−xwϕe−ϕw=ΔxC2ϕE+ϕC−2ϕC+ϕW=2ΔxCϕE−ϕW
这个常叫做“中心差分”，对于2维情况，与此类似，可得

(∂ϕ∂x)C=ϕE−ϕWxE−xW,(∂ϕ∂y)C=ϕN−ϕSxN−xS\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} (∂x∂ϕ)C=xE−xWϕE−ϕW,(∂y∂ϕ)C=xN−xSϕN−ϕS
当处理非正交网格或非结构网格时，情况就复杂得多了，这就需要用到更加通用方法，即Green-Gauss梯度法和最小二乘梯度法等。

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

这是计算梯度方法中应用最广的一个，即，单元形心处变量的梯度可以由面形心处的变量值与面积矢量复合后相加和除以单元体积来获取，用到的数学定理是Green-Gauss或梯度定理，即
∫V∇ϕdV=∮∂VϕdS⃗⇒∇ϕ‾V=∮∂VϕdS⃗⇒∇ϕC=1VC∑f−nb(C)ϕfS⃗f\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 ∫V∇ϕdV=∮∂VϕdS⇒∇ϕV=∮∂VϕdS⇒∇ϕC=VC1f−nb(C)∑ϕfSf
即
∇ϕC=1VC∑f−nb(C)ϕfS⃗f\nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_f \vec S_f ∇ϕC=VC1f−nb(C)∑ϕfSf
其中VCV_CVC为单元体积，fff代表单元的某个面，而S⃗f\vec S_fSf为该面的面积矢量，面形心的变量值ϕf\phi_fϕf是未知的，仍旧需要计算出来，不然上面这个公式是用不了的。ϕf\phi_fϕf的计算方法有两种，一是用紧致框架（compact stencil）由面左右单元（所属单元和邻居单元，即面两侧单元）形心值来计算，二是用扩展框架（extended stencil）先用面的角点周围单元上的值来得到面的角点值，然后再由角点值来平均得到面形心的值。

方法1：紧致框架

对于上图(a)和(b)所示的2维和3维情况，一种最简单的方法是直接由面两侧单元形心值平均来获得面上变量的值，即
ϕf=gCϕC+(1−gC)ϕF\phi_f=g_C\phi_C+(1-g_C)\phi_F ϕf=gCϕC+(1−gC)ϕF
其中gCg_CgC为几何权重系数，等于
gC=∣∣r⃗F−r⃗f∣∣∣∣r⃗F−r⃗C∣∣=dFfdFCg_C=\frac{||\vec r_F-\vec r_f||}{||\vec r_F-\vec r_C||}=\frac{d_{Ff}}{d_{FC}} gC=∣∣rF−rC∣∣∣∣rF−rf∣∣=dFCdFf
其中r⃗\vec rr为距离矢量，而ddd为两点之间的距离值，若该面恰好位于两单元中心的中间位置，则有
ϕf=ϕC+ϕF2\phi_f=\frac{\phi_C+\phi_F}{2} ϕf=2ϕC+ϕF
这个方法非常简便易行，然而遗憾的是，只有当线段CFCFCF和面S⃗f\vec S_fSf的交点与面S⃗f\vec S_fSf的中心重合时（即上图（a）（b）情况），该方法才能保证2阶精度，而大多数情况下，直接由上式算得的梯度精度是没法保证的。

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

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

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

选择1

点f′f'f′就选择在线段CFCFCF和面S⃗f\vec S_fSf的交点处，以n⃗\vec nn代表面积单位矢量，即，n⃗=S⃗f/∣∣S⃗f∣∣\vec n=\vec S_f/||\vec S_f||n=Sf/∣∣Sf∣∣，以e⃗\vec ee代表沿着CFCFCF的单位矢量，即e⃗=CF→/∣∣CF→∣∣\vec e=\overrightarrow{CF}/||\overrightarrow{CF}||e=CF/∣∣CF∣∣，则f′f'f′的位置可由几何关系算出，为
r⃗f′=(r⃗f−r⃗C)⋅n⃗e⃗⋅n⃗e⃗+r⃗C=r⃗f⋅n⃗e⃗⋅n⃗e⃗\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 rf′=e⋅n(rf−rC)⋅ne+rC=e⋅nrf⋅ne
其中，(r⃗f−r⃗C)⋅n⃗(\vec r_f - \vec r_C) \cdot \vec n(rf−rC)⋅n为点CCC到面S⃗f\vec S_fSf的垂直距离向量，分母e⃗⋅n⃗\vec e \cdot \vec ne⋅n为该垂直向量与CF→\overrightarrow{CF}CF夹角的余弦值cosθcos\thetacosθ，于是两者相除便得到了CCC到f′f'f′的向量Cf′→\overrightarrow{Cf'}Cf′，再与r⃗C\vec r_CrC相加便是点f′f'f′的位置矢量。
得到f′f'f′的位置后，可以直接算出gCg_CgC的值
gC=∣∣r⃗F−r⃗f′∣∣∣∣r⃗F−r⃗C∣∣=dFf′dFCg_C=\frac{||\vec r_F - \vec r_{f'}||}{||\vec r_F - \vec r_{C}||}=\frac{d_{Ff'}}{d_{FC}} gC=∣∣rF−rC∣∣∣∣rF−rf′∣∣=dFCdFf′
计算流程如下：

计算ϕf′\phi_{f'}ϕf′使用ϕf′=gCϕC+(1−gC)ϕF\phi_{f'}=g_C\phi_C+(1-g_C)\phi_Fϕf′=gCϕC+(1−gC)ϕF
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕf′S⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f∇ϕC=VC1f−nb(C)∑ϕf′Sf
接下来用下面步骤来修正梯度场
更新ϕf\phi_{f}ϕf使用ϕf=ϕf′+gC(∇ϕ)C⋅(r⃗f−r⃗C)+(1−gC)(∇ϕ)F⋅(r⃗f−r⃗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)ϕf=ϕf′+gC(∇ϕ)C⋅(rf−rC)+(1−gC)(∇ϕ)F⋅(rf−rF)
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f∇ϕC=VC1f−nb(C)∑ϕfSf
返回步骤3再迭代一次。

选择2

点f′f'f′就选择在线段CFCFCF的中点，相当于gC=0.5g_C=0.5gC=0.5，这就使得计算简单多了，其计算流程如下：

计算ϕf′\phi_{f'}ϕf′使用ϕf′=ϕC+ϕF2\displaystyle \phi_{f'}=\frac{\phi_C+\phi_F}{2}ϕf′=2ϕC+ϕF
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕf′S⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f∇ϕC=VC1f−nb(C)∑ϕf′Sf
接下来用下面步骤来修正梯度场
更新ϕf\phi_{f}ϕf使用ϕf=ϕf′+(∇ϕ)C+(∇ϕ)F2⋅(r⃗f−r⃗C+r⃗F2)\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)ϕf=ϕf′+2(∇ϕ)C+(∇ϕ)F⋅(rf−2rC+rF)
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f∇ϕC=VC1f−nb(C)∑ϕfSf
返回步骤3再迭代一次。

选择3

点f′f'f′的位置选择要保证距离ff′ff'ff′是最短的，即，ff′ff'ff′是垂直于CFCFCF的，这使得第1步迭代计算变得更加准确。f′f'f′的位置计算很简单，直接把向量Cf→\overrightarrow {Cf}Cf投影到CF→\overrightarrow {CF}CF上就妥了，即
r⃗f′=r⃗C+r⃗Cf⋅r⃗CFr⃗CF⋅r⃗CF(r⃗F−r⃗C)\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) rf′=rC+rCF⋅rCFrCf⋅rCF(rF−rC)
其计算流程如下：

计算r⃗f′\vec r_{f'}rf′使用r⃗f′=r⃗C+r⃗Cf⋅r⃗CFr⃗CF⋅r⃗CF(r⃗F−r⃗C)\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)rf′=rC+rCF⋅rCFrCf⋅rCF(rF−rC)
计算gCg_CgC使用gC=∣∣r⃗F−r⃗f′∣∣∣∣r⃗F−r⃗C∣∣\displaystyle g_C=\frac{||\vec r_F - \vec r_{f'}||}{||\vec r_F - \vec r_C||}gC=∣∣rF−rC∣∣∣∣rF−rf′∣∣
计算ϕf′\phi_{f'}ϕf′使用ϕf′=gCϕC+(1−gC)ϕF\phi_{f'}=g_C\phi_C+(1-g_C)\phi_Fϕf′=gCϕC+(1−gC)ϕF
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕf′S⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f∇ϕC=VC1f−nb(C)∑ϕf′Sf
接下来用下面步骤来修正梯度场
计算∇ϕf′\nabla\phi_{f'}∇ϕf′使用∇ϕf′=gC∇ϕC+(1−gC)∇ϕF\nabla\phi_{f'}=g_C\nabla\phi_C+(1-g_C)\nabla\phi_F∇ϕf′=gC∇ϕC+(1−gC)∇ϕF
更新ϕf\phi_{f}ϕf使用ϕf=ϕf′+∇ϕf′⋅(r⃗f−r⃗f′)\phi_f=\phi_{f'}+\nabla\phi_{f'}\cdot(\vec r_f-\vec r_{f'})ϕf=ϕf′+∇ϕf′⋅(rf−rf′)
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f∇ϕC=VC1f−nb(C)∑ϕfSf
返回步骤5再迭代一次。

例1

上图所示网格，单元形心CCC与其邻近单元形心F1F_1F1到F6F_6F6的坐标为
C(13,11),F1(4.5,9.5),F2(8,3),F3(17,3.5),F4(22,10),F5(16,20),F6(7,18)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) C(13,11),F1(4.5,9.5),F2(8,3),F3(17,3.5),F4(22,10),F5(16,20),F6(7,18)
角点n1n_1n1到n2n_2n2的坐标为
n1(9,14),n2(8,8),n3(12,5),n4(17,9),n5(17.5,14),n6(12,17)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) n1(9,14),n2(8,8),n3(12,5),n4(17,9),n5(17.5,14),n6(12,17)
变量ϕ\phiϕ在单元形心的值为
ϕC=167,ϕF1=56.75,ϕF2=35,ϕF3=80,ϕF4=252,ϕF5=356,ϕF6=151\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=167,ϕF1=56.75,ϕF2=35,ϕF3=80,ϕF4=252,ϕF5=356,ϕF6=151
C单元的邻近单元形心处变量ϕ\phiϕ的梯度值∇ϕ\nabla\phi∇ϕ为
∇ϕF1=10.5i+5.5j，∇ϕF2=4i+9j,∇ϕF3=4.5i+18j,∇ϕF4=11i+23j,∇ϕF5=21i+17j,∇ϕF6=19i+8j\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 ∇ϕF1=10.5i+5.5j，∇ϕF2=4i+9j,∇ϕF3=4.5i+18j,∇ϕF4=11i+23j,∇ϕF5=21i+17j,∇ϕF6=19i+8j
单元C的体积
VC=76V_C=76 VC=76
求解单元C形心处的梯度值∇ϕC\nabla\phi_C∇ϕC，使用如下方法：
a. Green-Gauss方法，不含修正
b. Green-Gauss方法，含skewness correction，f′f'f′选择为线段CF的中点

解法a. Green-Gauss方法求解∇ϕC\nabla\phi_C∇ϕC，不含修正
计算面形心（2维情况下就是面中心）坐标值
xf1=(xn1+xn2)/2=(9+8)/2=8.5yf1=(yn1+yn2)/2=(14+8)/2=11x_{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 xf1=(xn1+xn2)/2=(9+8)/2=8.5yf1=(yn1+yn2)/2=(14+8)/2=11
即
f1(8.5,11)f_1(8.5, 11) f1(8.5,11)
同样方法算得单元C的其余面的形心坐标
f2(10,6.5),f3(14.5,7),f4(17.25,11.5),f5(14.75,15.5),f6(10.5,15.5)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) f2(10,6.5),f3(14.5,7),f4(17.25,11.5),f5(14.75,15.5),f6(10.5,15.5)
计算面积矢量S⃗f1\vec S_{f_1}Sf1
S⃗f1=Δyi−Δxj=(yn2−yn1)i−(xn2−xn1)j=(8−14)i−(8−9)j=−6i+j\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 Sf1=Δyi−Δxj=(yn2−yn1)i−(xn2−xn1)j=(8−14)i−(8−9)j=−6i+j
同样方法算得其余面的面积矢量
S⃗f2=−3i−4jS⃗f3=4i−5jS⃗f4=5i−0.5jS⃗f5=3i+5.5jS⃗f6=−3i+3j\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 Sf2=−3i−4jSf3=4i−5jSf4=5i−0.5jSf5=3i+5.5jSf6=−3i+3j
计算插值系数gC1g_{C_1}gC1
gC1=F1f1F1f1+Cf1F1f1=(4.5−8.5)2+(9.5−11)2=4.272Cf1=(13−8.5)2+(11−11)2=4.5g_{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 gC1=F1f1+Cf1F1f1F1f1=(4.5−8.5)2+(9.5−11)2=4.272Cf1=(13−8.5)2+(11−11)2=4.5
算得
gC1=0.487g_{C_1}=0.487 gC1=0.487
同样，算得其它面的插值系数
gC2=0.427,gC3=0.502,gC4=0.538,gC5=0.492,gC6=0.455g_{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 gC2=0.427,gC3=0.502,gC4=0.538,gC5=0.492,gC6=0.455
计算面形心的变量值ϕf\phi_fϕf
ϕf1=gC1ϕC+(1−gC1)ϕF1=0.487∗167+(1−0.487)∗56.75=110.442\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 ϕf1=gC1ϕC+(1−gC1)ϕF1=0.487∗167+(1−0.487)∗56.75=110.442
同样算得其它面形心的变量值
ϕf2=91.364,ϕf3=123.674,ϕf4=206.27,ϕf5=263.012,ϕf6=158.28\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 ϕf2=91.364,ϕf3=123.674,ϕf4=206.27,ϕf5=263.012,ϕf6=158.28
计算单元形心C处的梯度值
∇ϕC=1VC∑f=16ϕfS⃗f=176{[110.442(−6i+j)+91.364(3i−4j)+123.674(4i−5j)+206.27(5i−0.5j)+263.012(3i+5.5j)+158.28(−3i+3j)]}=11.889i+12.433j\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 ∇ϕC=VC1f=1∑6ϕfSf=761{[110.442(−6i+j)+91.364(3i−4j)+123.674(4i−5j)+206.27(5i−0.5j)+263.012(3i+5.5j)+158.28(−3i+3j)]}=11.889i+12.433j

解法b. Green-Gauss方法，含skewness correction，f′f'f′选择为线段CF的中点
计算CF的中点f′f'f′处的变量值，使用公式ϕf′=(ϕC+ϕF)/2\displaystyle \phi_{f'}=(\phi_C+\phi_F)/2ϕf′=(ϕC+ϕF)/2，得
ϕf1=111.875,ϕf2=101ϕf3=123.5,ϕf4=209.5,ϕf5=261.5,ϕf6=159\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 ϕf1=111.875,ϕf2=101ϕf3=123.5,ϕf4=209.5,ϕf5=261.5,ϕf6=159
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕf′S⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f∇ϕC=VC1f−nb(C)∑ϕf′Sf，得
∇ϕC=1VC∑f=16ϕfS⃗f=176{[111.875(−6i+j)+101(3i−4j)+123.5(4i−5j)+209.5(5i−0.5j)+261.5(3i+5.5j)+159(−3i+3j)]}=11.510i+11.854j\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 ∇ϕC=VC1f=1∑6ϕfSf=761{[111.875(−6i+j)+101(3i−4j)+123.5(4i−5j)+209.5(5i−0.5j)+261.5(3i+5.5j)+159(−3i+3j)]}=11.510i+11.854j
接下来修正梯度值
计算d⃗f=r⃗f−r⃗C+r⃗F2\displaystyle \vec d_f = \vec r_f-\frac{\vec r_C+\vec r_F}{2}df=rf−2rC+rF，得
d⃗f1=−0.25i+0.75j,d⃗f2=−0.5i−0.5j,d⃗f3=−0.5i−0.25j,d⃗f4=−0.25i+j,d⃗f5=−0.25i,d⃗f6=0.5i+j\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 df1=−0.25i+0.75j,df2=−0.5i−0.5j,df3=−0.5i−0.25j,df4=−0.25i+j,df5=−0.25i,df6=0.5i+j
更新ϕf\phi_{f}ϕf使用ϕf=ϕf′+(∇ϕ)C+(∇ϕ)F2⋅(r⃗f−r⃗C+r⃗F2)\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)ϕf=ϕf′+2(∇ϕ)C+(∇ϕ)F⋅(rf−2rC+rF)，得
ϕf1=115.631,ϕf2=91.909ϕf3=115.766,ϕf4=224.113ϕf5=265.564,ϕf6=176.554\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 ϕf1=115.631,ϕf2=91.909ϕf3=115.766,ϕf4=224.113ϕf5=265.564,ϕf6=176.554
计算∇ϕC\nabla\phi_C∇ϕC使用∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f∇ϕC=VC1f−nb(C)∑ϕfSf，得
∇ϕC=11.594i+13.781j\nabla\phi_C = 11.594 \bold i + 13.781 \bold j ∇ϕC=11.594i+13.781j
再迭代一次修正过程，得
∇ϕC=11.652i+15.793j\nabla\phi_C = 11.652 \bold i + 15.793 \bold j ∇ϕC=11.652i+15.793j
如果继续迭代修正，你会发现，这个值压根不会收敛，而且会越来越大直至崩溃，也就是说，修正上一次两次就OK了，既不会浪费时间，也不会让值太离谱了。

方法2：扩展框架

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

用F1,F2,F3F_1, F_2, F_3F1,F2,F3处的值来算得n1n_1n1处的值，用F2,F3,F4F_2, F_3, F_4F2,F3,F4处的值来算得n2n_2n2处的值，最后再用n1,n2n_1,n_2n1,n2处的值来算得fff处的值，即可。

角点处的值由围绕角点的单元形心处的值加权平均算出，加权系数取为距离的倒数（距离越远，影响越小），即
ϕn=∑k=1NB(n)ϕFk∣∣r⃗n−r⃗Fk∣∣∑k=1NB(n)1∣∣r⃗n−r⃗Fk∣∣\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=k=1∑NB(n)∣∣rn−rFk∣∣1k=1∑NB(n)∣∣rn−rFk∣∣ϕFk
其中nnn代表角点，FkF_kFk代表邻近单元，NB(n)NB(n)NB(n)为围绕角点nnn的单元总数，∣∣r⃗n−r⃗Fk∣∣||\vec r_n-\vec r_{F_k}||∣∣rn−rFk∣∣为角点到邻近单元形心的距离。

角点值得到后，面形心的值也可得到，以2维问题为例，ϕf\phi_fϕf为
ϕf=ϕn1+ϕn22\phi_f=\frac{\phi_{n1}+\phi_{n2}}{2} ϕf=2ϕn1+ϕn2
紧接着，可算得单元形心的梯度∇ϕC\nabla\phi_C∇ϕC为
∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f ∇ϕC=VC1f−nb(C)∑ϕfSf
对于3维情况，面形心的值ϕf\phi_fϕf由角点值的距离加权平均计算，即
ϕf=∑k=1nb(f)ϕnk∣∣r⃗f−r⃗nk∣∣∑k=1nb(f)1∣∣r⃗f−r⃗nk∣∣\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}||}}} ϕf=k=1∑nb(f)∣∣rf−rnk∣∣1k=1∑nb(f)∣∣rf−rnk∣∣ϕnk
而单元形心梯度值的计算方法照旧不变。

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

用方法2：扩展框架计算单元中心梯度值
先计算角点与其周围单元形心的距离值dnFk=∣∣r⃗n−r⃗Fk∣∣d_{nF_k}=||\vec r_n-\vec r_{F_k}||dnFk=∣∣rn−rFk∣∣，得
dn1C=5,dn1F6=4.472,dn1F1=6.364d_{n_1C}=5, \quad d_{n_1F_6}=4.472, \quad d_{n_1F_1}=6.364dn1C=5,dn1F6=4.472,dn1F1=6.364
dn2C=5.831,dn2F1=3.808,dn2F2=5d_{n_2C}=5.831, \quad d_{n_2F_1}=3.808, \quad d_{n_2F_2}=5dn2C=5.831,dn2F1=3.808,dn2F2=5
dn3C=6.083,dn3F2=4.472,dn3F3=5.220d_{n_3C}=6.083, \quad d_{n_3F_2}=4.472, \quad d_{n_3F_3}=5.220dn3C=6.083,dn3F2=4.472,dn3F3=5.220
dn4C=4.472,dn4F3=5.5,dn4F4=5.099d_{n_4C}=4.472, \quad d_{n_4F_3}=5.5, \quad d_{n_4F_4}=5.099dn4C=4.472,dn4F3=5.5,dn4F4=5.099
dn5C=5.408,dn5F4=6.021,dn5F5=6.185d_{n_5C}=5.408, \quad d_{n_5F_4}=6.021, \quad d_{n_5F_5}=6.185dn5C=5.408,dn5F4=6.021,dn5F5=6.185
dn6C=6.083,dn6F5=5,dn6F6=5.099d_{n_6C}=6.083, \quad d_{n_6F_5}=5, \quad d_{n_6F_6}=5.099dn6C=6.083,dn6F5=5,dn6F6=5.099
用距离倒数作为权系数，由角点周围单元形心值获取角点值，即
ϕn=∑k=1NB(n)ϕFk∣∣r⃗n−r⃗Fk∣∣∑k=1NB(n)1∣∣r⃗n−r⃗Fk∣∣\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=k=1∑NB(n)∣∣rn−rFk∣∣1k=1∑NB(n)∣∣rn−rFk∣∣ϕFk
得
ϕn1=131.008,ϕn2=79.708ϕn3=87.317,ϕn4=168.416ϕn5=254.144,ϕn6=228.840\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 ϕn1=131.008,ϕn2=79.708ϕn3=87.317,ϕn4=168.416ϕn5=254.144,ϕn6=228.840
计算面形心处的变量值，用面角点值平均来获取，例如ϕf1=(ϕn1+ϕn2)/2\phi_{f1}=(\phi_{n1} + \phi_{n2})/2ϕf1=(ϕn1+ϕn2)/2，得
ϕf1=105.358,ϕf2=83.512ϕf3=127.866,ϕf4=211.280ϕf5=241.492,ϕf6=179.924\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 ϕf1=105.358,ϕf2=83.512ϕf3=127.866,ϕf4=211.280ϕf5=241.492,ϕf6=179.924
计算单元中心梯度值，用∇ϕC=1VC∑f−nb(C)ϕfS⃗f\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f∇ϕC=VC1f−nb(C)∑ϕfSf，得
∇ϕC=11.446i+11.767j\nabla\phi_C = 11.446 \bold i + 11.767 \bold j ∇ϕC=11.446i+11.767j

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

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

考虑上图，单元CCC有第1层邻近单元和第2层邻近单元，那么，如果单元形心的梯度∇ϕC\nabla\phi_C∇ϕC是精确的话，有
ϕF=ϕC+(∇ϕ)C⋅(r⃗F−r⃗C)\phi_F=\phi_C+(\nabla\phi)_C\cdot(\vec r_F - \vec r_C) ϕF=ϕC+(∇ϕ)C⋅(rF−rC)
在最小二乘法中，设法让上式算得的单元邻近单元值的加权加和值最小，即找到如下函数的最小值
GC=∑k=1NB(C){wk[ϕFk−(ϕC+(∇ϕ)C⋅r⃗CFk)]2}=∑k=1NB(C){wk[Δϕk−(Δxk(∂ϕ∂x)C+Δyk(∂ϕ∂y)C+Δzk(∂ϕ∂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⋅rCFk)]2}=k=1∑NB(C){wk[Δϕk−(Δxk(∂x∂ϕ)C+Δyk(∂y∂ϕ)C+Δzk(∂z∂ϕ)C)]2}
其中wkw_kwk为加权系数
Δϕk=ϕFk−ϕCΔxk=r⃗CFk⋅i⃗Δyk=r⃗CFk⋅j⃗Δzk=r⃗CFk⋅k⃗\Delta\phi_k=\phi_{F_k}-\phi_C \\ \quad \\ \Delta x_k=\vec r_{CF_k}\cdot\vec i \\ \quad \\ \Delta y_k=\vec r_{CF_k}\cdot\vec j \\ \quad \\ \Delta z_k=\vec r_{CF_k}\cdot\vec k Δϕk=ϕFk−ϕCΔxk=rCFk⋅iΔyk=rCFk⋅jΔzk=rCFk⋅k
GCG_CGC函数的最小值应满足如下条件
∂GC∂(∂ϕ∂x)=∂GC∂(∂ϕ∂y)=∂GC∂(∂ϕ∂z)=0\frac{\partial G_C}{\partial \displaystyle \left( \frac{\partial\phi}{\partial x} \right)} = \frac{\partial G_C}{\partial \displaystyle \left( \frac{\partial\phi}{\partial y} \right)} =\frac{\partial G_C}{\partial \displaystyle \left( \frac{\partial\phi}{\partial z} \right)} = 0 ∂(∂x∂ϕ)∂GC=∂(∂y∂ϕ)∂GC=∂(∂z∂ϕ)∂GC=0
得到了如下含三个未知量的三个方程
∑k=1NB(C){2wkΔxk[−Δϕk+Δxk(∂ϕ∂x)C+Δyk(∂ϕ∂y)C+Δzk(∂ϕ∂z)C]}=0∑k=1NB(C){2wkΔyk[−Δϕk+Δxk(∂ϕ∂x)C+Δyk(∂ϕ∂y)C+Δzk(∂ϕ∂z)C]}=0∑k=1NB(C){2wkΔzk[−Δϕk+Δxk(∂ϕ∂x)C+Δyk(∂ϕ∂y)C+Δzk(∂ϕ∂z)C]}=0\sum_{k=1}^{NB(C)}\left\{ 2w_k\Delta x_k\left[-\Delta\phi_k + \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\}=0 \\ \sum_{k=1}^{NB(C)}\left\{ 2w_k\Delta y_k\left[-\Delta\phi_k + \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\}=0 \\ \sum_{k=1}^{NB(C)}\left\{ 2w_k\Delta z_k\left[-\Delta\phi_k + \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\}=0 k=1∑NB(C){2wkΔxk[−Δϕk+Δxk(∂x∂ϕ)C+Δyk(∂y∂ϕ)C+Δzk(∂z∂ϕ)C]}=0k=1∑NB(C){2wkΔyk[−Δϕk+Δxk(∂x∂ϕ)C+Δyk(∂y∂ϕ)C+Δzk(∂z∂ϕ)C]}=0k=1∑NB(C){2wkΔzk[−Δϕk+Δxk(∂x∂ϕ)C+Δyk(∂y∂ϕ)C+Δzk(∂z∂ϕ)C]}=0
可以写成矩阵形式
[∑k=1NB(C)wkΔxkΔxk∑k=1NB(C)wkΔxkΔyk∑k=1NB(C)wkΔxkΔzk∑k=1NB(C)wkΔykΔxk∑k=1NB(C)wkΔykΔyk∑k=1NB(C)wkΔykΔzk∑k=1NB(C)wkΔzkΔxk∑k=1NB(C)wkΔzkΔyk∑k=1NB(C)wkΔzkΔzk][(∂ϕ∂x)C(∂ϕ∂y)C(∂ϕ∂z)C]=[∑k=1NB(C)wkΔxkΔϕk∑k=1NB(C)wkΔykΔϕk∑k=1NB(C)wkΔzkΔϕk]\begin{bmatrix} \displaystyle \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 \\ \displaystyle \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 \\ \displaystyle \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{bmatrix} \begin{bmatrix} \displaystyle \left(\frac{\partial\phi}{\partial x}\right)_C \\ \\ \displaystyle \left(\frac{\partial\phi}{\partial y}\right)_C \\ \\ \displaystyle \left(\frac{\partial\phi}{\partial z}\right)_C \end{bmatrix} = \begin{bmatrix} \displaystyle \sum_{k=1}^{NB(C)}w_k\Delta x_k\Delta \phi_k \\ \displaystyle \sum_{k=1}^{NB(C)}w_k\Delta y_k\Delta \phi_k \\ \displaystyle \sum_{k=1}^{NB(C)}w_k\Delta z_k\Delta \phi_k \end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡k=1∑NB(C)wkΔxkΔxkk=1∑NB(C)wkΔykΔxkk=1∑NB(C)wkΔzkΔxk∑k=1NB(C)wkΔxkΔyk∑k=1NB(C)wkΔykΔyk∑k=1NB(C)wkΔzkΔyk∑k=1NB(C)wkΔxkΔzk∑k=1NB(C)wkΔykΔzk∑k=1NB(C)wkΔzkΔzk⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡(∂x∂ϕ)C(∂y∂ϕ)C(∂z∂ϕ)C⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡k=1∑NB(C)wkΔxkΔϕkk=1∑NB(C)wkΔykΔϕkk=1∑NB(C)wkΔzkΔϕk⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
求解上述方程组，便可得到单元形心的梯度值∇ϕC\nabla\phi_C∇ϕC，只要左端项矩阵不奇异，即存在解。wkw_kwk的选择是至关重要的，比如把它设成1就不太合适了，常见的设置方法是距离的倒数，即
wk=1∣∣r⃗Fk−r⃗C∣∣=1ΔxFk2+ΔyFk2+ΔzFk2w_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}} wk=∣∣rFk−rC∣∣1=ΔxFk2+ΔyFk2+ΔzFk21
也可以用距离的nnn次幂的倒数来做加权系数，即
wk=1∣∣r⃗Fk−r⃗C∣∣nw_k=\frac{1}{||\vec r_{F_k}-\vec r_C||^n} wk=∣∣rFk−rC∣∣n1
幂次nnn可以是1,2,3等。

4 梯度插值到面

同样地，由面两侧的单元值插值获取，然后再做修正，修正的原则是让面梯度沿着CFCFCF方向的值与由CCC和FFF处ϕ\phiϕ值定义的局部梯度相等，如此，面梯度∇ϕf\nabla\phi_f∇ϕf为
∇ϕf=∇ϕf‾+[ϕF−ϕCdCF−(∇ϕf‾⋅e⃗CF)]e⃗CF\nabla\phi_f=\overline{\nabla\phi_f}+\left[\frac{\phi_F-\phi_C}{d_{CF}}-(\overline{\nabla\phi_f}\cdot\vec e_{CF})\right]\vec e_{CF} ∇ϕf=∇ϕf+[dCFϕF−ϕC−(∇ϕf⋅eCF)]eCF
其中
∇ϕf‾=gC∇ϕC+gF∇ϕF\overline{\nabla\phi_f}=g_C\nabla\phi_C+g_F\nabla\phi_F ∇ϕf=gC∇ϕC+gF∇ϕF
e⃗CF=d⃗CFdCF\vec e_{CF}=\frac{\vec d_{CF}}{d_{CF}} eCF=dCFdCF
d⃗CF=r⃗F−r⃗C\vec d_{CF}=\vec r_F - \vec r_C dCF=rF−rC

补充一些个人领悟：这里从单元梯度计算面梯度的过程，其思想是这样子的，在e⃗CF\vec{e}_{CF}eCF方向上的梯度由ϕF−ϕCdCFe⃗CF\displaystyle\frac{\phi_F-\phi_C}{d_{CF}}\vec e_{CF}dCFϕF−ϕCeCF给出；而在垂直于e⃗CF\vec{e}_{CF}eCF方向上的梯度则由原本两侧owner和neighbour单元（所属和邻近单元）的梯度值插值到面上的梯度值∇ϕf‾=gC∇ϕC+gF∇ϕF\overline{\nabla\phi_f}=g_C\nabla\phi_C+g_F\nabla\phi_F∇ϕf=gC∇ϕC+gF∇ϕF在垂直于e⃗CF\vec{e}_{CF}eCF方向上的分量来给出，即∇ϕf‾−(∇ϕf‾⋅e⃗CF)e⃗CF\overline{\nabla\phi_f}-(\overline{\nabla\phi_f}\cdot\vec e_{CF})\vec e_{CF}∇ϕf−(∇ϕf⋅eCF)eCF；如此一来，把e⃗CF\vec{e}_{CF}eCF方向和垂直于e⃗CF\vec{e}_{CF}eCF方向这两个方向上梯度值复合在一起，便是面上的梯度值了，即
∇ϕf=∇ϕf‾+[ϕF−ϕCdCF−(∇ϕf‾⋅e⃗CF)]e⃗CF\nabla\phi_f=\overline{\nabla\phi_f}+\left[\frac{\phi_F-\phi_C}{d_{CF}}-(\overline{\nabla\phi_f}\cdot\vec e_{CF})\right]\vec e_{CF} ∇ϕf=∇ϕf+[dCFϕF−ϕC−(∇ϕf⋅eCF)]eCF

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个，则提取网格信息theMeshtheMesh = cfdGetMesh;
endtheSize = 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=iFacesphiGrad(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:theNumberOfElementstheElement = fvmElements(iElement);phiGrad(:,iElement) = phiGrad(:,iElement)/theElement.volume;
end
%-----------------------------------------------------
% Set boundary Gradient equal to associated element
% Gradient
%-----------------------------------------------------
for iBPatch = 1:theNumberOfBElementsiBElement = 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;
elsedisp([theInterpolationScheme;grad;phi;mdot_f]);exit;
end

注意，面插值格式theInterpolationScheme若为"Average:Corrected"，则会使用更为准确的修正方式来计算面梯度值，即沿着CF方向的修正，即本章第4部分所讲内容，跟如下式子是完全对应的。
∇ϕf=∇ϕf‾+[ϕF−ϕCdCF−(∇ϕf‾⋅e⃗CF)]e⃗CF\nabla\phi_f=\overline{\nabla\phi_f}+\left[\frac{\phi_F-\phi_C}{d_{CF}}-(\overline{\nabla\phi_f}\cdot\vec e_{CF})\right]\vec e_{CF} ∇ϕf=∇ϕf+[dCFϕF−ϕC−(∇ϕf⋅eCF)]eCF
其中
∇ϕf‾=gC∇ϕC+gF∇ϕF\overline{\nabla\phi_f}=g_C\nabla\phi_C+g_F\nabla\phi_F ∇ϕf=gC∇ϕC+gF∇ϕF
e⃗CF=d⃗CFdCF\vec e_{CF}=\frac{\vec d_{CF}}{d_{CF}} eCF=dCFdCF
d⃗CF=r⃗F−r⃗C\vec d_{CF}=\vec r_F - \vec r_C dCF=rF−rC

5.2 OpenFOAM

to be continued…