一维黎曼问题数值解与计算程序Word格式文档下载.docx

1、El（E + p）u这里、u、p、E分别是流体的密度、速度、压力和单位体积总能。理想气体状态方程：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）需要

2、引入一个与-n n 1 n Q n nUi,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

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

4、侧压力是有间断的，而密度和速度是相等的。这个计算结果正确地反映了一维 Riema nn问题的物理特性，并被激波管实验所验证。A-2 一维Riemanr问题数值计算源程序1. C语言源程序/ MacCormackID.cpp :定义控制台应用程序的入口点。/* *利用MacCormac差分格式求解一维激波管问题（C语言版本）* */#in clude stdafx.h#in clude stdlib.hmath.h #define GAMA 1.4 气体常数#define PI 3.141592654#defi ne L 2.0计算区域#define TT 0.4/ 总时间#define Sf

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

6、=1.0 ,u1=0.0,p1=1.0; / 初始条件double rou2=0.125,u2=0.0,p2=0.1;dx=L/J;for（i=0;=J/2;Ui0=rou1;Ui1=rou1*u1; Ui2=p1/（GAMA-1）+rou1*u1*u1/2;for（i=J/2+1;=J+1;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;k3;k+）U0k=U1k;

7、/右边界for（k=0;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,时间步长； u, 计算得到的当前时刻 u 矢量。 */EfJ+23,double

8、void 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;q=fabs（fabs（ui+10-ui0）-fabs（ui0-ui-10） /（fabs（ui+10-ui0）+fabs（ui0-ui-10）+1e-100）; / 开关函数k+）Efik=uik+0.5*nu*q*（ui+1k-2*uik+ui-1k）;/ 人工黏性项 k+） for（i=1;i+）uik=Efik;i+）u2E（ui,Efi）;k+） ufik=uik-r

9、*（Efi+1k-Efik）; /u（n+1/2）（i+1/2）i+）u2E（ufi,Efi）; /E（n+1/2）（i+1/2）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）FILE *fp;double rou,u,p;fp=fopen（result.txt,w）;rou=Ui0;u=Ui1/rou;fprintf（fp,%20f%20.10e%20.10e%20

10、.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）;2. Fortran7语言源程序! MacCormack1D.for利用MacCormac差分格式求解一维激波管问题（Fortran77语言版本） */ program Mac

11、Cormack1D implicit 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:2） dimension Ef（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_Solv

12、er（U,Uf,Ef,dx,dt） call bound（U,dx）if（T.lt.TT）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:dmaxvel=1e-10do 10 i=1,Juu=U（i,1）/U（i,0） p=（GAMA-1）*（U（i,2）-0

13、.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入口 : 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： 初始条件rou1=1.0u1=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=0.1

14、dx=dL/Jdo 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 continue 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：

15、左边界do 30 k=0,2U（0,k）=U（1,k）30 continue 右边界do 31 k=0,2U（J+1,k）=U（J,k）31 continue 根据 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：2）,E（0：do 40 i=is,inuu=U（i,1）/U（i,0） p=（GAMA-1

16、）*（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 continue 一维MacCormac差分格式求解器 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

17、,TT,Sf dimension U（0: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）60continuedo 61 k=0,2do 61 i=1,JU（i,k）=Ef（i,k）61continuecall U2E（U,Ef,0,J+1）d

18、o 63 i=0,Jdo 63 k=0,2U（n+1/2）（i+1/2）Uf（i,k）=U（i,k）-r*（Ef（i+1,k）-Ef（i,k）63 continueE（n+1/2）（i+1/2）call U2E（Uf,Ef,0,J）do 64 i=1,Jdo 64 k=0,2U（n+1）（i）U（i,k）=0.5*（U（i,k）+Uf（i,k）-0.5*r*（Ef（i,k）-Ef（i-1,k）64 continue入口: 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: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