一维黎曼问题数值解与计算程序.docx

一维黎曼问题数值解与计算程序.docx

《一维黎曼问题数值解与计算程序.docx》由会员分享，可在线阅读，更多相关《一维黎曼问题数值解与计算程序.docx（15页珍藏版）》请在冰豆网上搜索。

一维黎曼问题数值解与计算程序

一维Riemann问题数值解与计算程序

一维Rieman"可题，即激波管问题，是一个典型的一维可压缩无黏气体动力

学问题，并有解析解。

对它采用二阶精度MacCormac两步差分格式进行数值求解。

同时，为了初学者入门和练习方便,这里给出了用C语言和Fortran77编写的计算一维Riemanri可题的计算程序，供大家学习参考。

A-1利用MacCormac两步差分格式求解一维Riemann问题

1.一维Riemanr问题

一维Riemand可题实际上就是激波管问题。

激波管是一根两端封闭、内部充

满气体的直管，如图A.1所示。

在直管中由一薄膜将激波管隔开，在薄膜两侧充

有均匀理想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。

当t乞0时，气体处于静止状态。

当t=0时，薄膜瞬时突然破裂，气体从高压端冲向低压端，同时在管内形成激波、稀疏波和接触间断等复杂波系。

2.基本方程组、初始条件和边界条件

//////////////////////////////////////

图A.1激波管问题示意图

设气体是理想气体。

一维Riemanri可题在数学上可以用一维可压缩无黏气体

Eulei方程组来描述。

在直角坐标系下量纲为一的一维Euler方程组为:

■^+f=0,—1兰x兰1（A.1）

.:

t:

x

rP、

「Pu、

其中

u=

Pu

f=

Pu2+p

（A.2）

l（E+p）u』

这里「、u、p、E分别是流体的密度、速度、压力和单位体积总能。

理想气体

状态方程：

p=-1◎-1E-Ju2v2（A.3）

初始条件：

J=1,5=0,山=1；=0.125,u2=0,p2=0.1。

边界条件：

X--1和X=1处为自由输出条件，UohU—UnhUn/。

3.二阶精度MacCormac差分格式

n-2

MacCormac两步差分格式:

2nn-n

Uj-Uj-rfj-fjj

其中r二卫。

计算实践表明，MacCormac两步差分格式不能抑制激波附近非物

Ax

理振荡。

因此在计算激波时，必须采用人工黏性滤波方法:

（A.5）

需要引入一个与

-nn1nQnn

Ui,j=Ui,j-U「1,j-2嘔-UiJ,j

为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零,

密度（或者压力）相关的开关函数:

由式（A.6）可以看出，在光滑区域，密度变化很缓，因此二值也很小；而在激波附近密度变化很陡，V值就很大。

带有开关函数的前置人工黏性滤波方法为:

（A.7）

_1

Uinj=unj+罗日（ut-2uinj+ulj

其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到:

（A.8）

由于声速不会超过3,所以取|a|=3，在本计算中取=0.25。

4.计算结果分析

计算分别采用标准的C语言和Fortran77语言编写程序。

计算中网格数取1000，计算总时间为T=0.4。

计算得到在T=0.4时刻的密度、速度和压力分布

如图A.2（C语言计算结果）和图A.3（Fortran77语言计算结果）所示。

采用两

种不同语言编写程序所得到的计算结果完全吻合。

从图A.2和图A.3中可以发现，MacCormac两步差分格式能很好地捕捉激波,计算得到的激波面很陡、很窄，计算激波精度是很高的。

采用带开关函数的前置

人工滤波法能消除激波附近的非物理振荡，计算效果很好

从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的，而密度和速度是相等的。

这个计算结果正确地反映了一维Riemann问题的物理特性，并被激波管实验所验证。

A-2一维Riemanr问题数值计算源程序

1.C语言源程序

//MacCormackID.cpp:

定义控制台应用程序的入口点。

//

/*

*利用MacCormac差分格式求解一维激波管问题（C语言版本）

*

*/

〃#include"stdafx.h"

#include

#include

#include#defineGAMA1.4〃气体常数

#definePI3.141592654

#defineL2.0〃计算区域

#defineTT0.4//总时间

#defineSf0.8//时间步长因子

#defineJ1000//网格数

//全局变量

doubleU[J+2][3],Uf[J+2][3],Ef[J+2][3];

/*

计算时间步长

入口：

U,当前物理量，dx，网格宽度；

返回:

时间步长。

*/

doubleCFL（doubleU[J+2][3],doubledx）

{

inti;

doublemaxvel,p,u,vel;

maxvel=1e-100;

for（i=1;i<=J;i++）

{

u=U[i][1]/U[i][0];

p=（GAMA-1）*（U[i][2]-0.5*U[i][0]*u*u）;

vel=sqrt（GAMA*p/U[i][0]）+fabs（u）;if（vel>maxvel）maxvel=vel;

}

returnSf*dx/maxvel;

}

