1. 源起

最近想用有限差分法計算二維的頂蓋驅動流，在推導過程中遇到了許多問題，例如對流項的離散、“線性化”這一在有限體積法中很常見的操作、\(1-\delta\)壓差在有限差分法中的實現、壓力項的離散等等。以上問題讓我感覺自己理論方面實在欠缺，所以還是從最簡單的問題練手。

這個系列來自B站Up主“Red-Green鯉魚”，視訊連結

2. 一維非穩態對流擴散方程

\[\frac{\partial f}{\partial t} + U\frac{\partial f}{\partial x} = D \frac{\partial^2 f}{\partial x^2} \]

2.1. 常係數：U和D為定值

方程特解：

\[f(x,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(x-1-Ut)^2}{D(4t+1)} \right) \]

2.2. 問題描述

用有限差分法求解上述方程，其中\(U=1.0, D=0.05, x\in[0,9]\)，並將數值解與\(t=2.5\)時的解析解對比

2.2.1. 初始條件

\[f(x,0)=exp \left( -\frac{(x-1)^2}{D} \right) \]

2.2.2. 邊界條件

直接給出邊界的值，為Dirichlet邊界條件

\[f(0,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(1+Ut)^2}{D(4t+1)} \right) \] \[f(9,t)=\frac{1}{\sqrt{4t+1}}exp \left( -\frac{(8-Ut)^2}{D(4t+1)} \right) \]

3. 方程離散和區域離散

3.1. 區域離散

有限差分法中，將連續的區域離散為一系列結點，每個結點上含有物理量的值。假設結點間等距，距離為\(\delta\)，那麼離散方式可以分為內外兩種：

1. 結點和邊界重合，為外結點法
2. 結點和邊界之間相差\(\delta /2\)，為內結點法

以下將比較兩種結點設定的方式。

1. 外結點法：設定\(N+1\)個結點，編號從0~N，結點間距\(\delta = L/N\)，結點\(i\)的座標為\(x_i=i\cdot \delta\)
2. 內結點法：設定\(N\)個結點，編號從1~N，結點間距\(\delta = L/N\)，結點\(i\)的座標為\(x_i=i\cdot \delta - \delta /2\)

3.2. 方程離散

3.2.1. 非穩態項

對於結點\(i\)，採用一階尤拉

\[\frac{\partial f}{\partial t}\bigg|_{i}=\frac{f^{n+1}_i-f^n_i}{\Delta t} \]

3.2.2. 對流項

一階迎風，由於給出對流項前置係數\(U=1.0\)，因此迎風格式為

\[\frac{\partial f}{\partial x}\bigg|_{i}=\frac{f_i-f_{i-1}}{\delta} \]

3.2.3. 擴散項

二階中心差分

\[\frac{\partial^2 f}{\partial x^2}\bigg|_{i}=\frac{f_{i+1}-2f_i+f_{i-1}}{\delta ^2} \]

注意：以上離散格式均是針對內結點。

3.3. 差分方程

將以上各項的離散形式組合起來得到差分方程。需要注意的是，如果對流項和擴散項中用到的物理量的值為下一時刻即\(n+1\)時間步的值，則為隱式離散；如果為當前時刻即\(n\)時間步的值，則為顯式離散。

為了表示方便，下文中\(i\)用\(P\)表示，\(i-1\)用\(W\)表示，\(i+1\)用\(E\)表示。

3.3.1. 顯式

\[f^{n+1}_P=a_P f_P + a_W f_W + a_E f_E \] \[a_P=1-\frac{U\Delta t}{\delta} -\frac{2D\Delta t}{\delta^2} \] \[a_W=\frac{U\Delta t}{\delta} + \frac{D\Delta t}{\delta^2} \] \[a_E=\frac{D\Delta t}{\delta^2} \]

根據Harten定理係數非負原則，該顯式格式穩定的條件是\(a_P\)、\(a_W\)、\(a_E\)同時大於0，可見，這三項中只有\(a_P\)存在小於0的可能，通過調整程式中\(\delta\)和\(\Delta t\)能夠發現這一現象。

3.3.2. 隱式

\[a_P f^{n+1}_P=a_W f^{n+1}_W + a_E f^{n+1}_E + b \] \[a_P=1+\frac{U\Delta t}{\delta}+\frac{2D\Delta t}{\delta^2} \] \[a_W=\frac{U\Delta t}{\delta} + \frac{D\Delta t}{\delta^2} \] \[a_E=\frac{D\Delta t}{\delta^2} \] \[b=f^n_P \]

3.4. 邊界離散

3.4.1. 外結點法

