计算流体力学CFD大作业

1 提出问题

[问题描述]

Sod激波管问题是典型的一类Riemann问题。如图所示，一管道左侧为高温高压气体，右侧为低温低压气体，中间用薄膜隔开。t=0 时刻，突然撤去薄膜，试分析其他的运动。

2,u2,p2

Sod模型问题：在一维激波管的左侧初始分布为： 1 1, p1 1, u1 0，右侧分布为： 2 0.125, p2 0.1, u2 0，两种状态之间有一隔膜位于x 0.5处。隔膜突然去掉，试给出在t 0.14时刻Euler方程的准确解，并给出在区间0 x 1这一时刻p, , u的分布图。

2 一维Euler方程组

分析可知，一维激波管流体流动符合一维Euler方程，具体方程如下： 矢量方程：

分量方程：

连续性方程、动量方程和能量方程分别是：

U f

0 t x

(0.1)

计算流体力学CFD大作业

u

0

x t

2

u u p

0

x x t

E u E p 0

x t

(0.2)

u2

其中 E cvT

2

对于完全气体，在量纲为一的形式下，状态方程为：

p T

Ma

2

(0.3)

在量纲为一的定义下，定容热容cv为：

cv

1

2Ma 1 (0.4)

联立(1.2)，(1.3)，(1.4)消去温度T和定容比热cv，得到气体压力公式为：

1 p 1 E u2

2

(0.5)

上式中 为气体常数，对于理想气体 1.4。

3 Euler方程组的离散

3.1 Jacibian矩阵特征值的分裂

Jacibian矩阵A的三个特征值分别是 1 u; 2 u c; 3 u c，依据如下算法将其分裂成正负特征值：

3.2 流通矢量的分裂

这里对流通矢量的分裂选用Steger-Warming分裂法，分裂后的流通矢量为

＋＋＋ 2 1 1 2 3

＋＋＋＋

f 2 1 1u 2 u c 3 u c

2 ＋＋

1 ＋u2 2u c2 3u c2 w＋

1 22

k

k k2 2

2

(0.6)

(0.7)

计算流体力学CFD大作业

－－－ 2 1 1 2 3

－－－－

f 2 1 1u 2 u c 3 u c

2 －－

22 1 －u2 2u c 3u c w－

1 22

2

3 c 23

w

2 1 (0.8)

其中：

c为量纲为一的声速：

c2

T

2

Ma

p

(0.9)

联立(1.3)，(1.9)式，消去来流马赫数得：c 3.3 一阶迎风显示格式离散Euler方程组 得到

算法如下：

n

Un+1j=Uj

Uin 1 Uin x fi x fi

0 t x

t

fj fj f f 1j 1j x

(0.10)

00

① 已知初始时刻t=0的速度、压力及密度分布u0j,Pj, j，则可得到特征值分裂

值 k0 ，从而求出流通矢量fj0 ；

② 应用一阶迎风显示格式可以计算出t 1 t时刻的组合变量U1j，从而得到

t 1 t时刻的速度、压力及密度分布u1j,Pj1, 1j；

1 ③ 利用t 1 t时刻的速度、压力及密度分布u1j,Pj1, 1j可得特征值分裂值 k，

从而求出流通矢量fj1 ；

22④ 按照步骤2的方法即可得到t 2 t时刻的速度、压力及密度分布u2j,Pj, j；

⑤ 循环以上过程即可得到t n 1 t时刻的速度、压力及密度分布

n+1un+1, n+1j,Pjj。

计算流体力学CFD大作业

4 计算结果分析

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

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

图4-1是0.14时刻的密度分布图，观察可知，在密度波的传播过程中，间断面上会出现了两次“沉降”，说明密度在沉降位置发生了剧烈变化。图4-2是0.14时刻的压力分布图，在压力波的传播过程中，在间断面上出现了一个“压力沉降”现象，说明压力在沉降位置突降。图4-3是0.14时刻的速度分布图，在间断面处产生一个向两边运动的速度，并且只有在隔膜附近才有气体流动，其他地方静止。

图4-1激波管内密度分布图(0.14s)

计算流体力学CFD大作业

图4-2激波管内压力分布图(0.14s)

图4-3激波管内速度分布图(0.14s)

源程序代码：

Ⅰ 分量和矩阵结合编写的源程序： 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_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:Nt

epso=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:3

lamta_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));

计算流体力学CFD大作业

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 ...…… 此处隐藏：2383字，全部文档内容请下载后查看。喜欢就下载吧 ……