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

1、一维 Riemann问题数值解与计算程序一维 Riemann问题，即激波管问题，是一个典型的一维可压缩无黏气体动力学问题，并有解析解。对它采用二阶精度MacCormack两步差分格式进行数值求解。同时，为了初学者入门和练习方便, 这里给出了用 C 语言和 Fortran77编写的计算一维 Riemann问题的计算程序，供大家学习参考。A-1 利用 MacCormack两步差分格式求解一维Riemann问题1. 一维 Riemann问题一维 Riemann问题实际上就是激波管问题。激波管是一根两端封闭、内部充满气体的直管， 如图 A.1 所示。在直管中由一薄膜将激波管隔开， 在薄膜两侧充有均匀理

2、想气体（可以是同一种气体，也可以是不同种气体），薄膜两侧气体的压力、密度不同。当 t 0 时，气体处于静止状态。当 t 0 时，薄膜瞬时突然破裂，气体从高压端冲向低压端，同时在管内形成激波、稀疏波和接触间断等复杂波系。2. 基本方程组、初始条件和边界条件图 A.1 激波管问题示意图设气体是理想气体。一维Riemann问题在数学上可以用一维可压缩无黏气体Euler方程组来描述。在直角坐标系下量纲为一的一维Euler方程组为：uf,1x1(A.1)tx0u其中uu , fu 2p(A.2)E(E p)u这里 、 u 、 p 、 E 分别是流体的密度、速度、压力和单位体积总能。理想气体状态方程：p1

3、e1E1u2v2（A.3）2初始条件： 1 1, u10, p1 1 ； 20.125, u20, p20.1。供参考边界条件： x1 和 x1 处为自由输出条件， u0u1 , uNuN 1 。3. 二阶精度 MacCormack差分格式MacCormack两步差分格式：n1unjr f jnf jn 1uj2(A.4)n 11nn 11n 1n122f j2uj2u juj2r f j 1其中 rt 。计算实践表明， MacCormack两步差分格式不能抑制激波附近非物x理振荡。因此在计算激波时，必须采用人工黏性滤波方法：uin,juin, j1uin1, j2uin, juin 1, j

4、(A.5)2为了在激波附近人工黏性起作用， 而在光滑区域人工黏性为零， 需要引入一个与密度（或者压力）相关的开关函数：i1iii1(A.6)i1iii1由式 (A.6) 可以看出，在光滑区域，密度变化很缓，因此 值也很小；而在激波附近密度变化很陡， 值就很大。带有开关函数的前置人工黏性滤波方法为：uin,juin, j1uin1, j 2uin, j uin 1, j(A.7)2其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到：t | a | 1t | a |(A.8)xx由于声速不会超过3，所以取 | a |3 ，在本计算中取0.25 。4. 计算结果分析计算分别采用标准的 C

5、语言和 Fortran77语言编写程序。计算中网格数取 1000 ，计算总时间为 T 0.4。计算得到在 T 0.4 时刻的密度、速度和压力分布如图 A.2 （ C 语言计算结果）和图 A.3 （ Fortran77语言计算结果）所示。采用两种不同语言编写程序所得到的计算结果完全吻合。从图 A.2 和图 A.3 中可以发现，MacCormack两步差分格式能很好地捕捉激波，计算得到的激波面很陡、 很窄，计算激波精度是很高的。 采用带开关函数的前置供参考人工滤波法能消除激波附近的非物理振荡，计算效果很好。从图 A.2 和图 A.3 中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布

6、中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的， 而密度和速度是相等的。 这个计算结果正确地反映了一维Riemann问题的物理特性，并被激波管实验所验证。图A.2采 用 C 语 言 程 序 得 到 的 一 维图 A.3采用 Fortran77语言程序得到的一维Riemann问题密度、速度和压力分布Riemann问题密度、速度和压力分布A-2 一维 Riemann问题数值计算源程序1.C 语言源程序/ MacCormack1D.cpp :定义控制台应用程序的入口点。/*-* 利用 MacCormack差分格式求解一维激波管问题（C 语言版本）*-*/#include

7、"stdafx.h"#include <stdio.h>#include <stdlib.h>#include <math.h>#define GAMA 1.4/ 气体常数#define PI3.141592654#define L2.0/ 计算区域供参考#define TT 0.4/ 总时间#define Sf 0.8/ 时间步长因子#define J1000/网格数/全局变量double UJ+23,UfJ+23,EfJ+23;/*-计算时间步长入口 : U ，当前物理量，dx，网格宽度；返回 : 时间步长。-*/double CFL

8、(double UJ+23,double dx)int i;double maxvel,p,u,vel;maxvel=1e-100;for(i=1;i<=J;i+)u=Ui1/Ui0;p=(GAMA-1)*(Ui2-0.5*Ui0*u*u);vel=sqrt(GAMA*p/Ui0)+fabs(u);if(vel>maxvel)maxvel=vel;return Sf*dx/maxvel;/*-初始化入口: 无；出口 : U ， 已经给定的初始值，dx, 网格宽度。-*/void Init(double UJ+23,double & dx)int i;double rou1=

