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$