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

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

2、想气体(可以是同一种气体,也可以是不同种气体),薄膜两侧气体Eulei方程组来描述.在直角坐标系卜晟纲为一的一维Eulei-方程组为：理+堂=0,-13 幻zA 1dt dx1§'' pllW= pll k= + P(A. 2)这里p、/(E分别是流体的密度、速度.压力和单位体积总能.理想气体状态方 程：初始条件:Pi = L W)=O,p, = 1 ； p2 =0.125, “2 = 0, P2 = 01o边界条件：x = -1和X = 1处为自由输出条件,“O =", ZZv = Wv3 .二阶精度MacConnack差分格式MacCoimack两步差分

3、格式：(A. 4)其中计算实践说明,MacConnack两步差分格式不能抑制激波附近非物Ax 理振荡.因此在计算激波时,必须采用人工黏性滤波方法：歌吆+如(心-2'W)(A. 5)为了在激波附近人工黏性起作用,而在光滑区域人工粘性为零,需要引入一个与密 度(或者压力)相关的开关函数：(A. 6)-二 IQ+II几】一加+IQ-久II由式(A.6)可以看出,在光滑区域,密度变化很缓,因此伙值也很小；而在激波附近 密度变化很陡,0值就很大.带有开关函数的前置人工黏性滤波方法为：磅F十如&仙-2炉血)(A. 7)其中参数往往需要通过实际试算来确定,也可采用线性近方法得到:(A. 8)

4、由于声速不会超过3所以取|=3,在本计算中取“二0.25.4 .计算结果分析计算分别采用标准的C语言和Foikan77语言编写程序.计算中网格数取1000,计算总时问为T = 0.4o计算得到在T = 0.4时刻的密度、速度和压力 分布如 图A. 2 C语言计算结果和图A.3 Fonian力语言计算结果所示.采用两种不 同语言编写程序所得到的计算结果完全吻合.从图A. 2和图A. 3中可以发现,MacCormack两步差分格式能很好地捕捉激 波, 计算得到的激波面很陡、很窄,计算激波精度是很高的.采用带开关函数的前置人工 滤波法能消除激波附近的非物理振荡,计算效果很好.从图A. 2和图A. 3

5、中可以看出通过激波后气体的密度压力和速度都是增加的； 在压力分布中存在第二个台阶,说明在这里存在一个接触间断,在接触间断两侧压力 是有间断的,而密度和速度是相等的.这个计算结果正确地反映了 一维Riemann问题 的物理特性,并被激波管实验所验证.值应力*i<K图A 2米用C语言程图A3采用Foilran77语言1. c语言源程序A-2 一维Riemaim问题数值计算源程序I I MacCormacklD.cpp :定义限制台应用程序的人 口点./*利用MacCoimack差分格式求解一维激波管问题C语言版本*/#include nstdafx.hH#include <stdio.

6、h>#include <stdlib.h>#include <math.h>#define GAMA 14 气体常数#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;计算时间步长人口当前物理量,dx,网格宽度;:时间步长.double CFL(double UJ+23,double dx)(int i;double maxvel , p , u , vel;

7、maxvel=le-100;for(i=l ;i<=J ;i+)u=Uil/Ui0;p=(GAMA-l)*(Ui2-0.5*Ui0*u*u);vel=sqrt(GAMAAp/UirO)+fabs(u);if(vel>maxvel)maxvel=vel;)return Sf*dx/maxvel;)初始化入口 ：无;出n : u,已经给定的初始值,dx,网格宽度勺void Init(double UJ+23,double & dx)(int i;double roul=1.0 9ul=0.0,pl=1.0; 初始条件double rou2=0.125,u2=0.0,p2=0.1

8、;dx=L/J;for(i=0;i<=J/2;i+)Ui0=roul;Uil=roulAul;Ui2=pV(GAMA-I)+roulAulAuV2;)for(i=J/2+l ;i<=J+l ;i+)(Ui0=rou2;Uil=rou2Au2;Ui2=p2/(GAMA-l)+rou2*u2*u2/2;)边界条件入口： dx,网格宽度出口 ： U,已经给定的边界.勺void bound(double UJ+239double dx)(int k;左边界for(k=0;k<3;k+)U0k=Ulk;右边界for(k=0;k<3;k+)UJ+lk=UJk;根据U计算E入口 当前

9、U矢量；出口： E,计算得到的E矢量,u、E的定义见Euler方程组.void U2E(double U3,double E3)(double u , p;u=UI/U0;p=(GAMA-l)*(U2-0.5*Ul*Ul/U0);EO=U1J;El=U0*u*u+p;E2=(U(2)+p)*u;*维MacConnack差分格式求解器An : 出上一时刻的IT矢量,uf> Ef,临时变量,出口 ： U,计算得到的当前时刻U矢量.dx.网格宽度dt,时间步长;勺 voidMacCormack_lDSolver(doubleUJ+2J3J4oubleUfJ+2(3), doubleEfJ+2

10、3,double dxAdouble dt)(int i , k;double r,nu,q;r=dt/dx;nu=0.25;for(i=l ;i<=J ;i+) q=fabs(fabs(Ui+lO-UiO)-fabs(UiO-Ui-l0)/(fabs(Ui+lO-UiO)+fabs(UiO-Ui-lO)+le-100); 开关函数for(k=0;k<3;k+)Efik=Uik+0.5*nu*q*(Ui+lk-2*Uik+Ui-lk);/A 工黏性项)for(k=0; k v3 ;k+)for(i=l ;i<=J ;i+)Ui k=Efi k;for(i=0;i<=J+

11、l;i+)U2E(Ui,Efi);for(i=0;i<=J ;i+)for (k=0 ;k<3 ;k+)Ufik=Uik-r*(Efi+lk-Efik);/U(n+V2)(i+V2)for(i=0;i<=J;i+)U2E(Ufi,Efi);/E(n+V2)(i+V2)for(i=l;i<=J;i+)for(k=0;k<3;k+)Uik=0.5*(UikWfik)-0.5*r*(Efik-Ef i-lk); /U(n+l)(i)输出结果,用6晌数据格式画图An : U,当前时刻u矢量,dx,网格宽度；出口 ：无.勺void Output(double UJ+239d

12、ouble dx) inti;FILE *fp;double rou5u,p;fp=fopen(ft result.txtfor(i=0;i<=J+l ;i+)rou=Ui0;u=Uril/rou;p=(GAMA-l)*(Ui2-0.5*Ui0*u*u);fprintf(fp,H%20f%20.10e%20.10e%20.10e%20.1 Oenn4 *dx,rou,u,p,U i 2);)fclose(fp);)/*主函数入口 ：无;出口 ：无.*/void main()(double T , dx , dt;Init(U,dx);T=0;whi!e(T<TT)dt=CFL(U,

13、dx);T+=dt;printf(nT= % 10g dt=%l 0gnM,dt);MacCormack_lDSolver(U9Uf,Ef?dx9dt); bound(U,dx);)Output(UAdx);2. Fortiun77语言源程序I MacCormacklD.for!利用MacCormack差分格式求解一维激波管问题(Fortian77语言版本)勺program MacCormacklDimplicit double precision (a-h,o-z)parameter (M=1000)common/G_de GAMA , PI , J , JJ , dL , TT , Sfd

14、imension U(O:M+1,0:2)JJf(0:M+l , 0(2) dimension Ef(0:M+l,0:2) !气体常数 GAMA=1.4PI=3.1415926J=M!计算区域dL=2.0!总时间TT=0.4!时间步长因子Sf=0.8call Init(U,dx)T=01 dt=CFL(U,dx)T=T+dtwrite(*A),T=T ; dt=dtcallMacCormack_lD_Solver(U?Uf,Ef9dx,dt) call bound(U,dx)if(T.lt.TT)goto 1 call Output(U,dx) end计算时间步长人口当前物理量,dx,网格宽度

15、；返回：时间步长.double precision function CFL(U,dx) implicit double precision (a-h,o-z) common/G_de GAMA , PI , J , JJ , dL , TT , Sfdimension U(0:J+M)：2) dmaxvel=le-10do 10 i=ljuu=U(i9l)/U(i,0)p=(GAMA-l)*(U(t2)-0.5*U(i,0)*uu*uu)vel=dsqrt(GAMAAp/U(i,O)+dabs(uu) if(veLgt.dmaxvel)dmaxvel=vel10 continueCFL=Sf

16、*dx/dmaxvelend初始化入口 ：无；出口已经给定的初始值,dx,网格subroutine Init(U,dx)implicit double precision (ah,oz)common/G_de GAMA , PI , J , JJ , dL , TT , Sfdimension U(0:J+l90:2)!初始条件roul=1.0ul=0vl=0pl=1.0rou2=0.125u2=0v2=0p2=0.1dx=dL/Jdo 20 i=0J/2U(i,0)=roulU(i,l)=roul*ulU(i,2)=pl/(GAMA-l)+0.5 A roul*ulA ul20 contin

17、uedo 21 i=J/2+l J+lU(i,0)=rou2U(i,l)=rou2*u2U(i,2)=p2/(GAMA-l)+0.5Arou2*u2Au221 continueend!边界条件!入口 ： dx,网格宽度;subroutine bound(U,dx) implicit double precision (a-h,o-z) common/G_de GAMA , PI , J , JJ , dL , TT , Sfdimension U(0:J+l,0:2)听边界do 30 k=0,2U(o , k)=U(l , k)30 continue!右边界do 31 k=0,2 U(J+l,

18、k)=U(J,k)31 continue end!根据U计算E!入口 ： U,当前U矢量；!出口 ：E,计算得到的e矢量,U、E定义见Euler方程组.subroutine LJ2E(U , E , is , in)implicit double precision (a-h,o-z) common/G_de GAMA , PI , J , JJ , dL , TT , Sfdimension U(0:J+l,0:2)9E(0:J+l,0:2)do 40 i=is9inuu=U(i4)/U(i,0) p=(GAMA-l)*(U(i,2)$ 0- 5*U(i, l)*U(i , 1)/U(i ,

19、 0)E(i,0)=U(i,l)E(i , l)=U(i , 0)*uu*iiu+pE(i , 2)=(U(i , 2)+p)*uu40 continue end!一维cCormack差分格式求解器！入II ： U,上一时亥U U 矢量,Uf>Ef,临时变量, ! dx,网格宽度,出时间步长； !出口： U,计算得到得当前时刻U矢量.subroutineMacCormack_lD_Solver(U?Uf,Ef9dx,dt) implicit double precision (a-h,o-z) common /G_defZGAMA , PI , J , JJ , dL , TT , Sf

20、dimension U(0: J+l90:2),Uf(0: J+l ,0:2) dimension Ef(0:J+l,0:2)r=dt/dx dnu=0.25do 60 i=lj do 60 k=0,2!开关函数q=dabs(dabs(U(i+hO)-U(i,O)-dabs(U(i,O)-U(i-l,O)$ /(dabs(U (i+1,0) -U (i,O)+dabs(U (i,0) -U (i-1,0) +le-10)!人工黏性项Ef(i, k)=U(i , k)+0- 5*dnu*q*(U(i+l , k)- 2*U(i , k)+U(i-hk)60 continuedo 61 k=0,2d