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$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}$(∂x∂ϕ​)e​=xE​−xC​ϕE​−ϕC​​=δxe​ϕE​−ϕC​​
单元形心$C$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}$(∂x∂ϕ​)C​=xe​−xw​ϕe​−ϕw​​=ΔxC​2ϕE​+ϕC​​−2ϕC​+ϕW​​​=2ΔxC​ϕE​−ϕW​​
这个常叫做“中心差分”，对于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}$(∂x∂ϕ​)C​=xE​−xW​ϕE​−ϕW​​,(∂y∂ϕ​)C​=xN​−xS​ϕN​−ϕ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$∫V​∇ϕdV=∮∂V​ϕdS⇒∇ϕ​V=∮∂V​ϕdS⇒∇ϕC​=VC​1​f−nb(C)∑​ϕf​Sf​
即
$\nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_f \vec S_f$∇ϕC​=VC​1​f−nb(C)∑​ϕf​Sf​
其中$V_C$VC​为单元体积，$f$f代表单元的某个面，而$\vec S_f$Sf​为该面的面积矢量，面形心的变量值$\phi_f$ϕf​是未知的，仍旧需要计算出来，不然上面这个公式是用不了的。$\phi_f$ϕf​的计算方法有两种，一是用紧致框架（compact stencil）由面左右单元（所属单元和邻居单元，即面两侧单元）形心值来计算，二是用扩展框架（extended stencil）先用面的角点周围角点上的值来得到面的角点值，然后再由角点值来平均得到面形心的值。

方法1：紧致框架

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

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

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

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

选择1

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

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

选择2

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

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

选择3

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

1. 计算$\vec r_{f'}$rf′​使用$\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​⋅rCF​rCf​⋅rCF​​(rF​−rC​)
2. 计算$g_C$gC​使用$\displaystyle g_C=\frac{||\vec r_F - \vec r_{f'}||}{||\vec r_F - \vec r_C||}$gC​=∣∣rF​−rC​∣∣∣∣rF​−rf′​∣∣​
3. 计算$\phi_{f'}$ϕf′​使用$\phi_{f'}=g_C\phi_C+(1-g_C)\phi_F$ϕf′​=gC​ϕC​+(1−gC​)ϕF​
4. 计算$\nabla\phi_C$∇ϕC​使用$\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f'} \vec S_f$∇ϕC​=VC​1​f−nb(C)∑​ϕf′​Sf​
   接下来用下面步骤来修正梯度场
5. 计算$\nabla\phi_{f'}$∇ϕf′​使用$\nabla\phi_{f'}=g_C\nabla\phi_C+(1-g_C)\nabla\phi_F$∇ϕf′​=gC​∇ϕC​+(1−gC​)∇ϕF​
6. 更新$\phi_{f}$ϕf​使用$\phi_f=\phi_{f'}+\nabla\phi_{f'}\cdot(\vec r_f-\vec r_{f'})$ϕf​=ϕf′​+∇ϕf′​⋅(rf​−rf′​)
7. 计算$\nabla\phi_C$∇ϕC​使用$\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$∇ϕC​=VC​1​f−nb(C)∑​ϕf​Sf​
8. 返回步骤5再迭代一次。

---

例1

上图所示网格，单元形心$C$C与其邻近单元形心$F_1$F1​到$F_6$F6​的坐标为
$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)
角点$n_1$n1​到$n_2$n2​的坐标为
$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$ϕ在单元形心的值为
$\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$∇ϕ为
$\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的体积
$V_C=76$VC​=76
求解单元C形心处的梯度值$\nabla\phi_C$∇ϕC​，使用如下方法：
a. Green-Gauss方法，不含修正
b. Green-Gauss方法，含skewness correction，$f'$f′选择为线段CF的中点

解法a. Green-Gauss方法求解$\nabla\phi_C$∇ϕ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$xf1​​=(xn1​​+xn2​​)/2=(9+8)/2=8.5yf1​​=(yn1​​+yn2​​)/2=(14+8)/2=11
即
$f_1(8.5, 11)$f1​(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)$f2​(10,6.5),f3​(14.5,7),f4​(17.25,11.5),f5​(14.75,15.5),f6​(10.5,15.5)
计算面积矢量$\vec S_{f_1}$Sf1​​
$\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
同样方法算得其余面的面积矢量
$\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
计算插值系数$g_{C_1}$gC1​​
$g_{C_1}=\frac{F_1 f_1}{F_1 f_1+Cf_1}$gC1​​=F1​f1​+Cf1​F1​f1​​
$F_1 f_1=\sqrt{(4.5-8.5)^2+(9.5-11)^2}=4.272$F1​f1​=(4.5−8.5)2+(9.5−11)2​=4.272
$Cf_1=\sqrt{(13-8.5)^2+(11-11)^2}=4.5$Cf1​=(13−8.5)2+(11−11)2​=4.5
算得
$g_{C_1}=0.487$gC1​​=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$gC2​​=0.427,gC3​​=0.502,gC4​​=0.538,gC5​​=0.492,gC6​​=0.455
计算面形心的变量值$\phi_f$ϕ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$ϕf1​​=gC1​​ϕC​+(1−gC1​​)ϕF1​​=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$ϕf2​​=91.364,ϕf3​​=123.674,ϕf4​​=206.27,ϕf5​​=263.012,ϕf6​​=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$∇ϕC​=VC​1​f=1∑6​ϕf​Sf​=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

to be continued …

---

方法2：扩展框架

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

用$F_1, F_2, F_3$F1​,F2​,F3​处的值来算得$n_1$n1​处的值，用$F_2, F_3, F_4$F2​,F3​,F4​处的值来算得$n_2$n2​处的值，最后再用$n_1,n_2$n1​,n2​处的值来算得$f$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​=k=1∑NB(n)​∣∣rn​−rFk​​∣∣1​k=1∑NB(n)​∣∣rn​−rFk​​∣∣ϕFk​​​​
其中$n$n代表角点，$F_k$Fk​代表邻近单元，$NB(n)$NB(n)为围绕角点$n$n的单元总数，$||\vec r_n-\vec r_{F_k}||$∣∣rn​−rFk​​∣∣为角点到邻近单元形心的距离。

角点值得到后，面形心的值也可得到，单元形心的梯度值随后可得出，以2维问题为例，$\phi_f$ϕf​为
$\phi_f=\frac{\phi_{n1}+\phi_{n2}}{2}$ϕf​=2ϕn1​+ϕn2​​
单元形心的梯度$\nabla\phi_C$∇ϕC​为
$\displaystyle \nabla\phi_C = \frac{1}{V_C}\sum_{f-nb(C)}\phi_{f} \vec S_f$∇ϕC​=VC​1​f−nb(C)∑​ϕf​Sf​
对于3维情况，面形心的值$\phi_f$ϕ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}||}}}$ϕf​=k=1∑nb(f)​∣∣rf​−rnk​​∣∣1​k=1∑nb(f)​∣∣rf​−rnk​​∣∣ϕnk​​​​
而单元形心梯度值的计算方法照旧不变。