結點編號從0到N，共N+1個結點，對於結點\(i\in [1,N-1]\)，均採用上述離散方程求解，對於\(i=0,N\)的兩個結點，由Dirichlet邊界條件直接獲得結點上的值。

3.4.2. 內結點法

結點編號從1到N，共N個結點，結點\(i\in [2,N-1]\)按照上述離散方程求解。結點\(i=1,N\)的對流項以及擴散項分別為

3.4.2.1. 對流項

\[\frac{\partial f}{\partial x}\bigg|_{i=1}=\frac{f_1-f_{x=0}}{\delta/2},\quad \frac{\partial f}{\partial x}\bigg|_{i=N}=\frac{f_N-f_{N-1}}{\delta} \]

3.4.2.2. 擴散項：

對於左邊界：

\[f_{x=0}=f_0=f_1-\frac{\delta}{2}f^{\prime}+\frac{\delta^2}{8}f^{\prime\prime} + \varepsilon(\delta^3) \\ f_2=f_1+\delta f^{\prime} + \frac{\delta^2}{2} f^{\prime\prime} +\varepsilon(\delta^3) \\ \]

消去\(f^{\prime}\)可得到二階導數

\[\frac{\partial^2 f}{\partial x^2}\bigg|_{i=1}=\frac{2f_0+f_{2}-3f_1}{\delta ^2 (3/4)}+\varepsilon(\delta) \]

注意：上式只有一階精度。同理，右邊界按類似的方式處理，令右邊界為\(N+1\)號結點，則結點\(N\)

\[\frac{\partial^2 f}{\partial x^2}\bigg|_{i=N}=\frac{2f_{N+1}+f_{N-1}-3f_N}{\delta ^2 (3/4)}+\varepsilon(\delta) \]

使用一階精度擴散項的顯式演算法：

\[f^{n+1}_1=f_1-\frac{2U\Delta t}{\delta}(f_1-f_0)+\frac{4D\Delta t}{3\delta^2}(2f_0+f_2-3f_1) \] \[f^{n+1}_N=f_N-\frac{U\Delta t}{\delta}(f_N-f_{N-1})+\frac{4D\Delta t}{3\delta^2}(2f_{N+1}+f_{N-1}-3f_N) \]

使用一階精度擴散項的隱式演算法：

\[\left( 1+\frac{2U\Delta t}{\delta}+\frac{4D\Delta t}{\delta ^2} \right)f_1^{n+1}=\frac{4D\Delta t}{3\delta^2}f_2^{n+1}+ \left( \frac{2U\Delta t}{\delta}+\frac{8D\Delta t}{3\delta^2} \right)f_0^{n+1} + f_1^n \] \[\left( 1+\frac{U\Delta t}{\delta}+\frac{4D\Delta t}{\delta^2} \right)f_N^{n+1}=\left( \frac{U\Delta t}{\delta}+\frac{4D\Delta t}{3\delta^2} \right)f_{N-1}^{n+1}+\frac{8D\Delta t}{3\delta^2}f_{N+1}^{n+1} + f_N^n \]

上式中的\(f_0^{n+1}\)以及\(f_{N+1}^{n+1}\)均為已知的邊界條件。

3.4.2.3. 擴散項二階精度的補充結點法

這種方法來自《數值傳熱學》例4-4

對於內結點法，\(i=1\)的結點和左邊界距離為\(\delta /2\)，我們在左邊界上設定一個結點，編號為\(i=0\)，在左邊界往左\(\delta/2\)處補充結點\(i=-1\)，對於結點\(i=1\)

\[\frac{\partial^2 f}{\partial x^2}\bigg|_{i=1}=\frac{f_{2}+f_{-1}-2f_1}{\delta ^2}+\varepsilon(\delta ^2) \]

對於結點\(i=0\)

\[\frac{\partial^2 f}{\partial x^2}\bigg|_{i=0}=\frac{f_{1}+f_{-1}-2f_{0}}{(\delta/2) ^2}+\varepsilon(\delta ^2) \]

由於結點0處於邊界，在滿足邊界條件的同時應儘量讓該結點滿足控制方程，所以對結點0的控制方程的顯式離散為

\[\frac{f_0^{n+1}-f_0^n}{\Delta t}+U\frac{f_0-f_{-1}}{\delta /2}=D\frac{f_{1}+f_{-1}-2f_{0}}{(\delta/2) ^2} \]

\(f_0^{n+1}\)以及\(f_0^n\)為Dirichlet邊界條件，均已知，因此可求出\(f_{-1}\)，將\(f_{-1}\)代入結點\(i=1\)的擴散項差分方程中，最終可得到滿足空間二階精度的擴散項差分方程。

4. 代數方程組求解方法