/*

初始化

入口：

无；

出口：

U，已经给定的初始值，

dx,网格宽度。

*/

voidInit（doubleU[J+2][3],double&dx）

{

inti;

doublerou1=1.0,u1=0.0,p1=1.0;//初始条件

doublerou2=0.125,u2=0.0,p2=0.1;

dx=L/J;

for（i=0;i<=J/2;i++）

{

U[i][0]=rou1;

U[i][1]=rou1*u1;U[i][2]=p1/（GAMA-1）+rou1*u1*u1/2;

}

for（i=J/2+1;i<=J+1;i++）

{

U[i][0]=rou2;

U[i][1]=rou2*u2;

U[i][2]=p2/（GAMA-1）+rou2*u2*u2/2;

}

}

/*

边界条件

入口：

dx，网格宽度；出口:

U，已经给定的边界。

*/voidbound（doubleU[J+2][3],doubledx）{

intk;

//左边界for（k=0;k<3;k++）U[0][k]=U[1][k];

//右边界

for（k=0;k<3;k++）U[J+1][k]=U[J][k];

}

/*

根据U计算E入口：

U，当前U矢量；出口：

E，计算得到的E矢量，

U、E的定义见Euler方程组。

*/voidU2E（doubleU[3],doubleE[3]）{

doubleu,p;

u=U[1]/U[0];p=（GAMA-1）*（U[2]-0.5*U[1]*U[1]/U[0]）;

E[0]=U[1];

E[1]=U[0]*u*u+p;E[2]=（U[2]+p）*u;

}

/*

一维MacCormac差分格式求解器

入口：

u,上一时刻的u矢量，Uf、Ef,临时变量，

dx，网格宽度，dt,时间步长；

出口：

u,计算得到的当前时刻u矢量。

*/

Ef[J+2][3],double

voidMacCormack_1DSolver（doubleu[J+2][3],doubleuf[J+2][3],doubledx,doubledt）

