計算流體力學——從原理到代碼（四）：有限體積法與Godunov Scheme 初步

05-23

來自專欄第三新東京市

太長不看

FVM的基本思想：

對於位置 $i$ 時刻 $n$ ，即每個格點上的數值解 $U_i^n$ ，都是x軸區間 $[x_{i-frac{1}{2}} , x_{i+frac{1}{2}}]$ 上的體積平均值： $U _ { i } ^ { n } = frac { 1} { Delta x } int _ { x _ { i - 1/ 2} } ^ { x _ { i + 1/ 2} } u left( x ,t _ { n } ight) d x$
每個x軸區間 $I_i = [x_{i-frac{1}{2}} , x_{i+frac{1}{2}}]$ 的交界處的流量數值值，時間t區間上的平均值， $F_{i-frac{1}{2}}^{n+frac{1}{2}} = frac { 1} { Delta t } int _ { t _ { n } } ^ { t _ { n + 1} } f left( u left( x _ { i - 1/ 2} ,t ight) ight) d t$
於是，更新下一個時間n+1上的數值解 $U_i^{n+1}$ , 即 $U _ { i } ^ { n + 1} = U _ { i } ^ { n } + frac { Delta t } { Delta x } left( F _ { i + frac{1}{2}} ^ { n + frac{1}{2}} - F _ { i - frac{1}{2}} ^ { n + frac{1}{2}} ight)$

Godunov Method的基本思想

在每個單元格點定義體積平均值 $U_i^n$ ，重建分段多項式函數 $ilde { u } ^ { n } ( x )$ （在 $[x_{i-1/2}, x_{i+1/2}]$ 上定義）（注意，上標n不是指數！代表的是時間為n的時候的分段函數！！）
最簡單的情況（分段常數函數）： $ilde { u } ^ { n } ( x ) = U _ { i } ^ { n } ext{ for } x _ { i - 1/ 2} leq x < x _ { i -1 / 2}$
使用n-階分段多項式函數 $ilde { u } ^ { n } ( x )$ 作為初始條件，在一個時間單位上 $Delta t$ ，對雙曲方程進行精確（或近似）來獲得下一步的解 $ilde { u } ^ { n + 1} ( x )$ 。
在每個單元上求平均值 $U _ { i } ^ { n + 1} = frac { 1} { Delta x } int _ { x _ { i - 1/ 2} } ^ { x _ { i + 1/ 2} } ilde { u } ^ { n + 1} ( x ) d x$ 以獲得新的單元平均值

一些老生常談的東西（了解的直接跳過好了）

CFD最困難的不是理論，也不是代碼，而是「辭彙」，即各種符號與術語，以及之間的關係。所以，請務必要清楚我們說的所謂「解」，「流量」，到底是什麼。
首先，時間-空間：t-x ，然後離散化，空間-時間的格點分別記為 $x_i, t_n$
小寫u，表示真解。大寫加上下標的 $U_i^n$ ，表示時空點 $(iDelta x, n Delta t )$ 上的數值解
解是什麼？一般來說，我們研究的u包含流體體積微元上的質量ρ，動量ρu（u代表速度）和能量E。
流量 F，是u的函數，一般來說表示單位時間內單位體積上解的進出值，一般由方程給出。
帶上下標的流量 $F_i^{n}$ 也是數值流量的意思。

一、流體方程的積分形式及弱解

之前，我們一直研究的是流體守恆方程的微分形式: $\ frac{partial }{partial { t }} u + A(u) abla_x u = 0$

上一節，我們討論了非線性的流體問題，例如，Burgers equation:

$frac{partial}{partial t} u + frac{ partial u^2/2}{partial x} =0 \Leftrightarrow frac{partial}{partial t} u + ufrac{partial}{partial x} u = 0$

在 $x=0$ 處的Riemann問題會產生解的非連續性。而非連續的地方，PDE則失效了。

但是這樣的解，可以用積分形式來表示——因為積分形式的方程是流體守恆的更基本的描述（之後會補充流體基本方程的推導過程，或者大家也可以參考z站上的其他答案）

$int _ { x _ { 1} } ^ { x _ { 2} } left[ u left( x ,t _ { 2} ight) - u left( x ,t _ { 1} ight) ight] d x = int _ { t _ { 1} } ^ { t _ { 2} } F left[ u left( x _ { 1} ,t ight) ight] - F left[ u left( x _ { 2} ,t ight) ight] d t\ ext{ for any } t _ { 1} ,t _ { 2} ,x _ { 1} ,x _ { 2}$

這樣的方程的解 $u(x,t)$ 被稱為弱解。有關弱解的唯一性暫時不深聊。

二、有限體積法思想

相比於差分格式，基於的是微分形式的方程，以至於處理初始條件間斷處的時候產生了不希望產生的不穩定效應。

所以，有限體積法 FInite Volume Method 的關鍵，就是為了處理弱解，而轉去基於方程的積分形式，從而來進行離散進而逼近真解。

對於位置 $i$ 時刻 $n$ ，即每個格點上的數值解 $U_i^n$ ，都是x軸區間 $[x_{i-frac{1}{2}} , x_{i+frac{1}{2}}]$ 上的體積平均值： $U _ { i } ^ { n } = frac { 1} { Delta x } int _ { x _ { i - 1/ 2} } ^ { x _ { i + 1/ 2} } u left( x ,t _ { n } ight) d x$
每個x軸區間 $I_i = [x_{i-frac{1}{2}} , x_{i+frac{1}{2}}]$ 的交界處的流量數值值，是時間t區間上的平均值， $F_{i-frac{1}{2}}^{n+frac{1}{2}} = frac { 1} { Delta t } int _ { t _ { n } } ^ { t _ { n + 1} } f left( u left( x _ { i - 1/ 2} ,t ight) ight) d t$
於是，更新下一個時間n+1上的數值解 $U_i^{n+1}$ , 即 $U _ { i } ^ { n + 1} = U _ { i } ^ { n } + frac { Delta t } { Delta x } left( F _ { i + frac{1}{2}} ^ { n + frac{1}{2}} - F _ { i - frac{1}{2}} ^ { n + frac{1}{2}} ight)$

