FVM in CFD 学习笔记_第8章_空间离散之扩散项

学习自F. Moukalled, L. Mangani, M. Darwish所著The Finite Volume Method in Computational Fluid Dynamics - An Advanced Introduction with OpenFOAM and Matlab
Chapter 8 Spatial Discretization: The Diffusion Term

本章将详细讲解由Laplace算子所表示的扩散项的空间离散方法。流体控制方程中，扩散项和对流项反映了两种大相径庭的物理现象，因此二者的离散方法也是不一样的，须要分开来说解。本章首先讲解二维矩形计算域在笛卡尔网格上的含源项的扩散方程的离散方法，而后讲解Dirichlet、Von Neumann、mixed和symmetry边界条件的添加方法。接下来，介绍下在非正交的笛卡尔网格上的离散方法，并详细讲讲非正交的结构和非结构网格上的离散方法、非正交交叉扩散，还会说起非各向同性扩散的处理方法、对高度非线性系统的欠松弛处理方法。最后，讲讲uFVM中对扩散项离散过程的代码实现。

1 在矩形域上的二维扩散

如上图所示，将简单的矩形区域用笛卡尔网格剖分，并在该网格上对下面的稳态扩散方程进行离散
$-\nabla\cdot(\Gamma^\phi\nabla\phi)=Q^\phi$
其中 $\phi$表明一个标量（好比温度、湍动能等）， $Q^\phi$为单位体积 $\phi$的生成项（源项）， $\Gamma^\phi$为扩散系数。该方程可写成更加通用的形式，定义diffusion flux扩散通量 $\bold J^{\phi,D}$为
$\bold J^{\phi,D}=-\Gamma^\phi\nabla\phi$
则原方程变为
$\nabla\cdot\bold J^{\phi,D}=Q^\phi$
根据前面章节提到的离散方式和Gauss梯度定理，可得单元 $C$的离散形式
$\sum_{f \sim nb(C)}(-\Gamma^\phi\nabla\phi)_f \cdot \bold S_f = Q_C^\phi V_C$
上式的展开形式为
$(-\Gamma^\phi\nabla\phi)_e \cdot \bold S_e + (-\Gamma^\phi\nabla\phi)_w \cdot \bold S_w + (-\Gamma^\phi\nabla\phi)_n \cdot \bold S_n + (-\Gamma^\phi\nabla\phi)_s \cdot \bold S_s = Q_C^\phi V_C$
对于图中的均匀笛卡尔网格，单元的构成面的面积矢量为
$\begin{aligned} \bold S_e &= +(\Delta y)_e \bold i=||\bold S_e|| \bold i=S_e \bold i, \quad\quad \bold S_w =-(\Delta y)_w \bold i=-||\bold S_w|| \bold i=-S_w \bold i \\ \bold S_n &=+(\Delta x)_n \bold j=||\bold S_n|| \bold j=S_n \bold j, \quad\quad \bold S_s =-(\Delta x)_s \bold j=-||\bold S_s|| \bold j=-S_s \bold j \end{aligned}$
那么东边面的扩散通量为
$\begin{aligned} J^{\phi,D}_e &= -(\Gamma^\phi\nabla\phi)_e \cdot \bold S_e \\ &= -\Gamma^\phi_e S_e \left( \frac{\partial\phi}{\partial x} \bold i+ \frac{\partial\phi}{\partial y} \bold j \right)_e \cdot \bold i \\ & = -\Gamma^\phi_e (\Delta y)_e \left( \frac{\partial\phi}{\partial x} \right)_e \end{aligned}$
该扩散通量的离散形式可写成
$J^{\phi,D}_e = FluxT_e=FluxC_e\phi_C+FluxF_e\phi_E+FluxV_e$
为了肯定系数 $FluxC_e,~FluxF_e,~FluxV_e$，须要给出变量\phi在该面两侧单元形心之间的变化规律，以便计算面上的梯度值。假设 $\phi$在两单元形心之间是线性变化的，如图