顯式離散能夠直接推匯出物理量的更新方程，因此只要給出\(\Delta x\)和\(\Delta t\)便能求解，但因為穩定性條件（如Harten定理）限制了所能採用的\(\Delta t\)

隱式離散的穩定性更好，但是需要聯立求解方程組。空間被離散為\(N\)個結點，每個結點都對應一個方程，即\(N\)個方程\(N\)個未知數。方程組求解的快慢直接決定了整個計算流程的效率，因此研究者開發了若干方程組求解方法，可分為兩大類，直接求解和迭代求解

4.1. 直接求解

當求解的未知數個數極少時，採用Cramer法則直接求解。本文討論的一維問題中，每個結點只和周圍兩個結點相關，即\(a_P\)只和\(a_W\)以及\(a_E\)相關，因此，求解的代數方程組只有三個對角上的元素非零，即三對角矩陣。下面介紹針對該矩陣發展的演算法，基於Gauss消元法的Thomas演算法：TDMA

4.1.1. TDMA: Tridiagonal matrix method

為了節省篇幅，以下只討論區域離散方式為外結點法的代數方程組求解。事實上有限體積法的離散方式更像內結點法，因為體心值作為網格平均值有二階精度（參考連結）。

矩陣形式為\(AX=B\)，由於使用了第一類邊界條件，\(x_1\)和\(x_N\)已知，因此\(a_{12}=a_{N,N-1}=0\)

\[\left[ \begin{matrix} a_{11} & a_{12} & & & & \\ a_{21} & a_{22} & a_{23} & & & & \\ & & & \cdots \\ & & a_{i,i-1} & a_{ii} & a_{i,i+1} & & \\ & & & \cdots \\ & & & & a_{N-1,N-2} & a_{N-1,N-1} & a_{N-1,N} \\ & & & & & a_{N,N-1} & a_{N,N} \\ \end{matrix} \right] \left[ \begin{matrix} x_{1} \\ x_{2} \\ \cdots \\ x_{i} \\ \cdots \\ x_{N-1} \\ x_{N} \end{matrix} \right]= \left[ \begin{matrix} b_{1} \\ b_{2} \\ \cdots \\ b_{i} \\ \cdots \\ b_{N-1} \\ b_{N} \end{matrix} \right] \]

TDMA演算法包含兩個步驟，消元和回代。消元指的是通過將上一行整體加到下一行從而消去下一行最左邊的元素，從矩陣形狀上來看就是將三對角矩陣變成了上對角矩陣，此做法的目的是將最後一行的\(a_{N,N-1}\)消去從而能直接計算出\(x_N\)，得到\(x_N\)後便逐一向上回代解出其他未知數。

為了討論方便，把通式寫成\(A_iT_i=B_iT_{i+1}+C_iT_{i-1}+D_i\)，消元的目的是把該式化為\(T_{i-1}=P_{i-1}T_i+Q_{i-1}\)，通過整理可以得到係數\(P_i,Q_i\)和\(B_i,C_i,D_i\)之間的關係，

\[P_i=\frac{B_i}{A_i-C_iP_{i-1}},\quad Q_i=\frac{D_i+C_iQ_{i-1}}{A_i-C_iP_{i-1}} \]

可以看出，\(P_i\)和\(Q_i\)都需要遞迴求解

\[P_1=\frac{B_1}{A_1},\quad Q_1=\frac{D_1}{A_1},\quad Q_N=T_N \]

綜上，TDMA計算流程為

1. 計算\(P_1,Q_1\)
2. 更新\(P_i,Q_i\)，得到\(P_N,Q_N\)，其中\(Q_N=T_N\)
3. 由\(T_{i-1}\)的表示式計算\(T_1-T_{N-1}\)

4.2. 迭代求解

待補充

5. 計算結果討論

5.1. 顯式格式

由上述離散過程可見，對流項採用一階迎風。一階格式的結果是截斷項首項為二階導數，而二階導數項有擴散性質，即物理量會向各個方向傳播。通常這種由二階導數引起的擴散現象稱為假擴散/人工粘性/數值粘性。具體討論可參考《數值傳熱學》第二版5.5節。

下圖為顯式離散在\(t=2.5s\)的計算結果，區域離散方式為外結點和內結點的差別不大。結點間距\(\Delta x=0.009\)，時間步\(\Delta t=0.0001s\)，該步長滿足Harten定理。可見，模擬結果與解析解相比，有矮胖特徵，意味著存在假擴散現象。

5.2 隱式格式

下圖為隱式離散在\(t=2.5s\)的計算結果，區域離散採用外結點，結點間距0.009，時間步\(\Delta t=0.001s\)，而在相同時間步下，顯式格式的計算結果會發散，比較來看，隱式格式能夠適用更大的時間步長。