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

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

2、理想气体（可以是同一种气体， 也可以是不同种气体），薄膜两侧气体 的压力、密度不同。当t乞0时，气体 处于静止状态。当t = 0时，薄膜瞬时 突然破裂，气体从高压端冲向低压端， 同时在管内形成激波、稀疏波和接触间 断等复杂波系。2.基本方程组、初始条件和边界条件/图A.1激波管问题示意图设气体是理想气体。一维Riema nri可题在数学上可以用一维可压缩无黏气体Eulei方程组来描述。在直角坐标系下量纲为一的一维 Euler方程组为:+f = 0, 1 兰 x 兰1 （A.1）.:t :xr P、Pu 、其中u =Pu,f =Pu2 + p(A.2)El(E + p)u这里、u、p、E分别是流

3、体的密度、速度、压力和单位体积总能。理想气体状态方程：p= -1 -1 E - J u2 v2 （ A.3）初始条件： J = 1,5 = 0,山=1 ； = 0.125, u2 = 0, p2 = 0.1。边界条件：X - -1和X=1处为自由输出条件，UohUUnhUn/。3.二阶精度MacCormac差分格式n -2MacCormac两步差分格式:2 n n - nUj -Uj -r fj - fjj其中r二卫。计算实践表明，MacCormac两步差分格式不能抑制激波附近非物Ax理振荡。因此在计算激波时，必须采用人工黏性滤波方法:(A.5)需要引入一个与-n n 1 n Q n nUi,

4、j =Ui,j - U1,j -2嘔-UiJ,j为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零,密度（或者压力）相关的开关函数:由式（A.6）可以看出，在光滑区域，密度变化很缓，因此 二值也很小；而在激波 附近密度变化很陡，V值就很大。带有开关函数的前置人工黏性滤波方法为:(A.7)_ 1Uinj =unj +罗日（ut - 2uinj + ulj其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到:(A.8)由于声速不会超过3,所以取|a|=3，在本计算中取 =0.25。4.计算结果分析计算分别采用标准的 C语言和Fortran77语言编写程序。计算中网格数取 1000，计算

5、总时间为T =0.4。计算得到在T =0.4时刻的密度、速度和压力分布如图A.2 （ C语言计算结果）和图A.3 （ Fortran77语言计算结果）所示。采用两种不同语言编写程序所得到的计算结果完全吻合。从图A.2和图A.3中可以发现，MacCormac两步差分格式能很好地捕捉激波, 计算得到的激波面很陡、很窄，计算激波精度是很高的。采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡，计算效果很好从图A.2和图A.3中可以看出通过激波后气体的密度、 压力和速度都是增加 的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断 两侧压力是有间断的，而密度和速度是相等的。这个

6、计算结果正确地反映了一维 Riema nn问题的物理特性，并被激波管实验所验证。A-2 一维Riemanr问题数值计算源程序1. C语言源程序/ MacCormackID.cpp :定义控制台应用程序的入口点。/* *利用MacCormac差分格式求解一维激波管问题（C语言版本）* */#in clude stdafx.h#in clude #in clude #in clude #define GAMA 1.4 气体常数#define PI 3.141592654#defi ne L 2.0计算区域#define TT 0.4/ 总时间#define Sf 0.8/ 时间步长因子#defin

7、e 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;imaxvel)maxvel=vel;return Sf*dx/maxvel;/* 初始化入口 ： 无；出口 ： U， 已经给定的初始值，dx, 网格宽度。 */void Init(double UJ+23,double & dx)int i;double

8、 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*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;/左边界 fo

9、r(k=0;k3;k+)U0k=U1k;/右边界for(k=0;k3;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;/* 一维MacCormac差分格式求解器入口 ： u, 上一时刻的 u矢量，Uf、Ef,临时变量，dx，网格宽度，dt,时间步长；出口

10、 ： u, 计算得到的当前时刻 u 矢量。 */EfJ+23,doublevoid MacCormack_1DSolver(double uJ+23,double ufJ+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;k3;k+)Efik=uik+0.5*nu*q*(ui+1k-2*uik+ui-1k

11、);/ 人工黏性项 for(k=0;k3;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;k3;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;k3;k+) uik=0.5*(uik+ufik)-0.5*r*(Efik-Efi-1k); /u(n+1)(i)/* 输出结果 , 用 Origin 数据格式画图

12、入口 ： 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.10e%20.10en,i*dx,rou,u,p,Ui2); fclose(fp);/* 主函数 入口 : 无； 出口 : 无。 */void main(

13、)double T,dx,dt;Init(U,dx);T=0; while(TTT) dt=CFL(U,dx);T+=dt;printf(T=%10g dt=%10gn,T,dt);MacCormack_1DSolver(U,Uf,Ef,dx,dt);bound(U,dx);Output(U,dx);2. Fortran7语言源程序! MacCormack1D.for!利用MacCormac差分格式求解一维激波管问题(Fortran77语言版本) */ program MacCormack1D implicit double precision (a-h,o-z) parameter (M=1

14、000) 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.4 PI=3.1415926! 网格数J=M! 计算区域dL=2.0! 总时间TT=0.4! 时间步长因子Sf=0.8call Init(U,dx)T=01 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.T

15、T)goto 1call 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-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(G

16、AMA*p/U(i,0)+dabs(uu) if(vel.gt.dmaxvel)dmaxvel=vel10 continueCFL=Sf*dx/dmaxvelend! !初始化!入口 : 无；!出口 ： 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.0u1=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=0.1dx=dL/Jd

17、o 20 i=0,J/2U(i,0)=rou1 U(i,1)=rou1*u1U(i,2)=p1/(GAMA-1)+0.5*rou1*u1*u120 continuedo 21 i=J/2+1,J+1U(i,0)=rou2U(i,1)=rou2*u2U(i,2)=p2/(GAMA-1)+0.5*rou2*u2*u2 21 continueend! !边界条件!入口 ： dx，网格宽度；! 出口 : U ， 已经给定边界。! subroutine bound(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,

18、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 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

19、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)*uu40 continueend! ! 一维MacCormac差分格式求解器! 入口 : U ， 上一时刻 U 矢量，! Uf 、 Ef ，临时变量，! dx，网格宽度，dt,时间步长；! 出口 : U ， 计算得到得当前时刻 U 矢量。! subroutine MacCormack_1D_So

20、lver(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/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) ! 人

21、工黏性项Ef(i,k)=U(i,k)+0.5*dnu*q*(U(i+1,k)-2*U(i,k)+U(i-1,k)60continuedo 61 k=0,2do 61 i=1,JU(i,k)=Ef(i,k)61continuecall U2E(U,Ef,0,J+1)do 63 i=0,Jdo 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,Jdo 64 k=0,2!U(n+1)(i)U(i,k)=0.5*(U(i,k

22、)+Uf(i,k)-0.5*r*(Ef(i,k)-Ef(i-1,k)64 continueend! !输出结果 , 用 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+1rou=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 continueclose(1)81 format(D20.10,D20.10,D20.10,D20.10,D20.10) end