那麼，如果我們只知道體積平均值，即數值解 $U_i^n$ 。為了獲得每個單元 $I_i^n$ 之間界面的通量，本質上就是需要知道解在 $x_{i pm 1/2}$ 的值。一種方法，通過某種對U的平均/插值來估計，或者另一種方法，找到直接逼近界面通量的方法。這是有限體積法的關鍵。

三、逼近界面流量：線性對流方程情形

我們先從最簡單的出發——

對於線性平流方程，來給出上述有限差分方法的解釋：

先回顧一下，對於方程 $u_t + Au_x = 0$

1、一階迎風格式（Upwind Scheme）如下：

$frac { U _ { i } ^ { n + 1} - U _ { i } ^ { n } } { Delta t } = left{ egin{array} { l } { - A left( U _ { i + 1} ^ { n } - U _ { i } ^ { n } ight) / Delta x } & { ( A < 0) } \ { - A left( U _ { i } ^ { n } - U _ { i - 1} ^ { n } ight) / Delta x } & { ( A geq 0) } end{array} ight.$

so，界面的流量即如下：

$F _ { i - 1/ 2} = left{ egin{array} { l l } { A U _ { i - 1} } & { ( A geq 0) } \ { A U _ { i } } & { ( A < 0) } end{array} ight.$

2、Lax-Friedreich

$frac { u _ { i } ^ { n + 1} - u _ { i } ^ { n } } { Delta t } = - A left( frac { u _ { i + 1} ^ { n } - u _ { i - 1} ^ { n } } { 2Delta x } ight) + frac { Delta x ^ { 2} } { 2Delta t } left( frac { u _ { i - 1} ^ { n } - 2u _ { i } ^ { n } + u _ { i + 1} ^ { n } } { Delta x ^ { 2} } ight)$

so， $U_i^{n+1}$ 根據前部分FVM更新規則，即 $U _ { i } ^ { n + 1} = U _ { i } ^ { n } + frac { Delta t } { Delta x } left( F _ { i + frac { 1} { 2} } ^ { n + frac { 1} { 2} } - F _ { i - frac { 1} { 2} } ^ { n + frac { 1} { 2} } ight)$ ，那麼對應這個格式裡面：

$F _ { i - 1/ 2} = frac { 1} { 2} A left( U _ { i - 1} + U _ { i } ight) - frac { Delta x } { 2Delta t } left( U _ { i } - U _ { i - 1} ight)$

3、Lax-Wendroff 格式

$frac { u _ { i } ^ { n + 1} - u _ { i } ^ { n } } { Delta t } = - A left( frac { u _ { i + 1} ^ { n } - u _ { i - 1} ^ { n } } { 2Delta x } ight) + frac { A ^ { 2} Delta t } { 2} left( frac { u _ { i + 1} ^ { n } - 2u _ { i } ^ { n } + u _ { i - 1} ^ { n } } { Delta x ^ { 2} } ight)$

$F _ { i - 1/ 2} = frac { 1} { 2} A left( U _ { i - 1} + U _ { i } ight) - frac { Delta t } { 2Delta x } A ^ { 2} left( U _ { i } - U _ { i - 1} ight)$

四、Godunov Method的基本思想

在每個單元格點定義體積平均值 $U_i^n$ ，重建分段多項式函數 $ilde { u } ^ { n } ( x )$ （在 $[x_{i-1/2}, x_{i+1/2}]$ 上定義）（注意，上標n不是指數！代表的是時間為n的時候的分段函數！！）
最簡單的情況（分段常數函數）： $ilde { u } ^ { n } ( x ) = U _ { i } ^ { n } ext{ for } x _ { i - 1/ 2} leq x < x _ { i -1 / 2}$
使用n-階分段多項式函數 $ilde { u } ^ { n } ( x )$ 作為初始條件，在一個時間單位上 $Delta t$ ，對雙曲方程進行精確（或近似）來獲得下一步的解 $ilde { u } ^ { n + 1} ( x )$ 。
在每個單元上求平均值 $U _ { i } ^ { n + 1} = frac { 1} { Delta x } int _ { x _ { i - 1/ 2} } ^ { x _ { i + 1/ 2} } ilde { u } ^ { n + 1} ( x ) d x$ 以獲得新的單元平均值

這是雛形，先說到這裡，之後會花大力氣細講

五、代碼的大致框架

基本上，框架非常簡單，困難就在於你選擇的格式，Your_Scheme_Type，這個函數的內部功能，我們接下來會花大工夫來說。

總體來說：

% 設置網格，初值% 略for n=1:nstep % overall time iteration for i=1:N-1 % N is the number of X-axis grids % computing F_i+1/2 F(i,:) = Your_Scheme_Type(U_left(i,:), U_right(i,:), dt); end for i=2:N-1 % computing U_i^(n+1) using FISRT ORDER time scheme U(i,:) = U(i,:)- dt / dx * (F(i,:) - F(i-1,:)); end % higher accuracy limiter % [U_left, U_right] = Your_Limiter( U ) % otherwise U_left = U(1:N-1,:); U_right = U(2:N,:); % boundary conditions % 略end