一维黎曼问题

《一维黎曼问题》由会员分享，可在线阅读，更多相关《一维黎曼问题（8页珍藏版）》请在金锄头文库上搜索。

1、图A.1激波管问题示意图其中du dfdt dx=0, -1 x 1/ p/pu 、u=pu,f=pu 2+pE J.(E + p)u J(A.1)(A.2)附录A 一维Riemann问题数值解与计算程序一维Riemarn问题，即激波管问题，是一个典型的一维可压缩无黏气体动力 学问题，并有解析解。对它采用二阶精度MacCormac两步差分格式进行数值求 解。同时，为了初学者入门和练习方便，这里给出了用C语言和Fortran7编写的 计算一维Riemamf问题的计算程序，供大家学习参考。A-1利用 MacCormack两步差分格式求解一维Rieman 问题1. 一维Riemam问题一维Riema

2、m问题实际上就是激波管问题。激波管是一根两端封闭、内部充 满气体的直管，如图A.1所示。在直管中由一薄膜将激波管隔开，在薄膜两侧充 有均匀理想气体（可以是同一种气体， 也可以是不同种气体），薄膜两侧气体 的压力、密度不同。当1 0时，气体 处于静止状态。当t = 0时，薄膜瞬时 突然破裂，气体从高压端冲向低压端， 同时在管内形成激波、稀疏波和接触间 断等复杂波系。2.基本方程组、初始条件和边界条件设气体是理想气体。一维Riemam问题在数学上可以用一维可压缩无黏气体Eulers程组来描述。在直角坐标系下量纲为一的一维Euler方程组为:这里p、u、e分别是流体的密度、速度、压力和单位体积总能。

3、理想气体状态方程:p =(Y -1)pe =(Y -1)(u 2 + v 2)（A.3）初始条件：p = 1, u = 0, p = 1 ； p = 0.125, u = 0, p = 0.1。111222边界条件：x = -1和x = 1处为自由输出条件，u =, UN = UN 13.二阶精度MacCormac差分格式MacCormac两步差分格式:1Un+ 2 = Un - r jj (1)j j )if n - fj1-r2n )_ j-1(1fn+2 - fj+1(A.4)其中r =翌。计算实践表明，MacCormac两步差分格式不能抑制激波附近非物 Ax理振荡。因此在计算激波时，必

4、须采用人工黏性滤波方法：(A.5)1 （,喝）U n = Un +q Wn - 2Un + Un 7i, ji, j 2i+1，ji, j i-1,j为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零，需要引入一个与密度（或者压力）相关的开关函数:EzP I-|P Ip -p.l+lp.i+1iiS_1i - Pi-1(A.6)由式（A.6）可以看出，在光滑区域，密度变化很缓，因此9值也很小；而在激波 附近密度变化很陡，9值就很大。带有开关函数的前置人工黏性滤波方法为:(A.7) 一 ，1，一）U n = Un +q9 Un 一 2U n + Un 7i , j i , j 2i +1,

5、j i , j i - 1, j其中参数门往往需要通过实际试算来确定，也可采用线性近似方法得到:(A.8)At I a I 七 At I a l) 1Ax vAx )由于声速不会超过3，所以取I a l= 3，在本计算中取q = 0.25。4.计算结果分析计算分别采用标准的C语言和Fortran77语言编写程序。计算中网格数取 1000，计算总时间为T = 0.4。计算得到在T = 0.4时刻的密度、速度和压力分布 如图A.2（ C语言计算结果）和图A.3（ Fortran77语言计算结果）所示。采用两 种不同语言编写程序所得到的计算结果完全吻合。从图A.2和图A.3中可以发现，MacCorm

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

7、iemann问题数值计算源程序1. C语言源程序/ MacCormacklD.cpp :定义控制台应用程序的入口点。/*利用MacCormack差分格式求解一维激波管问题（C语言版本）*/#include stdafx.h#include #include #include #define GAMA 1.4/气 体常数#define 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;/*计算时间步长入

8、口： 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 rou1=1.0 ,u1=0.0,p1=1.0; /初始条件 double rou2=0.125,u2=0.0,p2=0.1;d

9、x=L/J;for(i=0;i=J/2;i+)Ui0=rou1;Ui1=rou1*u1;Ui=p1/(GAMA-1)+rou1*u1*u1/2;for(i=J/2+1;i=J+1;i+)Ui0=rou2;Ui1=rou2*u2;Ui H2 =p2/(GAMA-1)+rou2*u2*u2/2;/*边界条件入口： dx,网格宽度；出口： U,已经给定的边界。*/void bound(double UJ+23,double dx)int k;/佐边界for(k=0;k3;k+)U0k=U1k;/右边界for(k=0;k3;k+)UJ+1k=UJk;/*根据U计算E入口： U，当前U矢量；出口： E，

10、计算得到的E矢量，U、E的定义见Euler方程组。*/void U2E(double U3,double E3)double u,p;u=U1/U0;p=(GAMA-1)*(U -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,double UfJ+23,double EfJ+23,double dx,double

11、 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);/人工黏性项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+)U

12、fik=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数据格式画图入口： U，当前时刻U矢量，dx，网格宽度； 出口：无。*/void Output(double UJ+23,double dx)int i;FILE *fp;double rou,u,p;fp=fopen(result.txt

13、,w);for(i=0;i=J+1;i+)rou=Ui0;u=Ui1/rou;p=(GAMA-1)*(Ui-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(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);