由於多數偏微分方程無法求得解析解，為了解決這個問題，學者們發展出不同的數值方法來求得偏微分方程的近似解，其中一種方法稱為有限體積法(Finite Volume Method)。
透過有限體積法，我們能將偏微分方程轉換成線性代數方程式並使用電腦求解，而將偏微分方程轉換為線性代數方程式的步驟稱為離散，今天就來與各位聊聊有限體積法離散原理！
《延伸閱讀：CFD離散程序簡介，3個不可不知的離散步驟！》
傳輸方程有限體積法離散
STEP1：運用散度定理將傳輸方程積分式轉換為環積分
為了解釋有限體積法的離散概念，我們以二維傳輸方程作為例子來介紹，而傳輸方程如下：
\frac{\partial (\rho \phi)}{\partial t} + \triangledown \cdot(\rho \mathbf{v} \phi) = \triangledown \cdot (\Gamma \triangledown\phi) + Q - (1)
其中ø為純量變數，t為時間，ρ為密度，v為速度向量，Γ為擴散係數，Q為源項。若上述傳輸方程式處於穩態，意即不考慮時間項變化，可將方程式改寫成式(2)。
\triangledown \cdot(\rho \mathbf{v} \phi) = \triangledown \cdot (\Gamma \triangledown\phi) + Q - (2)
在圖1中的元素C區域內將式(2)積分，即可得到式(3)。接著透過散度定理(Divergence theorem)將式(3)對流項與擴散項的體積分轉換為環積分，即可得式(4)。
\int_{V_C}\triangledown \cdot(\rho \mathbf{v} \phi) dV= \int_{V_C}\triangledown \cdot (\Gamma \triangledown\phi) dV+ \int_{V_C}Q dV- (3)
\oint_{\partial Vc}(\rho \mathbf{v} \phi)dS = \oint_{\partial Vc} (\Gamma \triangledown\phi)dS + \int_{V_C}QdV - (4)
STEP2：控制容積的面通量總和
我們將傳輸方程中的擴散項與對流項轉換成環積分形式，環積分代表擴散項通量與對流項通量沿著封閉曲線積分。對於元素C而言，即是將各個面上的通量相加，因此可將擴散項與對流項環積分以式(5) – 式(8)表示。
J^{\phi,c} = \rho \mathbf{v}\phi - (5)
J^{\phi,D} = - \Gamma \triangledown \phi - (6)
\oint_{{\partial V_c}}J^{\phi,c} \cdot dS = {\sum\limits_{faces(V_c)}} \left ( \int_{f} (\rho \mathbf{v}\phi) \cdot dS\right ) - (7)
\oint_{{\partial V_c}}J^{\phi,D} \cdot dS = {\sum\limits_{faces(V_c)}} \left ( \int_{f} (\Gamma \triangledown \phi) \cdot dS\right ) - (8)
其中J^{\phi,c}為對流項通量，J^{\phi,D}為擴散項通量。
儘管式(7)與式(8)將環積分轉換成元素C各個面上的通量積分總和，卻依然有著積分項，對於電腦而言並不是個方便計算的形式，因此我們需要進一步離散積分項，導入高斯積分點的概念，將面通量的積分項簡化並以式(9)來表示。
\int_{f} J^{\phi} \cdot d\textbf{S} = \int_{f} (J^{\phi} \cdot \textbf{n}) dS = \sum\limits_{ip(f)}(J^{\phi} \cdot \textbf{n})_{ip} \omega_{ip} S_f - (9)
其中ip指的是積分點(integration point)，而ip(f)為積分點數量，\omega_{ip}則是權重函數。若取表面中心點作為高斯積分點，則\omega_{ip}等於1。基於式(9)的概念，可以將式(7)與式(8)表示成式(11)與式(12)，即控制容積的面通量總和。
\oint_{{\partial V_c}} \left ( \int_{f} (\rho \mathbf{v}\phi) \cdot dS\right )= {\sum\limits_{faces(V_c)}}\sum\limits_{ip(f)}(\rho \mathbf{v}\phi)_{ip} \omega_{ip} \textbf{S}_f - (11)
\oint_{{\partial V_c}} \left ( \int_{f}(-\Gamma \triangledown \phi) \cdot dS\right )= {\sum\limits_{faces(V_c)}}\sum\limits_{ip(f)}(-\Gamma \triangledown \phi)_{ip} \omega_{ip} \textbf{S}_f - (12)
STEP3：控制容積總通量線性化
求得控制容積的面通量總和之後，對於圖2的元素C而言，在不考慮源項且表面中心作為高斯積分點的情形下，半離散穩態方程式可被簡化成下列形式，如式(13)所示。
\sum\limits_{f-nb(C)}(\rho \textbf{v}\phi - \Gamma^{\phi }\triangledown \phi) \cdot \textbf{S}_f = 0- (13)
由於一般物理量的變數都是儲存於各個元素中心，如圖2中的點C、點F1至點F6。因此，需要將式(13)的面通量以元素中心變數表示，此步驟稱為通量線性化(Flux linearlization)。若元素C中心點物理量為 \phi_C且周遭元素中心點物理量為 \phi_F，則線性化後的式子可表示成式(14)的形式。
\textbf{J}_f \cdot \textbf{S}_f = FluxC_f \ \phi_C + FluxF_f \ \phi_F + FluxVf- (14)
其中FluxCf與FluxFf分別為元素C中心點與周遭元素中心點線性化後的係數項，FluxVf為非線性化項。要特別注意的是，式(14)為線性化後的通式，無論是哪種離散項與離散格式，最終都能被表示成式(14)的形式。若將總通量以各面通量總和的形式來表達，式(14)可改寫成式(15)。
\sum\limits_{f-nb(C)}(\textbf{J}_f \cdot \textbf{S}_f) = \sum\limits_{f-nb(C)}(FluxC_f \ \phi_C + FluxF_f \ \phi_F + FluxVf)- (15)
我們進一步將式(15)化簡並以更精簡的符號表示，結果如式(16) – 式(19)所示。至此完成以元素中心變數來表示離散式的通式。
a_C \phi_C + \sum\limits_{f-nb(c)}(a_F \phi F) = b_C- (16)
a_C = \sum\limits_{f-nb(c)}FluxCf- (17)
a_F = FluxFf- (18)
b_C = -\sum\limits_{f-nb(c)}FluxVf- (19)
主廚結語
本次與各位介紹穩態傳輸方程式的離散化，其基本概念為將偏微方程轉換為積分式後，運用散度定理將積分式轉換成環積分，由於環積分為沿著封閉曲線的積分，對控制容積而言，即為各表面的積分項總和。
有了各表面的通量積分項總和後，對於單一表面的通量積分而言，可利用高斯積分點將其轉換成線性代數方程式，便能求得各表面通量總和的線性代數方程。然而，由於各變數儲存於元素中心點，需要進一步將各表面的通量總和以元素中心變數來表示，整個離散流程如圖3所示。
對於初次接觸有限體積法的人而言想必有些難懂，但只要釐清每個步驟背後的意義，其實並沒有想像中困難。如果對這類文章有興趣的話，可以按讚科技雞湯Facebook，持續追蹤最新消息喔！