---

例2

---

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

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

考虑上图，单元$C$C由第1层邻近单元和第2层邻近单元，那么，如果单元形心的梯度$\nabla\phi_C$∇ϕC​是精确的话，有
$\phi_F=\phi_C+(\nabla\phi)_C\cdot(\vec r_F - \vec r_C)$ϕF​=ϕC​+(∇ϕ)C​⋅(rF​−rC​)
在最小二乘法中，设法让上式算得的单元邻近单元值的加权加和值最小，即找到如下函数的最小值
$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}
其中$w_k$wk​为加权系数
$\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
$G_C$GC​函数的最小值应满足如下条件
$\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
得到了如下含三个未知量的三个方程
$\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
可以写成矩阵形式
$\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​Δxk​k=1∑NB(C)​wk​Δyk​Δxk​k=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​Δϕk​k=1∑NB(C)​wk​Δyk​Δϕk​k=1∑NB(C)​wk​Δzk​Δϕk​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤​
求解上述方程组，便可得到单元形心的梯度值$\nabla\phi_C$∇ϕC​，只要左端项矩阵不奇异，即存在解。$w_k$wk​的选择是至关重要的，比如把它设成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}}$wk​=∣∣rFk​​−rC​∣∣1​=ΔxFk​2​+ΔyFk​2​+ΔzFk​2​​1​
也可以用距离的$n$n次幂的倒数来做加权系数，即
$w_k=\frac{1}{||\vec r_{F_k}-\vec r_C||^n}$wk​=∣∣rFk​​−rC​∣∣n1​
幂次$n$n可以是1,2,3等。

4 梯度插值到面

同样地，由面两侧的单元值插值获取，然后再做修正，修正的原则是让面梯度沿着$CF$CF方向的值与由$C$C和$F$F处$\phi$ϕ值定义的局部梯度相等，如此，面梯度$\nabla\phi_f$∇ϕf​为
$\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​
其中
$\overline{\nabla\phi_f}=g_C\nabla\phi_C+g_F\nabla\phi_F$∇ϕf​​=gC​∇ϕC​+gF​∇ϕF​
$\vec e_{CF}=\frac{\vec d_{CF}}{d_{CF}}$eCF​=dCF​dCF​​
$\vec d_{CF}=\vec r_F - \vec r_C$dCF​=rF​−rC​

5 代码讲解