9、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+)供参考Ui0=rou1;Ui1=rou1*u1;Ui2=p1/(GAMA-1)+rou1*u1*u1/2;for(i=J/2+1;i<=J+1;i+)Ui0=rou2;Ui1=rou2*u2;Ui2=p2/(GAMA-1)+rou2*u2*u2/2;/*-边界条件入口 : dx，网格宽度；出口 : U ， 已经给定的边界。-*/void bound(double UJ+23,double dx)int k;/左边界for

10、(k=0;k<3;k+)U0k=U1k;/右边界for(k=0;k<3;k+)UJ+1k=UJk;/*-根据 U计算 E入口:U，当前 U矢量；出口 : E， 计算得到的E 矢量，U、 E 的定义见Euler 方程组。-*/void U2E(double U3,double E3)double u,p;u=U1/U0;p=(GAMA-1)*(U2-0.5*U1*U1/U0);E0=U1;E1=U0*u*u+p;E2=(U2+p)*u;供参考/*-一维 MacCormack差分格式求解器入口 : U ， 上一时刻的U 矢量， Uf 、 Ef ，临时变量，dx，网格宽度， dt, 时间

11、步长；出口 : U ， 计算得到的当前时刻U 矢量。-*/void MacCormack_1DSolver(double UJ+23,double UfJ+23,double EfJ+23,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(Ui+10-Ui0)-fabs(Ui0-Ui-10)/(fabs(Ui+10-Ui0)+fabs(Ui0-Ui-10)+1e-100); /开关函数for(k=0;k<3;k+)Efik=Uik+0.5*nu*q*(Ui+1k

12、-2*Uik+Ui-1k);/人工黏性项for(k=0;k<3;k+)for(i=1;i<=J;i+)Uik=Efik;for(i=0;i<=J+1;i+)U2E(Ui,Efi);for(i=0;i<=J;i+)for(k=0;k<3;k+)Ufik=Uik-r*(Efi+1k-Efik); /U(n+1/2)(i+1/2)for(i=0;i<=J;i+)U2E(Ufi,Efi); /E(n+1/2)(i+1/2)for(i=1;i<=J;i+)for(k=0;k<3;k+)Uik=0.5*(Uik+Ufik)-0.5*r*(Efik-Efi-1

13、k); /U(n+1)(i)/*-输出结果 , 用 Origin 数据格式画图入口 : U ， 当前时刻U 矢量， dx， 网格宽度；出口: 无。供参考-*/void Output(double UJ+23,double dx)int i;FILE *fp;double rou,u,p;fp=fopen("result.txt","w");for(i=0;i<=J+1;i+)rou=Ui0;u=Ui1/rou;p=(GAMA-1)*(Ui2-0.5*Ui0*u*u);fprintf(fp,"%20f%20.10e%20.10e%20.10

14、e%20.10en",i*dx,rou,u,p,Ui2);fclose(fp);/*-主函数入口: 无；出口: 无。-*/void main()double T,dx,dt;Init(U,dx);T=0;while(T<TT)dt=CFL(U,dx);T+=dt;printf("T=%10gdt=%10gn",T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);-供参考2.Fortran77语言源程序! MacCormack1D.for-! 利用 MacCormack差分格式求解

15、一维激波管问题（ Fortran77语言版本）-*/program MacCormack1Dimplicit double precision (a-h,o-z)parameter (M=1000)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:M+1,0:2),Uf(0:M+1,0:2)dimension Ef(0:M+1,0:2)! 气体常数GAMA=1.4PI=3.1415926! 网格数J=M! 计算区域 dL=2.0! 总时间TT=0.4! 时间步长因子Sf=0.8call Init(U,dx)T=01 dt=CFL(U,dx)T=

16、T+dtwrite(*,*)'T=',T,'dt=',dtcall MacCormack_1D_Solver(U,Uf,Ef,dx,dt)call bound(U,dx)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,Sfdimension U(0:J+1,0

17、:2)dmaxvel=1e-10do 10 i=1,Juu=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)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,Sfdimension U(0:J+1,0

18、:2)! 初始条件 rou1=1.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=0.1dx=dL/Jdo 20 i=0,J/2U(i,0)=rou1U(i,1)=rou1*u1U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120 continuedo 21 i=J/2+1,J+1 U(i,0)=rou2供参考U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221 continue end!-! 边界条件! 入口 : dx ，网格宽度；! 出口 : U ， 已经给定边界。!-subroutine

19、 bound(U,dx)implicit double precision (a-h,o-z)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:J+1,0:2)! 左边界do 30 k=0,2U(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

20、 precision (a-h,o-z)common /G_def/ GAMA,PI,J,JJ,dL,TT,Sfdimension U(0:J+1,0:2),E(0:J+1,0:2)do 40 i=is,inuu=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+pE(i,2)=(U(i,2)+p)*uu40continue供参考end!-! 一维 MacCormack差分格式求解器!入口 : U， 上一时刻U 矢量，! Uf 、 Ef ，临时变量，! dx，网格宽

21、度， 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,Sfdimension U(0:J+1,0:2),Uf(0:J+1,0:2)dimension Ef(0:J+1,0:2)r=dt/dxdnu=0.25do 60 i=1,Jdo 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)