計算流體力學(xué)大作業(yè)

1、1提出問題問題描述Sod激波管問題是典型的一類 Riemann問題。如圖所示，一管道左側(cè) 為高溫高壓氣體，右側(cè)為低溫低壓氣體，中間用薄膜隔開。 t=0時刻，突然撤去 薄膜，試分析其他的運動。Sod模型問題：在一維激波管的左側(cè)初始分布為： 咯=1, pi =1, ui =。，右側(cè)分 布為：p2 =0.125, p2=0.1, u2 =0，兩種狀態(tài)之間有一隔膜位丁 X = O.5處。隔膜 突然去掉，試給出在t=0.14時刻Euler方程的準確解，并給出在區(qū)間°主x<1這 一時刻P，P，u的分布圖。2 一維Euler方程組分析可知，一維激波管流體流動符合一維 Euler方程，具體方程

2、如下： 欠量方程：:U: f二 一=0(0.1):t :X分量方程：連續(xù)性方程、動量方程和能量方程分別是:Lap a(pu)0十=0 ctGX一 r 2"u " u 里二0ctexex虎雄(E + p)o；-t: x(0.2)/2、其中 E = P cvT +L <2 J對丁完全氣體，在量綱為一的形式下，狀態(tài)方程為:p =仃Ma2：(0.3)在量綱為一的定義下，定容熱容cv為:Cv =9Ma 2-1(0.4)聯(lián)立(1.2), (1.3), (1.4)消去溫度T和定容比熱Cv,得到氣體壓力公式為:(0.5)上式中丫為氣體常數(shù)，對丁理想氣體丫=1.4。3 Euler方程組

3、的離散3.1 Jacibian矩陣特征值的分裂Jacibian矩陣A的三個特征值分別是X, = u;人2 = u+c; = u-c ,依據(jù)如下算法將其分裂成正負特征值:(0.6)3.2流通欠量的分裂這里對流通欠量的分裂選用Steger-Warming分裂法,分裂后的流通欠量為(0.7)2(了-1)妃 +# + 萬|-| +32(丫 一1)& u + '；(u_c)+7；(u+c)fi +1 + + 22 A a2+('T )，- u + (uc)+ (u+c)+w 22J2( ? 1 )兀 + A2 + A_3(0.8)2 (了 一 1)孫 u + Z2 (u c )

4、+ Z3 (u + c )n 2/32-(u c) +-( u +c)十w其中:c為量綱為一的聲速:3 ii，2一 ,與一 cw =2-1-第7頁(0.9)T2Ma 二聯(lián)立(1.3), (1.9)式，消去來流馬赫數(shù)得：c = , P、3.3 一階迎風(fēng)顯示格式離散Euler方程組(0.10)Uin 1 -Uin'" '"一八 =0Lt匚X得到Un+1=U； -三_ f-fj頃-丫算法如下： 已知初始時刻t=0的速度、壓力及密度分布u0,Pj°,P0,則可得到特征值分裂值沱土，從而求出流通欠量f土;應(yīng)用一階迎風(fēng)顯示格式可以計算出t =1約 時刻的組合變

