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

一维Riemann 问题数值解与计算程序
一维Riemann 问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有 解析解。

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

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

A-1利用MacCormack 两步差分格式求解一维Riemann 问题
1.一维Riemann 问题
一维Riemann 问题实际上就是激波管问题。

激波管是一根两端封闭、内部充满气体的直管，如图A.1所示。

在直管中由一薄膜将激波管隔开，在薄膜两侧充有均匀理想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。

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

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

2.基本方程组、初始条件和边界条件
设气体是理想气体。

一维Riemann 问题在数学上可以用一维可压缩无黏气体
Euler 方程组来描述。

在直角坐标系下量纲为一的一维Euler 方程组为：
,11x t x
∂∂+=-≤≤∂∂0u f
(A.1) 其中 2,()u u u p E E p u ρρρρ⎛⎫⎛⎫ ⎪ ⎪
==+ ⎪ ⎪ ⎪ ⎪+⎝⎭⎝⎭
u f (A.2)
这里ρ、u 、p 、E 分别是流体的密度、速度、压力和单位体积总能。

理想气体状态方程：
()()()221112p e E u v γργρ⎡⎤
=-=--+⎢⎥⎣⎦
（A.3）
初始条件：1111,0,1u p ρ===；2220.125,0,0.1u p ρ===。

图A.1 激波管问题示意图
边界条件：1x =-和1x =处为自由输出条件，01u u =,1N N u u -=。

3.二阶精度MacCormack 差分格式
MacCormack 两步差分格式：
()
12
1111
1222111
22n n n n
j
j j j n n n n n j
j j j j r r +
-+++++=--⎛⎫⎛⎫=+-- ⎪ ⎪⎝⎭⎝⎭
u u f f u
u u f f (A.4)
其中t
r x
∆=
∆。

计算实践表明，MacCormack 两步差分格式不能抑制激波附近非物 理振荡。

因此在计算激波时，必须采用人工黏性滤波方法：
(),,1,,1,122
n n
n n n i j i j
i j i j i j η+-=+-+u u u u u (A.5) 为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零，需要引入一个与密度（或者压力）相关的开关函数：
11
11
i i i i i i i i ρρρρθρρρρ+-+----=-+- (A.6)
由式(A.6) 可以看出，在光滑区域，密度变化很缓，因此θ值也很小；而在激波附近密度变化很陡，θ值就很大。

带有开关函数的前置人工黏性滤波方法为：
(),,1,,1,122
n n
n n n i j i j
i j i j i j ηθ+-=+-+u u u u u (A.7) 其中参数η往往需要通过实际试算来确定，也可采用线性近似方法得到：
||||1t a t a x x η∆∆⎛⎫
=
- ⎪∆∆⎝⎭
(A.8)
由于声速不会超过3，所以取||3a =，在本计算中取0.25η=。

4.计算结果分析
计算分别采用标准的C 语言和Fortran77语言编写程序。

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

计算得到在0.4T =时刻的密度、速度和压力分布
如图A.2（C 语言计算结果）和图A.3（Fortran77语言计算结果）所示。

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

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

采用带开关函数的前置
人工滤波法能消除激波附近的非物理振荡，计算效果很好。

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

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

A-2 一维Riemann问题数值计算源程序
1. C语言源程序
// MacCormack1D.cpp : 定义控制台应用程序的入口点。

