计算流体力学大作业

1、1提出问题问题描述Sod激波管问题是典型的一类 Riemann问题。如图所示，一管道左侧 为高温高压气体，右侧为低温低压气体，中间用薄膜隔开。 t=0时刻，突然撤去 薄膜，试分析其他的运动。Sod模型问题：在一维激波管的左侧初始分布为： 咯=1, pi =1, ui =。，右侧分 布为：p2 =0.125, p2=0.1, u2 =0，两种状态之间有一隔膜位丁 X = O.5处。隔膜 突然去掉，试给出在t=0.14时刻Euler方程的准确解，并给出在区间°主x<1这 一时刻P，P，u的分布图。2 一维Euler方程组分析可知，一维激波管流体流动符合一维 Euler方程，具体方程

2、如下： 欠量方程：:U: f二 一=0(0.1):t :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对丁完全气体，在量纲为一的形式下，状态方程为:p =仃Ma2：(0.3)在量纲为一的定义下，定容热容cv为:Cv =9Ma 2-1(0.4)联立(1.2), (1.3), (1.4)消去温度T和定容比热Cv,得到气体压力公式为:(0.5)上式中丫为气体常数，对丁理想气体丫=1.4。3 Euler方程组

3、的离散3.1 Jacibian矩阵特征值的分裂Jacibian矩阵A的三个特征值分别是X, = u;人2 = u+c; = u-c ,依据如下算法将其分裂成正负特征值:(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 二联立(1.3), (1.9)式，消去来流马赫数得：c = , P、3.3 一阶迎风显示格式离散Euler方程组(0.10)Uin 1 -Uin'" '"一八 =0Lt匚X得到Un+1=U； -三_ f-fj顷-丫算法如下： 已知初始时刻t=0的速度、压力及密度分布u0,Pj°,P0,则可得到特征值分裂值沱土，从而求出流通欠量f土;应用一阶迎风显示格式可以计算出t =1约 时刻的组合变

5、量U 1 ,从而得到t =1割时刻的速度、压力及密度分布u1j,P1,P：; 利用t =1或时刻的速度、压力及密度分布u1, Pj1, P；可得特征值分裂值从而求出流通欠量f ； £ 按照步骤2的方法即可得到t=2约时刻的速度、压力及密度分布u2, Pj2f2； 循环以上过程即可得到t=(n+1陪t时刻的速度、压力及密度分布n+1n+1- n+1uj , Pj,Pj。4计算结果分析实际编程中，空间步长取0.001,空间网格数为1001,时间步长取0.00001, 计算到终点时刻0.14s耗费机时137s,计算时间还是可以接受的。分析图4-14-3,可以观察到在隔膜附近流动参数变化剧烈

6、，与初始条件相比，可以看出激波的影响范围有限，始终在xw 0.3,0.75区间内变化。图4-1是0.14时刻的密度分布图，观察可知，在密度波的传播过程中，问断 面上会出现了两次“沉降”，说明密度在沉降位置发生了剧烈变化。图4-2是0.14时刻的压力分布图，在压力波的传播过程中，在问断面上出现了一个“压力沉降” 现象，说明压力在沉降位置突降。图4-3是0.14时刻的速度分布图，在问断面处 产生一个向两边运动的速度，并且只有在隔膜附近才有气体流动，其他地方静止。图4-1激波管内密度分布图(0.14s)11 i i i i 一 i i i o.i n? n3 cu ns ob o? ns g i图4

7、-2激波管内压力分布图(0.14s)in0/06S 0$瑕0.40.30.20.1ai 0.20.306 0E 0? DO D9X图4-3激波管内速度分布图(0.14s)源程序代码：I分量和矩阵结合编写的源程序：function sobtubing_SW() tic;close all ee=1e-8;%网格 %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);%翳数 %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,&