{

inti,k;

doubler,nu,q;

r=dt/dx;

nu=0.25;

for（i=1;i<=J;i++）

{

q=fabs（fabs（u[i+1][0]-u[i][0]）-fabs（u[i][0]-u[i-1][0]））/（fabs（u[i+1][0]-u[i][0]）+fabs（u[i][0]-u[i-1][0]）+1e-100）;//开关函数

for（k=0;k<3;k++）

Ef[i][k]=u[i][k]+0.5*nu*q*（u[i+1][k]-2*u[i][k]+u[i-1][k]）;//人工黏性项}

for（k=0;k<3;k++）for（i=1;i<=J;i++）u[i][k]=Ef[i][k];

for（i=0;i<=J+1;i++）u2E（u[i],Ef[i]）;

for（i=0;i<=J;i++）

for（k=0;k<3;k++）uf[i][k]=u[i][k]-r*（Ef[i+1][k]-Ef[i][k]）;//u（n+1/2）（i+1/2）

for（i=0;i<=J;i++）u2E（uf[i],Ef[i]）;//E（n+1/2）（i+1/2）

for（i=1;i<=J;i++）

for（k=0;k<3;k++）u[i][k]=0.5*（u[i][k]+uf[i][k]）-0.5*r*（Ef[i][k]-Ef[i-1][k]）;//u（n+1）（i）

}

/*

输出结果,用Origin数据格式画图

入口：

u,当前时刻u矢量,dx,网格宽度；

出口：

无。

*/

voidOutput（doubleU[J+2][3],doubledx）

{

inti;

FILE*fp;

doublerou,u,p;

fp=fopen（"result.txt","w"）;

for（i=0;i<=J+1;i++）

{

rou=U[i][0];

u=U[i][1]/rou;

p=（GAMA-1）*（U[i][2]-0.5*U[i][0]*u*u）;

fprintf（fp,"%20f%20.10e%20.10e%20.10e%20.10e\n",i*dx,rou,u,p,U[i][2]）;}fclose（fp）;

}

/*

主函数入口:

无；出口:

无。

*/

voidmain（）

{

doubleT,dx,dt;

Init（U,dx）;

T=0;while（T

{dt=CFL（U,dx）;

T+=dt;

printf（"T=%10gdt=%10g\n",T,dt）;

MacCormack_1DSolver（U,Uf,Ef,dx,dt）;

bound（U,dx）;

}

Output（U,dx）;

}

2.Fortran7语言源程序

!

MacCormack1D.for

!

利用MacCormac差分格式求解一维激波管问题（Fortran77语言版本）*/programMacCormack1Dimplicitdoubleprecision（a-h,o-z）parameter（M=1000）common/G_def/GAMA,PI,J,JJ,dL,TT,Sf

dimensionU（0:

M+1,0:

2）,Uf（0:

M+1,0:

2）dimensionEf（0:

M+1,0:

2）

!

气体常数

GAMA=1.4PI=3.1415926

!

网格数

J=M

!

计算区域

dL=2.0

!

总时间

TT=0.4

!

时间步长因子

Sf=0.8

callInit（U,dx）

T=0

1dt=CFL（U,dx）

T=T+dtwrite（*,*）'T=',T,'dt=',dtcallMacCormack_1D_Solver（U,Uf,Ef,dx,dt）callbound（U,dx）

if（T.lt.TT）goto1

callOutput（U,dx）end

!

!

计算时间步长

!

入口:

U，当前物理量，dx,网格宽度；

!

返回:

时间步长。

!

doubleprecisionfunctionCFL（U,dx）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0:

J+1,0:

2）

dmaxvel=1e-10

do10i=1,J

uu=U（i,1）/U（i,0）p=（GAMA-1）*（U（i,2）-0.5*U（i,0）*uu*uu）vel=dsqrt（GAMA*p/U（i,0））+dabs（uu）if（vel.gt.dmaxvel）dmaxvel=vel

10continue

CFL=Sf*dx/dmaxvel

end

!

!

初始化

!

入口:

无；

!

出口：

U,已经给定的初始值，dx,网格宽度。

!

subroutineInit（U,dx）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0：

J+1,0：

2）

!

初始条件

rou1=1.0

u1=0

v1=0

p1=1.0

rou2=0.125

u2=0

v2=0

p2=0.1

dx=dL/J

do20i=0,J/2

U（i,0）=rou1U（i,1）=rou1*u1

U（i,2）=p1/（GAMA-1）+0.5*rou1*u1*u1

20continue

do21i=J/2+1,J+1

U（i,0）=rou2

U（i,1）=rou2*u2

U（i,2）=p2/（GAMA-1）+0.5*rou2*u2*u221continue

end

!

!

边界条件

!

入口：

dx，网格宽度；

!

出口:

U，已经给定边界。

!

subroutinebound（U,dx）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0：

J+1,0：

2）

!

左边界

do30k=0,2

U（0,k）=U（1,k）

30continue

!

右边界

do31k=0,2

U（J+1,k）=U（J,k）

31continue

end

!

!

根据U计算E

!

入口：

U，当前U矢量；

!

出口：

E，计算得到的E矢量，

!

U、E定义见Euler方程组。

!

subroutineU2E（U,E,is,in）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0：

J+1,0：

2）,E（0：

J+1,0：

2）

do40i=is,in

uu=U（i,1）/U（i,0）p=（GAMA-1）*（U（i,2）

$-0.5*U（i,1）*U（i,1）/U（i,0））

E（i,0）=U（i,1）

E（i,1）=U（i,0）*uu*uu+p

E（i,2）=（U（i,2）+p）*uu

40continue

end

!

!

一维MacCormac差分格式求解器

!

入口:

U，上一时刻U矢量，

!

Uf、Ef，临时变量，

!

dx，网格宽度，dt,,时间步长；

!

出口:

U，计算得到得当前时刻U矢量。

!

subroutineMacCormack_1D_Solver（U,Uf,Ef,dx,dt）implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0:

J+1,0:

2）,Uf（0:

J+1,0:

2）dimensionEf（0:

J+1,0:

2）

r=dt/dx

dnu=0.25

do60i=1,J

do60k=0,2

!

开关函数

q=dabs（dabs（U（i+1,0）-U（i,0））-dabs（U（i,0）-U（i-1,0）））

$/（dabs（U（i+1,0）-U（i,0））+dabs（U（i,0）-U（i-1,0））+1e-10）!

人工黏性项

Ef（i,k）=U（i,k）+0.5*dnu*q*（U（i+1,k）-2*U（i,k）+U（i-1,k））

60continue

do61k=0,2

do61i=1,J

U（i,k）=Ef（i,k）

61continue

callU2E（U,Ef,0,J+1）

do63i=0,J

do63k=0,2

!

U（n+1/2）（i+1/2）

Uf（i,k）=U（i,k）-r*（Ef（i+1,k）-Ef（i,k））

63continue

!

E（n+1/2）（i+1/2）

callU2E（Uf,Ef,0,J）

do64i=1,J

do64k=0,2

!

U（n+1）（i）

U（i,k）=0.5*（U（i,k）+Uf（i,k））-0.5*r*（Ef（i,k）-Ef（i-1,k））

64continue

end

!

!

输出结果,用Origin数据格式画图

!

入口:

U，当前时刻U矢量，

!

dx，网格宽度；

!

出口:

无。

!

subroutineOutput（U,dx）

implicitdoubleprecision（a-h,o-z）common/G_def/GAMA,PI,J,JJ,dL,TT,SfdimensionU（0:

J+1,0:

2）

open（1,file='result.txt',status='unknown'）

do80i=0,J+1

rou=U（i,0）

uu=U（i,1）/roup=（GAMA-1）*（U（i,2）-0.5*U（i,0）*uu*uu）write（1,81）i*dx,rou,uu,p,U（i,2）

80continue

close

（1）

81format（D20.10,D20.10,D20.10,D20.10,D20.10）end