//
/*-----------------------------------------------------------------------------------------------------
*利用MacCormack差分格式求解一维激波管问题（C语言版本）
*
-------------------------------------------------------------------------------------------------------*/
//#include "stdafx.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define GAMA 1.4//气体常数
#define PI 3.141592654
图 A.2采用C语言程序得到的一维
Riemann问题密度、速度和压力分布
图 A.3采用Fortran77语言程序得到的一维
Riemann问题密度、速度和压力分布
#define L 2.0//计算区域
#define TT 0.4//总时间
#define Sf 0.8//时间步长因子
#define J 1000//网格数
//全局变量
double U[J+2][3],Uf[J+2][3],Ef[J+2][3];
/*------------------------------------------------------- 计算时间步长
入口: U，当前物理量，dx，网格宽度；
返回: 时间步长。

---------------------------------------------------------*/ double CFL(double U[J+2][3],double dx)
{
int i;
double maxvel,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;
}
return Sf*dx/maxvel;
}
/*------------------------------------------------------- 初始化
入口: 无；
出口: U，已经给定的初始值，
dx, 网格宽度。

---------------------------------------------------------*/ void Init(double U[J+2][3],double & dx)
{
int i;
double rou1=1.0 ,u1=0.0,p1=1.0; //初始条件
double rou2=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，已经给定的边界。

---------------------------------------------------------*/ void bound(double U[J+2][3],double dx)
{
int k;
//左边界
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方程组。

---------------------------------------------------------*/ void U2E(double U[3],double E[3])
{
double u,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;
}
/*-------------------------------------------------------
一维MacCormack差分格式求解器
入口: U，上一时刻的U矢量，Uf、Ef，临时变量，
dx，网格宽度，dt, 时间步长；
出口: U，计算得到的当前时刻U矢量。

---------------------------------------------------------*/
void MacCormack_1DSolver(double U[J+2][3],double Uf[J+2][3],double Ef[J+2][3],double dx,double dt)
{
int i,k;
double r,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，网格宽度；
出口: 无。

---------------------------------------------------------*/
void Output(double U[J+2][3],double dx)
{
int i;
FILE *fp;
double rou,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);
}
/*-------------------------------------------------------
主函数
入口: 无；
出口: 无。

---------------------------------------------------------*/
void main()
{
double T,dx,dt;
Init(U,dx);
T=0;
while(T<TT)
{
dt=CFL(U,dx);
T+=dt;
printf("T=%10g dt=%10g\n",T,dt);
MacCormack_1DSolver(U,Uf,Ef,dx,dt);
bound(U,dx);
}
Output(U,dx);
}
--------------------------------------------------------------------
--------------------------------
----------------------------------------------------------------------------------------------------
2. Fortran77语言源程序
! MacCormack1D.for
------------------------------------------------------------------------------------------------------------
!利用MacCormack差分格式求解一维激波管问题（Fortran77语言版本）
--------------------------------------------------------------------------------------------------------------*/
program MacCormack1D
implicit double precision (a-h,o-z)
parameter (M=1000)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:M+1,0:2),Uf(0:M+1,0:2)
dimension Ef(0:M+1,0:2)
!气体常数
GAMA=1.4
PI=3.1415926
!网格数
J=M
!计算区域
dL=2.0
!总时间
TT=0.4
!时间步长因子
Sf=0.8
call Init(U,dx)
T=0
1 dt=CFL(U,dx)
T=T+dt
write(*,*)'T=',T,'dt=',dt
call MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
call bound(U,dx)
if(T.lt.TT)goto 1
call Output(U,dx)
end
!------------------------------------------------------- !计算时间步长
!入口: U，当前物理量，dx, 网格宽度；
!返回: 时间步长。

!------------------------------------------------------- double precision function CFL(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
dmaxvel=1e-10
do 10 i=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
10 continue
CFL=Sf*dx/dmaxvel
end
!------------------------------------------------------- !初始化
!入口: 无；
!出口: U, 已经给定的初始值，dx，网格宽度。

!------------------------------------------------------- subroutine Init(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(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
do 20 i=0,J/2
U(i,0)=rou1
U(i,1)=rou1*u1
U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u1
20 continue
do 21 i=J/2+1,J+1
U(i,0)=rou2
U(i,1)=rou2*u2
U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2
21 continue
end
!------------------------------------------------------- !边界条件
!入口: dx，网格宽度；
!出口: U，已经给定边界。

!------------------------------------------------------- subroutine bound(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
!左边界
do 30 k=0,2
U(0,k)=U(1,k)
30 continue
!右边界
do 31 k=0,2
U(J+1,k)=U(J,k)
31 continue
end
!------------------------------------------------------- !根据U计算E
!入口: U，当前U矢量；
!出口: E，计算得到的E矢量，
! U、E定义见Euler方程组。

!------------------------------------------------------- subroutine U2E(U,E,is,in)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2),E(0:J+1,0:2)
do 40 i=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
40 continue
end
!-------------------------------------------------------
!一维MacCormack差分格式求解器
!入口: U，上一时刻U矢量，
! Uf、Ef，临时变量，
! dx，网格宽度，dt,，时间步长；
!出口: U，计算得到得当前时刻U矢量。

!-------------------------------------------------------
subroutine MacCormack_1D_Solver(U,Uf,Ef,dx,dt)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2),Uf(0:J+1,0:2)
dimension Ef(0:J+1,0:2)
r=dt/dx
dnu=0.25
do 60 i=1,J
do 60 k=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)) 60 continue
do 61 k=0,2
do 61 i=1,J
U(i,k)=Ef(i,k)
61 continue
call U2E(U,Ef,0,J+1)
do 63 i=0,J
do 63 k=0,2
!U(n+1/2)(i+1/2)
Uf(i,k)=U(i,k)-r*(Ef(i+1,k)-Ef(i,k))
63 continue
!E(n+1/2)(i+1/2)
call U2E(Uf,Ef,0,J)
do 64 i=1,J
do 64 k=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))
64 continue
end
!-------------------------------------------------------
!输出结果, 用Origin数据格式画图
!入口: U，当前时刻U矢量，
! dx，网格宽度；
!出口: 无。

!-------------------------------------------------------
subroutine Output(U,dx)
implicit double precision (a-h,o-z)
common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf
dimension U(0:J+1,0:2)
open(1,file='result.txt',status='unknown')
do 80 i=0,J+1
rou=U(i,0)
uu=U(i,1)/rou
p=(GAMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)
write(1,81)i*dx,rou,uu,p,U(i,2)
80 continue
close(1)
81 format(D20.10,D20.10,D20.10,D20.10,D20.10)
end
------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------
------------------------------------------------------------------------------------------------------------------。