計算流體力學——從原理到代碼（一）：對流方程的格式與實現

05-23

來自專欄第三新東京市

先放最重要的，以後再慢慢修稿。

一般來說，流體方程（組）都是雙曲型偏微分方程（組），不了解的可以看這個答案：

流體力學中的歐拉方程並沒有二階偏導項，為什麼還是雙曲型的？?

www.zhihu.com

對於線性對流方程（linear advection equation）：

$\partial _ { t } u + A partial _ { x } u = 0$

解析解：

$\u ( x ,t ) = u _ { 0} ( x - A t )$

這樣的解沿著x-t平面上的直線（如下）是常數：

$\ X ( t ) = X _ { 0} + A t$

我們把這樣的線稱為方程的特徵線，特徵線上，原來的偏微分方程退化成一個常微分方程。

（證明：

$\ left.egin{aligned} frac { d } { d t } u ( X ( t ) ,t ) & = partial _ { t } u ( X ( t ) ,t ) + X ^ { prime } ( t ) partial _ { x } u ( X ( t ) ,t ) \ & = partial _ { t } u ( X ,t ) + A partial _ { x } u ( X ,t ) = 0end{aligned} ight.$

特徵線的示意圖：

理解得不太明白的可以看這個：

怎麼理解偏微分方程中的特徵線，什麼樣的線可以稱為特徵線？?

www.zhihu.com

離散化：

目的：給計算機流體的初始條件，讓程序解這個方程，得到我們想要的某時刻之後的近似解。

流體問題 -> PDE求解(連續的數學問題) -> 演算法（離通過散方法，逼近連續問題的解）-> 流體問題近似的解

最簡單的考慮：空間[0,1] ，變成[0,1]之間100，1000，10000個等距格點：0.001，0.002，…

每一個格點上計算流體的量

時間也類似，0.001s, 0.002s, ……，不斷向前

如果解 $f(x)$ 在 $x_0$ 處足夠的光滑，由Taylor Expansion： $f left( x _ { 0} + Delta x ight) = f left( x _ { 0} ight) + f _ { x } ^ { prime } left( x _ { 0} ight) Delta x + frac { Delta x ^ { 2} } { 2} f _ { x x } ^ { prime prime } left( x _ { 0} ight) + frac { Delta x ^ { 3} } { 6} f _ { x x x } ^ { prime prime prime } left( x _ { 0} ight) + dots$

那 $f(x)$ 在 $x_0$ 處的估計，結果是一階精確： $frac { f left( x _ { 0} + Delta x ight) - f left( x _ { 0} ight) } { Delta x } = f _ { x } ^ { prime } left( x _ { 0} ight) + [ frac { Delta x } { 2} f _ { x x } ^ { prime prime } left( x _ { 0} ight) + frac { Delta x ^ { 2} } { 6} f _ { x x x } ^ { prime prime prime } left( x _ { 0} ight) + dots approx f_x(x_0) + o(Delta x)$

或者，使用中心差分，結果是二階精確：

$frac { f left( x _ { 0} + Delta x ight) - f left( x _ { 0} - Delta x ight) } { 2Delta x } = f _ { x } ^ { prime } left( x _ { 0} ight) + | frac { Delta x ^ { 2} } { 6} f _ { x x x } ^ { prime prime prime } left( x _ { 0} ight) + dotsapprox f_x(x_0) + o(Delta x^2)$

所以，最直接的格式FTCS（Forward time central space）:

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

這個格式，空間上 2階精度，時間上 1階。

那我們怎麼實現呢？

對於方程 $partial _ { t } u + partial _ { x } u = 0$ ，將區間[ 0， 1]，劃分成101個格點，前後再多加2個格點，為了實現outlet邊界條件（以後會講）

首先是設置變數

% This is a demo of Forward Time Central Space for Advection Equation clearclc% Set the parametersxmin = 0;xmax = 1;N = 100;dx = (xmax-xmin)/N;tmin = 0;tmax = 0.5;dt = 0.001;v = 1;% Set the initial conditionx = linspace(xmin-dx, xmax+dx, N+3);u0 = [ones(1,5), zeros(1,103-5)];utemp = u0;uplus = u0;nstep = tmax/dt;t = 0;

然後編寫主循環的思路是：

對於時間步數循環

計算邊界點的值
根據差分格式，更新每個點的值
時間和空間+1

for n = 1:nstep % calculate the boundary condition utemp(1) = utemp(3); utemp(N+3) = utemp(N+1); % calculate the new value for i = 2:N+2 uplus(i) = utemp(i) - v*dt*(utemp(i+1) - utemp(i-1))/(2*dx); % YOU ONLY NEED TO CHANGE THIS LINE!!! end % update time step and u t = t + dt; utemp = uplus; % visualization plot(x,utemp,bo-,markerfacecolor,b); shg; pause(dt); end

運行結果

這是一個非常不穩定的格式！

所以，一切沒有想像中那麼簡單，那麼，為什麼這麼一個看似能二階精度的格式會不穩定呢？

（不想搞懂理論解釋的，可以跳過下面這一部分）

von Neumann stability analysis

用有限傅立葉級數表示離散解（嚴格來說，這僅適用於周期性域中的線性問題）：

$\ u _ { i } ^ { n } = sum a e ^ { i k ( i Delta x ) }$

下一個時刻的解： $u _ { i } ^ { n + 1} =sum A _ { k } u _ { i } ^ { n }$ ，這個 $A_k$ 表示對於第k個成分的復振幅因子

我們根據FTCS的格式 $frac { u _ { i } ^ { n + 1} - u _ { i } ^ { n } } { Delta t } = - A left( frac { u _ { i + 1} ^ { n } - u _ { i - 1} ^ { n } } { 2Delta x } ight)$ ，可以推出： $\A _ { k } = 1- ext{i} frac { A Delta t } { Delta x } sin ( k Delta x )$

我們看到，無論選擇什麼時間步， $|A_k|>1$ 。這意味著解決方案在時間上呈指數級放大，無條件不穩定！

所以，早些年，大量研究在尋找好的格式，和格式穩定性理論，其中著名的格式：

Lax-Friedrichs (LF) method

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

代碼：

uplus(i) = (utemp(i+1) + utemp(i-1))/2 - v*dt*(utemp(i+1) - utemp(i-1))/(2*dx);

（做gif太煩了！老子不幹了！）

同樣的，馮諾依曼分析一下，格式： $frac { u _ { i } ^ { n + 1} - left( u _ { i - 1} ^ { n } + u _ { i + 1} ^ { n } ight) / 2} { Delta t } = - A left( frac { u _ { i + 1} ^ { n } - u _ { i - 1} ^ { n } } { 2Delta x } ight)$ ，

我們將解視為成傅里葉級數 $u _ { i } ^ { n } = a e ^ { i k ( i Delta x ) }$ 組成

對於時間n+1步， $u _ { i } ^ { n + 1} = A _ { k } u _ { i } ^ { n }$

代入原格式，化簡得到: $A _ { k } - cos ( k Delta x ) = - ext{i} frac { A Delta t } { Delta x } sin ( k Delta x )$

所以，如果希望 $|A_k|<1$ ，則滿足 $frac { A Delta t } { Delta x } < 1$ 即可

其中， $frac { A Delta t } { Delta x }$ 就是大名鼎鼎的，由計算數學三個祖師爺名字構成的——庫朗-弗雷德里奇-李維條件是 Courant-Friedrichs-Lewy (CFL) condition。

從圖上結果，發現有耗散（就是物理中的能量耗散）

為什麼有耗散呢？

$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)$

