一维问题激波管问题数值解与其计算程序流程Word文件下载.doc

1、 （A.5）为了在激波附近人工黏性起作用，而在光滑区域人工黏性为零，需要引入一个与密度（或者压力）相关的开关函数： （A.6）由式（A.6） 可以看出，在光滑区域，密度变化很缓，因此值也很小；而在激波附近密度变化很陡，值就很大。带有开关函数的前置人工黏性滤波方法为： （A.7）其中参数往往需要通过实际试算来确定，也可采用线性近似方法得到：（A.8）由于声速不会超过3，所以取，在本计算中取。4.计算结果分析计算分别采用标准的语言和语言编写程序。计算中网格数取，计算总时间为。计算得到在时刻的密度、速度和压力分布如图A.2（语言计算结果）和图A.3（语言计算结果）所示。采用两种不同语言编写程序所得到

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

3、言源程序/ MacCormack1D.cpp : 定义控制台应用程序的入口点。/*-*利用差分格式求解一维激波管问题（语言版本）* -*/#include stdafx.h #include stdlib.hmath.h#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;/*-计算时间步长入口: U，当前物理量，dx，网格宽度；返回: 时间步长。-*

4、/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=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;

5、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; /右边界 k+）UJ+1k=UJk;根据U计算E U， 当前U矢量； E， 计算得到的E矢量，U、E的定义见Euler方程组。void U2E（double U3,double E3）

6、 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;一维差分格式求解器 U， 上一时刻的U矢量，Uf、Ef，临时变量， dx，网格宽度，dt, 时间步长； U， 计算得到的当前时刻U矢量。void MacCormack_1DSolver（double UJ+23,double UfJ+23,double EfJ+23,double dx,double dt） int i,k; double r,nu,q; r=dt/dx; nu=0.25; q=fabs（fabs（Ui+10-U

7、i0）-fabs（Ui0-Ui-10） /（fabs（Ui+10-Ui0）+fabs（Ui0-Ui-10）+1e-100）; /开关函数 for（k=0;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*（Efi+1k-Efik）; /U（n+1/2）（i+1/2）i+）U2E（Ufi,Efi）; /E（n+1/2）（i+1/2） for（k=0; Uik=0.5*（Uik+Ufik）-0.5*r*（Efik-Efi-1k）; /U（n+1

8、）（i） 输出结果, 用数据格式画图 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.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+=d

9、t; printf（T=%10g dt=%10gn,T,dt）; MacCormack_1DSolver（U,Uf,Ef,dx,dt）; bound（U,dx）; Output（U,dx）;-2. 语言源程序! MacCormack1D.for -利用差分格式求解一维激波管问题（语言版本） -*/ program MacCormack1D implicit double precision （a-h,o-z） parameter （M=1000） common /G_def/ GAMA,PI,J,JJ,dL,TT,Sf dimension U（0:M+1,0:2）,Uf（0:2） dimens

10、ion Ef（0:气体常数 GAMA=1.4 PI=3.1415926网格数 J=M计算区域 dL=2.0总时间 TT=0.4时间步长因子 Sf=0.8 call 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.TT）goto 1 call Output（U,dx） end- U， 当前物理量，dx, 网格宽度； double precision function CFL（U,dx）J+1,0: d

11、maxvel=1e-10 do 10 i=1,J uu=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 continue CFL=Sf*dx/dmaxvel U, 已经给定的初始值，dx，网格宽度。 subroutine Init（U,dx）初始条件 rou1=1.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=0.1 dx=dL/J do 20 i=0,J/2 U（i,0）

12、=rou1 U（i,1）=rou1*u1 U（i,2）=p1/（GAMA-1）+0.5*rou1*u1*u120 continue do 21 i=J/2+1,J+1 U（i,0）=rou2 U（i,1）=rou2*u2 U（i,2）=p2/（GAMA-1）+0.5*rou2*u2*u221 continue U， 已经给定边界。 subroutine bound（U,dx）左边界 do 30 k=0,2 U（0,k）=U（1,k）30 continue右边界 do 31 k=0,2 U（J+1,k）=U（J,k）31 continue U，当前U矢量； E，计算得到的E矢量， U、E定义见E

13、uler方程组。 subroutine U2E（U,E,is,in）2）,E（0: do 40 i=is,in 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+p E（i,2）=（U（i,2）+p）*uu40 continue U， 上一时刻U矢量， Uf、Ef，临时变量， dx，网格宽度，dt,，时间步长； U， 计算得到得当前时刻U矢量。 subroutine MacCormack_1D_Solver（U,Uf,Ef,dx,dt） r=dt/dx dnu=0.25 do 60

14、i=1,J do 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,J U（i,k）=Ef（i,k）61 continue call U2E（U,Ef,0,J+1） do 63 i=0,J do 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 continue E（n+1/2）（i+1/2） call U2E（Uf,Ef,0,J） do 64 i=1,J do 64 k=0,2U（n+1）（i） U（i,k）=0.5*（