5、量U 1 ,從而得到t =1割時刻的速度、壓力及密度分布u1j,P1,P：; 利用t =1或時刻的速度、壓力及密度分布u1, Pj1, P；可得特征值分裂值從而求出流通欠量f ； £ 按照步驟2的方法即可得到t=2約時刻的速度、壓力及密度分布u2, Pj2f2； 循環(huán)以上過程即可得到t=(n+1陪t時刻的速度、壓力及密度分布n+1n+1- n+1uj , Pj,Pj。4計算結(jié)果分析實際編程中，空間步長取0.001,空間網(wǎng)格數(shù)為1001,時間步長取0.00001, 計算到終點時刻0.14s耗費機時137s,計算時間還是可以接受的。分析圖4-14-3,可以觀察到在隔膜附近流動參數(shù)變化劇烈

6、，與初始條件相比，可以看出激波的影響范圍有限，始終在xw 0.3,0.75區(qū)間內(nèi)變化。圖4-1是0.14時刻的密度分布圖，觀察可知，在密度波的傳播過程中，問斷 面上會出現(xiàn)了兩次“沉降”，說明密度在沉降位置發(fā)生了劇烈變化。圖4-2是0.14時刻的壓力分布圖，在壓力波的傳播過程中，在問斷面上出現(xiàn)了一個“壓力沉降” 現(xiàn)象，說明壓力在沉降位置突降。圖4-3是0.14時刻的速度分布圖，在問斷面處 產(chǎn)生一個向兩邊運動的速度，并且只有在隔膜附近才有氣體流動，其他地方靜止。圖4-1激波管內(nèi)密度分布圖(0.14s)11 i i i i 一 i i i o.i n? n3 cu ns ob o? ns g i圖4

7、-2激波管內(nèi)壓力分布圖(0.14s)in0/06S 0$瑕0.40.30.20.1ai 0.20.306 0E 0? DO D9X圖4-3激波管內(nèi)速度分布圖(0.14s)源程序代碼：I分量和矩陣結(jié)合編寫的源程序：function sobtubing_SW() tic;close all ee=1e-8;%網(wǎng)格 %delta_t=0.00001;Nt=round(0.14/delta_t);delta_x=0.001;N_left=round(0.5/delta_x+1);N_right=round(0.5/delta_x+1) N=N_left+N_right-1;%1%P=ones(1,N_

8、left-1) 0.1*ones(1,N_right);Den=ones(1,N_left-1) 0.125*ones(1,N_right);u=zeros(1,N);gama=1.4;Den_u=Den.*u;E=P./(gama-1) +(0.5*u.A2).*Den;%值分裂 %for j=1:Ntepso=1e-8*ones(1,N);gama=1.4;C=sqrt(gama*P./Den);lamta=ones(3,N);lamta_p=ones(3,N);lamta_n=ones(3,N);lamta(1,:)=u; lamta(2,:)=u-C; lamta(3,:)=u+C;

9、for i=1:3lamta_p(i,:)=0.5*(lamta(i,:)+sqrt(lamta(i,:).A2+epso.A2);lamta_n(i,:)=0.5*(lamta(i,:)-sqrt(lamta(i,:).A2+epso.A2);end % % gama=1.4;C=sqrt(gama*P./Den);Ftran_p=ones(3,N);f1_Pos=0.5/gama*Den.*(2*(gama-1)*lamta_p(1,:)+lamta_p(2,:)+lamta_p(3,：)；f2_Pos=0.5/gama*Den.*(2*(gama-1)*lamta_p(1,:).*u+l

10、amta_p(2,:).*(u-C)+ lamta_p(3,:).*(u+C);f3_Pos=0.5/gama*Den.*(gama-1)*lamta_p(1,:).*u.A2+0.5*lamta_p(2,:).*(u-C).A2 .+0.5*lamta_p(3,:).*(u+C).A2+(0.5*(3-gama)/(gama-1)*(lamta_p(2,:)+lamta_p(3,:).*C.A2);% % gama=1.4;C=sqrt(gama*P./Den);f1_Neg=0.5/gama*Den.*(2*(gama-1)*lamta_n(1,:)+lamta_n(2,:)+lamta_

11、n(3, :);f2_Neg=0.5/gama*Den.*(2*(gama-1)*lamta_n(1,:).*u+lamta_n(2,:).*(u-C)+ lamta_n(3,:).*(u+C);f3_Neg=0.5/gama*Den.*(gama-1)*lamta_n(1,:).*u.A2+0.5*lamta_n(2,:).*(u-C).A2 .+0.5*lamta_n(3,:).*(u+C).A2+(0.5*(3-gama)/(gama-1)*(lamta_n(2,:)+lamta_n(3,:).*C.A2);%翳數(shù) %for i=2:N-1% 密度計算temp1(1,i) = (f1_P

12、os(1,i) - f1_Pos(1,i-1) / delta_x;temp2(1,i) = (f1_Neg(1,i+1) - f1_Neg(1,i) / delta_x;Den(1,i) = Den(1,i) - delta_t*(temp1(1,i) + temp2(1,i);% 密度、速度乘積計算temp1(1,i) = (f2_Pos(1,i) - f2_Pos(1,i-1) / delta_x;temp2(1,i) = (f2_Neg(1,i+1) - f2_Neg(1,i) /delta_x;Den_u(1,i) = Den_u(1,i) - delta_t*(temp2(1,i)

13、 + temp1(1,i);%速度計算u(1,i) = Den_u(1,i) / Den(1,i);%能量計算弗，Ui, Pitemp1(1,i) = (f3_Pos(1,i) - f3_Pos(1,i-1) /delta_x;temp2(1,i) = (f3_Neg(1,i+1) - f3_Neg(1,i) /delta_x;E(1,i) = E(1,i) - delta_t*(temp2(1,i) + temp1(1,i);% 壓強計算P(1,i) = (gama - 1)*(E(1,i) - 0.5*Den(1,i)*u(1,i)A2);endend%吉果顯示x=0:0.001:1;axis(0 1 0 1);plot(x,u,'LineWidth',2);xlabel('fontsize14x');ylabel('fontsize14 速度V');figure (2) ;plot(x,P,&