而這個最後一項

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

帶來了耗散！

還有一種格式

Upwind method

格式：

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

該方法穩定且擴散性較小，但仍然只有一階準確。

可以使用高階空間插值方案進行改進。

這是大名鼎鼎的Godnouv的前身

我知道你們要代碼

% This is a demo of First Order Upwind for Advection Equation clearclc% Set the parametersxmin = 0;xmax = 1;N = 100;dx = (xmax-xmin)/N;tmin = 0;tmax = 0.5;dt = 0.009;v = 1;% Set the initial conditionx = linspace(xmin-dx, xmax+dx, N+3);u0 = exp( -200*(x-0.25).^2 );% u0 = [ones(1,5), zeros(1,103-5)];utemp = u0;uplus = u0;nstep = tmax/dt;t = 0;for n = 1:nstep % calculate the boundary condition utemp(1) = utemp(3); utemp(N+3) = utemp(N+1); % calculate the new value for i = 2:N+2 uplus(i) = utemp(i) - v*dt*(utemp(i) - utemp(i-1))/dx; end % update time step and u t = t + dt; utemp = uplus; % calculate exact solution exact = exp(-200*(x-0.25-v*t).^2); % visualization plot(x,exact,r-) hold on plot(x,utemp,bo-,markerfacecolor,b); hold off shg; pause(dt); end

Lax-Wendroff method

$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)$

這個格式的動機，對於時間我們希望更高的精度：

$left.egin{aligned} u ( x ,t + Delta t ) & = u ( x ,t ) + Delta t partial _ { t } u ( x ,t ) + frac { 1} { 2} Delta t ^ { 2} partial _ { t t } u ( x ,t ) + ldots \ & = u ( x ,t ) - A Delta t partial _ { x } u ( x ,t ) + frac { 1} { 2} A ^ { 2} Delta t ^ { 2} partial _ { x x } u ( x ,t ) + dots end{aligned} ight.$

然後呢，因為時間的差分可以近似時間的偏導:

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

整理一下，就可以得到最上面的格式

（我建議一定要手推一下，記住推導的思想，以後就不用翻書，可以自己現推，而且……不難……）

LW格式的結果很有意思，方法是穩定的，但是：

1）不連續處的解會震蕩。

2）平滑區域的相移。

和cyf跑步去了

預告：

二里講有限體積法

三里講Riemann 問題求解器

四補充具體實現中關於邊界條件，量綱，物理參數等問題

五開始講N-S方程的數值方法

有心情就補充：

磁流體力學MHD方程的建立

保守場的數值方法（重力場、電磁場）

正壓模式和協壓模式

次網格方法和雷諾平均