计算流体力学大作业

1 提出问题[问题描述]Sod 激波管问题是典型的一类Riemann 问题。

如图所示，一管道左侧为高温高压气体，右侧为低温低压气体，中间用薄膜隔开。

t=0 时刻，突然撤去薄膜，试分析其他的运动。

Sod 模型问题：在一维激波管的左侧初始分布为：0 ,1 ,1111===u p ρ，右侧分布为：0 ,1.0 ,125.0222===u p ρ，两种状态之间有一隔膜位于5.0=x 处。

隔膜突然去掉，试给出在14.0=t 时刻Euler 方程的准确解，并给出在区间10≤≤x 这一时刻u p , ,ρ的分布图。

2 一维Euler 方程组分析可知，一维激波管流体流动符合一维Euler 方程，具体方程如下： 矢量方程：0U ft x∂∂+=∂∂ (0.1)分量方程：连续性方程、动量方程和能量方程分别是：222,,p u ρ()()()()2000u tx u u pt x x u E p E tx ρρρρ∂⎧∂+=⎪∂∂⎪⎪∂∂∂⎪++=⎨∂∂∂⎪⎪∂+⎡⎤∂⎣⎦+=⎪∂∂⎪⎩ (0.2)其中 22v u E c T ρ⎛⎫=+ ⎪⎝⎭对于完全气体，在量纲为一的形式下，状态方程为：()2p T Ma ργ∞=(0.3)在量纲为一的定义下，定容热容v c 为：()211v c Ma γγ∞=- (0.4)联立(1.2)，(1.3)，(1.4)消去温度T 和定容比热v c ，得到气体压力公式为：()2112p E u γρ⎛⎫=-- ⎪⎝⎭(0.5)上式中γ为气体常数，对于理想气体4.1=γ。

3 Euler 方程组的离散3.1 Jacibian 矩阵特征值的分裂Jacibian 矩阵A 的三个特征值分别是123;;u u c u c λλλ==+=-，依据如下算法将其分裂成正负特征值：()12222k k k λλελ±±+=(0.6)3.2 流通矢量的分裂这里对流通矢量的分裂选用Steger-Warming 分裂法，分裂后的流通矢量为()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭＋＋＋＋＋＋＋＋＋＋＋(0.7)()()()()()()()12312322232121212122f u u c u c u u c u c w γλλλργλλλγλλγλ⎛⎫⎪-++ ⎪=-+-++ ⎪ ⎪ ⎪-+-+++ ⎪⎝⎭－－－－－－－－－－－(0.8)其中：()()()223321c w γλλγ±±±-+=- c 为量纲为一的声速：22Tc Ma ∞=(0.9)联立(1.3)，(1.9)式，消去来流马赫数得：ργp c =3.3 一阶迎风显示格式离散Euler 方程组 10n n i i x i x i U U f f t xδδ+-++--++=∆∆ (0.10)得到()()n+1nj j 11U =U j j j j t f f f f x++---+∆⎡⎤--+-⎣⎦∆ 算法如下：① 已知初始时刻t=0的速度、压力及密度分布000,,j j j u P ρ，则可得到特征值分裂值0k λ±，从而求出流通矢量0j f ±；② 应用一阶迎风显示格式可以计算出1t t =∆时刻的组合变量1j U ，从而得到1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ；③ 利用1t t =∆时刻的速度、压力及密度分布111,,j j j u P ρ可得特征值分裂值1k λ±，从而求出流通矢量1j f ±；④ 按照步骤2的方法即可得到2t t =∆时刻的速度、压力及密度分布222,,j j j u P ρ；⑤ 循环以上过程即可得到()1t n t =+∆时刻的速度、压力及密度分布n+1n+1n+1,,j j j u P ρ。

4 计算结果分析实际编程中，空间步长取0.001，空间网格数为1001，时间步长取0.00001，计算到终点时刻0.14s耗费机时137s，计算时间还是可以接受的。

分析图4-1~4-3，可以观察到在隔膜附近流动参数变化剧烈，与初始条件相比，可以看出激波的影响范围有限，始终在[0.3,0.75]x 区间内变化。

图4-1是0.14时刻的密度分布图，观察可知，在密度波的传播过程中，间断面上会出现了两次“沉降”，说明密度在沉降位置发生了剧烈变化。

图4-2是0.14时刻的压力分布图，在压力波的传播过程中，在间断面上出现了一个“压力沉降”现象，说明压力在沉降位置突降。

图4-3是0.14时刻的速度分布图，在间断面处产生一个向两边运动的速度，并且只有在隔膜附近才有气体流动，其他地方静止。

图4-1激波管内密度分布图(0.14s)图4-2激波管内压力分布图(0.14s)图4-3激波管内速度分布图(0.14s)源程序代码：Ⅰ分量和矩阵结合编写的源程序：function sobtubing_SW()tic;close allee=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_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.^2).*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;for i=1:3lamta_p(i,:)=0.5*(lamta(i,:)+sqrt(lamta(i,:).^2+epso.^2)); lamta_n(i,:)=0.5*(lamta(i,:)-sqrt(lamta(i,:).^2+epso.^2));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+lamta_p(2,:).*(u-C)+ lamta_p(3,:).*(u+C));f3_Pos=0.5/gama*Den.*((gama-1)*lamta_p(1,:).*u.^2+0.5*lamta_p(2,:).*( u-C).^2 ...+0.5*lamta_p(3,:).*(u+C).^2+(0.5*(3-gama)/(gama-1)*(lamta_p(2,:)+lamt a_p(3,:)).*C.^2));%%%%%%%%%%%%%%计算负通量%%%%%%%%%%%%%%%gama=1.4;C=sqrt(gama*P./Den);f1_Neg=0.5/gama*Den.*(2*(gama-1)*lamta_n(1,:)+lamta_n(2,:)+lamta_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.^2+0.5*lamta_n(2,:).*( u-C).^2 ...+0.5*lamta_n(3,:).*(u+C).^2+(0.5*(3-gama)/(gama-1)*(lamta_n(2,:)+lamt a_n(3,:)).*C.^2));%%%%%%%%%%%%%%计算流动参数%%%%%%%%%%%%%%for i=2:N-1%密度计算temp1(1,i) = (f1_Pos(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) + temp1(1,i));% 速度计算u(1,i) = Den_u(1,i) / Den(1,i);% 能量计算temp1(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)^2);endend%结果显示x=0:0.001:1;axis([0 1 0 1]);plot(x,u,'LineWidth',2);xlabel('\fontsize{14}x');ylabel('\fontsize{14}速度V');figure(2);plot(x,P,'LineWidth',2);xlabel('\fontsize{14}x');ylabel('\f ontsize{14}压力P')figure(3);plot(x,Den,'LineWidth',2);xlabel('\fontsize{14}x');ylabel(' \fontsize{14}密度ρ')Calculate_time=toc。