計算流體力學——從原理到演算法（六）：可壓縮無粘流的計算方法（2）一維歐拉方程Riemann問題的實現

09-04

來自專欄第三新東京市4 人贊了文章

結論與分析

Demo: 一維Euler方程的Riemann問題的初始條件，間斷點位置x=0.3

從上到下：密度，速度，壓力隨著時間變化的情況。（網格點數：100，截止時間0.2s，CFL條件數0.2）

一、一維歐拉方程的黎曼問題：背景

考慮完全氣體的一維歐拉方程：

$left[ egin{array} { l } { ho } \ { ho u } \ { E } end{array} ight] + left[ egin{array} { c c c } { 0} & { 0} & { 1} & { 0} \ { frac { 1} { 2} ( gamma - 3) u ^ { 2} } & { ( 3- gamma ) u } & { gamma - 1} \ { u left[ frac { 1} { 2} ( gamma - 1) u ^ { 2} - H ight] } & { H ( gamma - 1) u ^ { 2} } & { gamma u } end{array} ight] left[ egin{array} { l } { ho } \ { ho u } \ { E } end{array} ight] _ { x } = 0$

的黎曼問題：

$left{ egin{array} { l } { U _ { t } + F ( U ) _ { x } = 0} \ { U ( x ,0) = U ^ { ( 0) } ( x ) } = left{ egin{array} { l l } { U _ { L } ,} & { x < 0} \ { U _ { R } ,} & { x > 0} end{array} ight. end{array} ight.$

相當於實際氣動力學問題中，單個激波shock的理想模型（近似在一個（不管多少維度）的管道中的激波運行問題）

（方程不理解，符號看不懂的參考計算流體力學——從原理到代碼（五）：可壓縮無粘流的計算方法（1）1D 歐拉方程開頭）

這個黎曼方程的相似解 $U ( xi ) = U ( x / t )$ ，見下圖，分別有三個波與三個特徵場相對應聯繫，分別對應右特徵向量 $R^{i}, i= 1,2,3$ 。三個波將 $t>0$ 的上半平面劃分為四個區域，或者說四個常數狀態 $U _ { L } ,U _{* ,L} ,U_{ * ,R} ,U _ { R }$

（雙曲型PDE特徵系統不清楚的可以參考計算流體力學——從原理到代碼（二）：從線性對流方程到黎曼問題）

實驗激波管：兩種不同的氣體狀態，在x = 0 =>特殊的黎曼問題時，被ul = ur = 0的膜隔開。

解通常會給一個衝擊波shock+接觸間斷Contact discontinuity+稀疏波Rarefaction，（對於其他初始條件，可以獲得兩個衝擊或兩個稀疏）但是中間一定是接觸連續——這也為後面Godunov Scheme埋下伏筆

一維Shock Tube理論上解是可以解析求出來的，但是太耗計算量了

理論上，一維黎曼問題的解

感覺講得不是很清楚，如果有什麼問題，請評論區告訴我，我會不斷的更新改進這些。

二、代碼實現部分

1、初始設置部分

input 完全空氣泊松比和理想氣體常數

clcclear% global gamma R;gamma = 1.4;R = 287;

網格初始設置

N = 100; % grid numberCFL = 0.2; % CFL numberx = linspace(0,1,N); % griddx = (x(end)-x(1))/(N-1);% grid spacingx_dis = 0.3; % discontinuity surface locationn_dis = x_dis*N;

設置每個點的初值

$vec{w_0}(x_i) = [ ho, u, p] ^T$

當然，考慮最簡單的Riemann問題，即間斷點 x_dis 左右分別為常數向量，記作 $W_l, W_r$ ? .

W_l = [1,0.75,1];W_r = [0.125,0,0.1];?W = zeros(N,3);for i=1:N if i<=n_dis W(i,:) = W_l; else W(i,:) = W_r; endend

繼續設置初始值

U = W2U(W); % conservation formF = zeros(N-1,3); % flux vector?steps = 0;flag = 1;?current_time = 0;t_max = 0.2;maxSteps = 1e4;

W2U 表示從非守恆量to守恆量，這個函數參考

徐政豪：計算流體力學——從原理到代碼（五）：可壓縮無粘流的計算方法（1）1D 歐拉方程?

zhuanlan.zhihu.com

的最後，包括後面出現的，U2W 和 U2F 。

2、主循環部分

主循環停機條件，達到最大時間或者迭代最大次數。

注意時間步長公式: $dt = ext{CFL} cdot frac{dx} {left| u ight| +sqrt{frac{gamma p}{ ho}}}$ , 其實分母就是激波波速： $lambda_3 = {left| u ight| +sqrt{frac{gamma p}{ ho}}} = |u| + a_{sound}$

note：

CFL條件不熟悉的，參考維基百科吧
這裡我們的CFL設置成了0.2
並不是說CFL越小越好，也不是說小於一就可以，這一部分我其實也不太懂……之後細講

while (flag) steps = steps+1; dt = CFL*dx/max(abs(W(:,2))+sqrt(gamma*W(:,3)./W(:,1))); % timestep if steps > maxSteps % stoping criteria disp(WARNING: maxSteps reached!); break; end if current_time+dt >= t_max % stoping criteria dt = t_max-current_time; flag = 0; end current_time = current_time+dt;

下面的代碼結構在計算流體力學——從原理到代碼（四）：有限體積法與Godunov Scheme 初步裡面有詳細解釋

for i = 1:N-1 % computing F_i+1/2 F(i,:) = godunov_scheme(U_l(i,:), U_r(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? % boundary conditions U(1,:) = U(2,:); U(N,:) = U(N-1,:); W = U2W(U); % visualization subplot(3,1,1); plot(x,W(:,1),bo-,markerfacecolor,r) xlabel(x);ylabel(density ho); title(sprintf(x discontinue=%0.2f,current time=%.3f,CFL=%.2f,N=%d,x_dis,current_time,CFL,N)) subplot(3,1,2); plot(x,W(:,2),bo-,markerfacecolor,b); xlabel(x);ylabel(velosity u); subplot(3,1,3); plot(x,W(:,3),bo-,markerfacecolor,g); xlabel(x);ylabel(pressure p); frame = getframe(h); im = frame2im(frame); [imind,cm] = rgb2ind(im,256); if steps==1 imwrite(imind,cm,1d_euler.gif,gif,Loopcount,Inf,DelayTime,0.0002); else imwrite(imind,cm,1d_euler.gif,gif,WriteMode,append,DelayTime,0.0002); endend

次回預告

我會先把沒寫完（三）Riemann問題Theoretical Analysis補上，有點難度，請大家做好準備

然後是（七）godunov_scheme 的原理（側重Numerical Analysis）與實現