那么在面 $e$上沿着 $\bold i$方向的梯度能够写成
$\left( \frac{\partial\phi}{\partial x} \right)_e=\frac{\phi_E-\phi_C}{(\delta x)_e}$
那么面 $e$处扩散通量的最终离散形式为
$\begin{aligned} FluxT_e &=-\Gamma^\phi_e (\Delta y)_e\frac{\phi_E-\phi_C}{(\delta x)_e} \\ &=-\Gamma^\phi_e \frac{(\Delta y)_e}{(\delta x)_e} (\phi_E-\phi_C) \\ &=FluxC_e\phi_C+FluxF_e\phi_E+FluxV_e \end{aligned}$
令
$gDiff_e=\frac{(\Delta y)_e}{(\delta x)_e} = \frac{||\bold S_e||}{||\bold d_{CE}||} =\frac{S_e}{d_{CE}}$
其中 $\bold d_{CE}$为单元 $C$和 $E$形心之间的距离矢量，系数为
$\begin{aligned} FluxC_e &= \Gamma^\phi_e gDiff_e \\ FluxF_e &= - \Gamma^\phi_e gDiff_e \\ FluxV_e &= 0 \end{aligned}$
采用一样的方法，也能够推得w界面的扩散通量离散形式
$FluxT_w =FluxC_w\phi_C+FluxF_w\phi_W+FluxV_w$
其中
$\begin{aligned} FluxC_w &= \Gamma^\phi_w gDiff_w \\ FluxF_w &= - \Gamma^\phi_w gDiff_w \\ FluxV_w &= 0 \end{aligned}$
而
$gDiff_w=\frac{(\Delta y)_w}{(\delta x)_w} = \frac{||\bold S_w||}{||\bold d_{CW}||} =\frac{S_w}{d_{CW}}$
一样也可获得 $n$和 $s$面上的扩散通量离散形式
$\begin{aligned} FluxT_n &= FluxC_n\phi_C+FluxF_n\phi_N+FluxV_n \\ FluxT_s &= FluxC_s\phi_C+FluxF_s\phi_S+FluxV_s \end{aligned}$
其中
$\begin{aligned} FluxC_n &= \Gamma^\phi_n gDiff_n, \quad FluxF_n = - \Gamma^\phi_n gDiff_n, \quad & FluxV_n = 0 \\ FluxC_s &= \Gamma^\phi_s gDiff_s, \quad FluxF_s = - \Gamma^\phi_s gDiff_s, \quad & FluxV_s = 0 \end{aligned}$
而
$gDiff_n=\frac{(\Delta x)_n}{(\delta y)_n} = \frac{||\bold S_n||}{||\bold d_{CN}||} =\frac{S_n}{d_{CN}} \\ \quad \\ gDiff_s=\frac{(\Delta x)_s}{(\delta y)_s} = \frac{||\bold S_s||}{||\bold d_{CS}||} =\frac{S_s}{d_{CS}}$
将各个面的对流通量离散形式汇总代入到原方程的半离散格式，可得
$a_C\phi_C + a_E\phi_E + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=FluxF_e = - \Gamma^\phi_e gDiff_e \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_e+FluxC_w+FluxC_n+FluxC_s \\ &=-(a_E+a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-(FluxV_e+FluxV_w+FluxV_n+FluxV_s) \end{aligned}$
或者，写成紧致格式
$a_C\phi_C+\sum_{F\sim NB(C)} a_F \phi_F=b_C$
其中
$\begin{aligned} a_F &=FluxF_f = - \Gamma^\phi_f gDiff_f \\ a_C &=\sum_{f\sim nb(C)} FluxC_f \\ b_C &= Q_C^\phi V_C-\sum_{f\sim nb(C)} FluxV_f \end{aligned}$
其中的 $F$表明单元 $C$的邻近单元 $E,W,N,S$，小写 $f$表明单元 $C$的构成面（邻近面） $e,w,n,s$。

2 离散方程的解释

合适的离散方法获得的离散代数方程应该能反映原守恒方程的特性。因此推出的离散方程的系数须要知足以下两个准则。

2.1 零加和准则

即无源项的时候，各系数加和应为0。
$a_C+\sum_{F\sim NB(C)} a_F =0$
或
$a_C=-\sum_{F\sim NB(C)} a_F$
或
$\sum_{F\sim NB(C)} \frac{a_F}{a_C} =-1$

2.2 反号准则

即 $a_C$和 $a_F$是反号的，这样才能保证当 $\phi_F$增长/减少的时候， $\phi_C$是减少/增长的，从而保证了有界性boundedness。

3 边界条件

边界条件的添加相当重要，尽管控制方程相同，然而不一样的边界条件所获得的解是不一样的，对应的物理意义也不同，添加错误的边界条件会致使没法获得解，或者获得的解与原物理问题绝不相干。

扩散方程，即，拉普拉斯方程的边界条件一般有Dirichlet、Von Neumann、mixed和symmetry，边界条件是施加在边界单元上的，它们有1个或多个面位于边界上，变量 $\phi$既存储在单元形心上，也存储在面形心上。

如图， $C$表明边界单元的形心，该单元有1个边界面， $b$为该边界面的形心， $\bold S_b$为该边界面的面积矢量。跟前面同样地，也能够推出单元 $C$的离散格式
$\sum_{f \sim nb(C)} (\bold J^{\phi,D} \cdot \bold S)_f = Q_C^\phi V_C$
内部面的通量还跟以前同样离散，而边界通量的离散则是要构建跟 $\phi_C$相关的线性格式
$\begin{aligned} \bold J^{\phi,D}_b \cdot \bold S_b &=FluxT_b \\ &= -\Gamma_b^\phi (\nabla \phi)_b \cdot \bold S_b \\ &= FluxC_b\phi_b+FluxV_b \end{aligned}$
边界条件的添加，要么给定边界值 $\phi_b$，要么是给出边界通量 $\bold J^{\phi,D}_b$。不一样边界条件的添加方法以下。

3.1 Dirichlet边界条件

Dirichlet边界条件是直接给定边界上的变量值 $\phi$，即
$\phi_b=\phi_{specified}$
那么
$\begin{aligned} FluxT_b &= -\Gamma_b^\phi (\nabla \phi)_b \cdot \bold S_b \\ &= -\Gamma_b^\phi \frac{||\bold S_b||}{||\bold d_{Cb}||} (\phi_b - \phi_C) \\ &= FluxC_b\phi_C+FluxV_b \end{aligned}$
推得
$\begin{aligned} FluxC_b &= \Gamma_b^\phi gDiff_b=a_b \\ FluxV_b &= -\Gamma_b^\phi gDiff_b \phi_b = -a_b \phi_b \end{aligned}$
其中
$gDiff_b=\frac{S_b}{d_{Cb}}$
对于图中所示的 $C$单元， $a_E$系数变成0了（这个边界单元也没有 $E$单元，因此系数确定是0了），其离散方程变为
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_b+\sum_{f\sim nb(C)}FluxC_f\\ &=FluxC_b+(FluxC_w+FluxC_n+FluxC_s) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$
注意：

1. 系数 $a_b$比其它邻居系数都要大，由于 $b$距离 $C$更近，因此其对 $\phi_C$的影响更大。
2. 系数 $a_C$仍旧是全部邻居系数的加和，包括 $a_b$。这也就是说，边界单元的 $\sum_{f\sim nb(C)}|a_F|/|a_C|<1$，不能保证离散方程知足零系数加和准则，有 $|a_C|>\sum_{f\sim nb(C)}|a_F|$（严格对角占优了），从而在使用迭代方法求解线性方程的时候，是能得到收敛解的。
3. $a_b\phi_b$项产生了 $FluxV_b$，如今位于方程的右端项，是 $b_C$的一部分，由于其不含未知量。

3.2 Von Neumann边界条件

如上图，若是边界上给定的是 $\phi$的通量（或垂直于面的梯度值），那么这种边界条件就是Neumann边界条件。此时，这个给定的通量是
$-(\Gamma^\phi\nabla\phi)_b\cdot\bold i=q_b$
这样的话，通量 $\bold J_b^{\phi,D}$可得
$\bold J_b^{\phi,D}\cdot\bold S_b=-(\Gamma^\phi\nabla\phi)_b \cdot ||\bold S_b|| \bold i=q_b||\bold S_b||=FluxC_b\phi_C+FluxV_b$
其中
$\begin{aligned} FluxC_b &=0 \\ FluxV_b &=q_bS_b=q_b(\Delta y)_C \end{aligned}$
当与所使用的坐标系统方向相同的时候，通量份量为正值。

这样，单元 $C$的离散方程为：
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=\sum_{f\sim nb(C)}FluxC_f=-(a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$
注意：

1. Von Neumann边界条件生成的系数 $a_C$并非占优的（也就意味着迭代解法是不靠谱的了）。
2. 若是 $q_b$和 $S_C^\phi$为0，即边界不流入通量，且单元内部没有源项，则 $\phi_C$的值是彻底由其周围邻近单元所限定的。而在其它状况下（ $q_b$或 $S_C^\phi$不0，即有边界通量流入，或者是单元内部有源项生成的时候）， $\phi_C$则能够超过（或是低于）周围邻近值 $\phi$，这是合情合理的。举个例子，若是 $\phi$是温度，则 $q_b$表明的就是施加在边界上的热流通量，那么若是边界上有热量流入的话，则靠近边界区域的温度是会比内部区域的温度要高的。
3. 关于边界上 $\phi$值的计算，一旦单元 $\phi_C$算出来后，边界的 $\phi_b$能够用下面的式子来计算
   $\phi_b=\frac{\Gamma_b^{\phi}gDiff_b\phi_C-q_b}{\Gamma_b^{\phi}gDiff_b}$
4. Von Neumann边界条件能够认为是FVM的天然边界条件，由于当某个边界面所给定的通量为0时，离散后跟该面相关的啥项也没有了。但是，Dirichlet边界条件中，即使边界面的给定 $\phi$为0，依然在离散方程中出现了相关项。

3.3 混合边界条件

混合边界条件如上图所示，边界条件是经过对流传热系数（ $h_\infin$）和 $\phi$的周围环境值 $\phi_\infin$来给出的
$\bold J_b^{\phi,D}\cdot\bold S_b=-(\Gamma^\phi\nabla\phi)_b \cdot ||\bold S_b|| \bold i = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C$
从新写成以下形式
$-\Gamma_b^\phi S_b \left( \frac{\phi_b - \phi_C}{\delta x_b} \right) = - h_\infin(\phi_\infin-\phi_b)S_b$
从这个方程能够解出 $\phi_b$
$\phi_b=\frac{h_\infin\phi_\infin + (\Gamma_b^\phi/\delta x_b)\phi_C}{h_\infin+(\Gamma_b^\phi/\delta x_b)}$
将 $\phi_b$代入到 $\bold J_b^{\phi,D}\cdot\bold S_b = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C$中，可得
$\begin{aligned} \bold J_b^{\phi,D}\cdot\bold S_b & = -h_\infin(\phi_\infin-\phi_b)(\Delta y)_C \\ & = - \underbrace{\left[ \frac{h_\infin (\Gamma_b^\phi/\delta x_b)}{h_\infin+(\Gamma_b^\phi/\delta x_b)} S_b\right]}_{R_{eq}}(\phi_\infin-\phi_C) \\ & = FluxC_b\phi_C + FluxV_b \end{aligned}$
其中
$\begin{aligned} FluxC_b &=R_{eq} \\ FluxV_b &=-R_{eq}\phi_{\infin} \end{aligned}$
那么，对于单元 $C$的离散方程，其最终形式为
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=FluxC_b+\sum_{f\sim nb(C)}FluxC_f=FluxC_b+(FluxC_w+FluxC_n+FluxC_s) \\ b_C &= Q_C^\phi V_C-(FluxV_b+\sum_{f\sim nb(C)}FluxV_f) \end{aligned}$

3.4 对称边界条件

沿着对称边界条件，标量 $\phi$的法向通量为0，所以，对称边界条件等效于Neumann边界条件（将通量值设为0），即， $FluxC_b=FluxV_b=0$，那么直接把以前推出来的Neumann边界条件离散格式中的 $q_b$设成0就好了，即
$a_C\phi_C + a_W \phi_W + a_N \phi_N + a_S \phi_S = b_C$
其中
$\begin{aligned} a_E &=0 \\ a_W &=FluxF_w = - \Gamma^\phi_w gDiff_w \\ a_N &=FluxF_n = - \Gamma^\phi_n gDiff_n \\ a_S &=FluxF_s = - \Gamma^\phi_s gDiff_s \\ a_C &=\sum_{f\sim nb(C)}FluxC_f=-(a_W+a_N+a_S) \\ b_C &= Q_C^\phi V_C-\sum_{f\sim nb(C)}FluxV_f \end{aligned}$

4 界面扩散系数

前面的离散方程中，还须要用到扩散系数 $\Gamma_e^\phi,~\Gamma_w^\phi,~\Gamma_n^\phi,~\Gamma_s^\phi$，它们分别表明着扩散系数 $\Gamma^\phi$在 $e,~w,~n,~s$面处的值。当扩散系数 $\Gamma^\phi$随空间位置变化时（好比在可压缩流动问题中，粘性系数是温度的函数，随温度不一样粘性也不相同），扩散系数一样也是存储在单元形心 $E,~W,~N,~S...$上的，那么就须要去计算界面上的扩散系数值了。实际上，湍流粘性、热传导系数等都是非均匀分部的。固然，若是扩散系数整场是定值的话，就不用这么折腾了。

计算界面扩散系数最简单的方法莫过于直接用两侧单元值作线性插值了，即
$\Gamma_e^\phi=(1-g_e)\Gamma_C^\phi+g_e\Gamma_E^\phi$
其中插值因子 $g_e$是距离的比值，即
$g_e=\frac{d_{Ce}}{d_{Ce}+d_{eE}}$
若是是笛卡尔网格，且界面恰好位于两侧单元中间的话， $g_e=0.5$，那么 $\Gamma_e^\phi=(\Gamma_C^\phi+\Gamma_E^\phi)/2$。一样也能够定义其它插值系数
$g_w=\frac{d_{Cw}}{d_{Ce}+d_{wW}} \\ \quad \\ g_n=\frac{d_{Cn}}{d_{Cn}+d_{nN}} \\ \quad \\ g_s=\frac{d_{Cs}}{d_{Cs}+d_{sS}}$
然而，在某些状况下，这种处理方法会致使错误，好比当复合材料中的传导突变的状况下。幸运的是，有个更好的且相对简单易于实现的方法。在该方法中，界面的局部传导系数不是重点，它主要关注的是如何得到界面上的扩散通量 $\bold J^{\phi,D}$。

如图所示的1维问题，假设单元 $C$是由导热系数 $\Gamma_C^\phi$的材料构成，而单元 $E$是由导热系数为 $\Gamma_E^\phi$的材料构成，则对于C和E之间的非均匀分片，稳态1维分析（不含源项）有（界面e的任何一边的通量假设是相同的）
$\begin{aligned} \bold J_e^{\phi,D}\cdot \bold S_e &=\frac{\phi_C-\phi_e}{\displaystyle \frac{(\delta x)_{Ce}}{\Gamma_C^\phi}} = \frac{\phi_e-\phi_E}{\displaystyle \frac{(\delta x)_{eE}}{\Gamma_E^\phi}} \\ &=\frac{\phi_C-\phi_E}{\displaystyle \frac{(\delta x)_{Ce}}{\Gamma_C^\phi} + \frac{(\delta x)_{eE}}{\Gamma_E^\phi}} = \frac{\phi_C-\phi_E}{\displaystyle \frac{(\delta x)_{CE}}{\Gamma_e^\phi}} \end{aligned}$
因而推得了该片的导热系数
$\frac{(\delta x)_{CE}}{\Gamma_e^\phi} = \frac{(\delta x)_{Ce}}{\Gamma_C^\phi} + \frac{(\delta x)_{eE}}{\Gamma_E^\phi} \Rightarrow \frac{1}{\Gamma_e^\phi} = \frac{1-g_e}{\Gamma_C^\phi} + \frac{g_e}{\Gamma_E^\phi}$
当界面恰好位于 $C$和 $E$的中间时， $g_e=0.5$，上式变为
$\Gamma_e^\phi=\frac{2\Gamma_C^\phi\Gamma_E^\phi}{\Gamma_C^\phi+\Gamma_E^\phi}$
这是 $\Gamma_C^\phi$和 $\Gamma_E^\phi$的调和平均（harmonic mean），而非算术平均。

须要注意的是，对于不连续扩散系数的调和平均插值方法，仅对1维扩散问题是精确的。然而，其对高维问题的应用也颇有优点，好比流体和固体区域能够不加区分，当成一个区域来计算传热系数。

---

例1

如图热传导问题，2维矩形区域由两种材料构成，以下偏微分方程为其控制方程
$\nabla \cdot (k \nabla T) = 0$
其中 $T$表明温度，热传导系数 $k$和边界条件如图所示
a.推导出全部单元的代数方程
b. 使用Gauss-Seidel迭代方法求解系统方程得到单元 $T$值
c. 计算底部、顶部、右侧的 $T$值
d. 计算经过底部、顶部、左侧边界的净热流通量，并检查能量守恒是否知足。

解
求解过程较为繁琐，再也不详述。

---

5 非笛卡尔正交网格

如图，网格仍是正交的，可是它和 $x$与 $y$轴不是对齐的，这样的网格至关因而把笛卡尔网格转了个角度。

这种网格上的离散方程和笛卡尔网格上的离散方程式彻底同样的，边界条件也是同样的，解法也是同样的。

好比
$\nabla\cdot\bold J^{\phi,D}=Q^\phi$
半离散形式
$\sum_{f \sim nb(C)}(-\Gamma^\phi\nabla\phi)_f \cdot \bold S_f = Q_C^\phi V_C$
考虑到界面 $e$，有
$J^{\phi,D}_e = -\Gamma^\phi_e (\nabla\phi \cdot \bold n)_e S_e = -\Gamma^\phi_e \left( \frac{\partial\phi}{\partial n} \right)_e S_e$
其中
$(\nabla\phi \cdot \bold n)_e = \left( \frac{\partial\phi}{\partial n} \right)_e$
为变量 $\phi$在界面 $e$上沿着 $\bold n$方向的梯度，再次对 $\phi$使用线性假设，该梯度可得
$\left( \frac{\partial\phi}{\partial n} \right)_e=\frac{\phi_E-\phi_C}{d_{CE}}$
其它项也将得出和笛卡尔网格一样的离散方程。

6 非正交非结构网格

6.1 非正交

前面所讲的状况下，通量是与面垂直的。然而，在通常状况下，结构化曲线网格和非结构网格是非正交的，即，面矢量 $\bold S_f$和链接面两侧单元形心的矢量 $\bold {CF}$是不共线的，如上图。这样，垂直于面的梯度就没法写成 $\phi_F$和 $\phi_C$的函数，而它实际上有了个垂直于 $\bold {CF}$的份量。

这样，在正交网格上梯度沿着界面法向为
$(\nabla\phi\cdot\bold n)_f=\left(\frac{\partial\phi}{\partial n}\right)_f=\frac{\phi_F-\phi_C}{||\bold r_F - \bold r_C||}=\frac{\phi_F-\phi_C}{d_{CF}}$
由于 $\bold {CF}$和 $\bold n$（面单位法矢量）是共线的。而在非结构网格上，上面这个式子所表示的梯度方向是包含 $\phi_F$和 $\phi_C$且沿着链接点 $C$和 $F$的连线方向的（跟面法向的梯度是不一样的，它仅仅是真正的面法向梯度的一个份量而已）。

若是 $\bold e$表明沿着链接节点 $C$和 $F$单位矢量，那么
$\bold e=\frac{\bold r_F - \bold r_C}{||\bold r_F - \bold r_C||}=\frac{\bold d_{CF}}{d_{CF}}$
这样，在 $\bold e$方向的梯度能够写成
$(\nabla\phi\cdot\bold e)_f = \left(\frac{\partial\phi}{\partial e}\right)_f=\frac{\phi_F-\phi_C}{||\bold r_F - \bold r_C||}=\frac{\phi_F-\phi_C}{d_{CF}}$
这样为了得到非正交网格上的线性化通量，面矢量 $\bold S_f$应该写成两个矢量 $\bold E_f$和 $\bold T_f$的加和形式，即
$\bold S_f = \bold E_f + \bold T_f$
其中 $\bold E_f$是沿着 $\bold {CF}$方向的，以保证扩散通量中的一部分能够写成节点值 $\phi_C$和 $\phi_F$的函数，即
$\begin{aligned} (\nabla\phi)_f\cdot\bold S_f &= \underbrace{(\nabla\phi)_f\cdot\bold E_f}_{orthogonal-like~contributation} + \underbrace{(\nabla\phi)_f\cdot\bold T}_{non-orthogonal~like~contributation} \\ &=E_f\left(\frac{\partial\phi}{\partial e}\right)_f + (\nabla\phi)_f\cdot\bold T_f \\ &=E_f\frac{\phi_F-\phi_C}{d_{CF}} + (\nabla\phi)_f\cdot\bold T_f \end{aligned}$
上式中右侧第一项与正交网格上的式子同样，其表明了正交贡献量。而第2项则是因为网格的非正交性所产生的交叉扩散或非正交扩散。对 $\bold S_f$的分解能够有不一样的选择，以下所述。

6.2 最小修正方法Minimum Correction Approach

如图所示， $\bold S_f$的分解方法是要保证非正交修正尽量的小，即让 $\bold E_f$和 $\bold T_f$正交。伴随着非正交修正的增长，从 $\phi_F$到 $\phi_C$的扩散贡献将减少。矢量 $\bold E_f$是这样计算的
$\bold E_f=(\bold e \cdot \bold S_f) \bold e=(S_f \cos \theta)\bold e$

6.3 正交修正方法Orthogonal Correction Approach

这种方法让包含 $\phi_F$到 $\phi_C$的项和它们在正交网格上的贡献是一样的，即， $\bold S_f$和 $\bold E_f$的幅值是同样的，那么矢量 $\bold E_f$选择为
$\bold E_f=S_f\bold e$

6.4 过松弛方法 Over-Relaxed Approach

该方法中，包含 $\phi_F$到 $\phi_C$的项将随着非正交性的增长而增长，如上图所示。经过将 $\bold T_f$选择为垂直于 $\bold S_f$来实现，这样 $\bold E_f$是这样计算的
$\bold E_f=\left( \frac{S_f}{\cos \theta} \right)\bold e=\left( \frac{S_f^2}{S_f \cos \theta} \right) \bold e = \frac{\bold S_f \cdot \bold S_f}{\bold e \cdot \bold S_f}\bold e$