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

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

2、/Z/7/77/检，%文档边界条件：x 1和x 1处为自由输出条件，UoU1 , UN UN 1。3.二阶精度MacCormacki分格式MacCormac I断步差分格式:n1Uj2n11nuj-ujj2jnUj r1 n _uj2fjnfjnnfj11211 n _ fj2(A.4)其中r 0计算实践表明，MacCormac新步差分格式不能抑制激波附近非物 x理振荡。因此在计算激波时，必须采用人工黏性滤波方法：一n nUUui ,jui ,j为了在激波附近人工黏性起作用,1 n 。n n-ui 1,j 2ui,jui 1,j而在光滑区域人工黏性为零,(A.5)需要引入一个与密度（或者压力）

3、相关的开关函数:(A.6)由式（A.6）可以看出，在光滑区域，密度变化很缓，因此值也很小；而在激波附近密度变化很陡，值就很大。带有开关函数的前置人工黏性滤波方法为：(A.7)(A.8)-nn 1n Q n nui,jui,j -ui 1,j 2ui,jui 1,j其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到:t|a | 1 t|a|由于声速不会超过3,所以取|a| 3,在本计算中取0.254.计算结果分析计算分别采用标准的 C语言和Fortran77g言编写程序。计算中网格数取 1000,计算总时间为T 0.4。计算得到在T 0.4时刻的密度、速度和压力分布 如图A.2 （C语

4、言计算结果）和图A.3 （ Fortran7M言计算结果）所示。采用两 种不同语言编写程序所得到的计算结果完全吻合。从图A.2和图A.3中可以发现，MacCormac新步差分格式能很好地捕捉激波， 计算得到的激波面很陡、很窄，计算激波精度是很高的。采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡，计算效果很好。从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断 两侧压力是有间断的，而密度和速度是相等的。这个计算结果正确地反映了一维 Riemann诃题的物理特性，并被激波管实验所验证。图A.2采用C

5、语言程序得到的一维 Riemann题密度、速度和压力分布图A.3采用Fortran7旃言程序得到的一维 Riemann可题密度、速度和压力分布A-2 一维Riemann题数值计算源程序1. C语言源程序/ MacCormackID.cpp :定义控制台应用程序的入口点。/*# 利用MacCormackE分格式求解一维激波管问题（C语言版本）*/#include "stdafX.h"#include <stdio.h>#include <stdlib.h>#include <math.h># define GAMA 1.4 气体常数#def

6、ine PI 3.141592654# define L 2.0计算区域# define TT 0.4/总时间# define Sf 0.8时间步长因子# define J 1000网格数/全局变量double UJ+23,UfJ+23,EfJ+23;/*计算时间步长入口： U ,当前物理量，dx,网格宽度； 返回：时间步长。*/double CFL(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);

7、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=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*u

8、1*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(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

9、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,时间步长； 出口： U ,计算得到的当前时刻 U矢量。*/void MacCormack_1DSolver(double UJ+23,doubleEfJ+23,double dx,double dt) int i,k;double r,nu,q;r=dt/dx;nu=0.25;for(i=1;i<=J;i+)

10、 q=fabs(fabs(Ui+10-Ui0)-fabs(Ui0-Ui-1。) /(fabs(Ui+10-Ui0)+fabs(Ui0-Ui-10)+1e-100); for(k=0;k<3;k+)Efik=Uik+0.5*nu*q*(Ui+1k-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)f

11、or(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-1k); U(n+1)(i) /*UfJ+23,double开关函数人工黏性项输出结果，用Origin数据格式画图入口： U,当前时刻U矢量，dx,网格宽度；出口 ：无。*/void Output(double UJ+23,double dx) int i;FILE *fp;double rou,u,p;fp=fopen("result.txt&qu

12、ot;,"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.10e%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=%10g dt=%10gn&q

13、uot;,T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);2. Fortran77s言源程序! MacCormackID.forFortran77§言版本)!利用MacCormacki分格式求解一维激波管问题（*/program MacCormackIDimplicit 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)d

14、imension 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=T+dtwrite(*,*)'T=',T,'dt=',dtcall MacCormack_1D_Solver(U,Uf,Ef,dx,dt) call bound(U,dx)if(T.lt.TT)goto 1call Output(U,dx)end计算时间步长入口： U,当前物理量，dx,网格宽度;返回：时间步长。doubl

15、e 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-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) if(vel.gt.dmaxvel)dmaxvel=vel10 continueCFL=Sf*dx/dmaxvelend!初始化

16、!入口 ：无；!出口： 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:2)!初始条件rou1=1.0u1=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=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+

17、1,J+1U(i,0)=rou2U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u221 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,2U(0,k)=U(1,k)30 continue!右边界do 31 k=0,2U(J+1,k尸U(J,k)31

18、 continueend根据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,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*

19、uu+pE(i,2)=(U(i,2)+p)*uu40 continueend一维MacCormackE分格式求解器入口： 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,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) 60 continue do 61 k=0,2 do 61 i=1,JU(i,k尸Ef(i,k