《一维黎曼问题数值解与计算程序.doc》由会员分享，可在线阅读，更多相关《一维黎曼问题数值解与计算程序.doc（25页珍藏版）》请在三一办公上搜索。
1、一维黎曼问题数值解与计算程序，计算总时间为。计算得到在时刻的密度、速度和压力分布如图A.2（语言计算结果）和图A.3（语言计算结果）所示。采用两种不同语言编写程序所得到的计算结果完全吻合。从图A.2和图A.3中可以发现，两步差分格式能很好地捕捉激波，计算得到的激波面很陡、很窄，计算激波精度是很高的。采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡，计算效果很好。从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的；在压力分布中存在第二个台阶，表明在这里存在一个接触间断，在接触间断两侧压力是有间断的，而密度和速度是相等的。这个计算结果正确地反映了一维问题的物理特性，并
2、被激波管实验所验证。 图A.2采用语言程序得到的一维问题密度、速度和压力分布图A.3采用语言程序得到的一维问题密度、速度和压力分布A-2 一维问题数值计算源程序1. 语言源程序/ MacCormack1D.cpp : 定义控制台应用程序的入口点。/*-*利用差分格式求解一维激波管问题（语言版本）* -*/#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/时间步长因子#d
3、efine 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;
4、 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;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) in
5、t k; /左边界 for(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;/*-一维差分格式求解器入口: U， 上一时刻的U矢量，Uf、Ef，临时变量， dx，网格宽度，dt, 时间步长；出口: U， 计算得到
6、的当前时刻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; 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);/人
7、工黏性项 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) /*-输出结果, 用数据格式画图入口: U，
8、当前时刻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() doubl
9、e 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. 语言源程序! MacCormack1D.for -!利用差分格式求解一维激波管问题（语言版本） -*/ program MacCormack1D implicit double precision (a-h,o-z) parameter (M=1000) common /G_def/ G
10、AMA,PI,J,JJ,dL,TT,Sf dimension 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.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 Out
11、put(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-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(u
12、u) if(vel.gt.dmaxvel)dmaxvel=vel10 continue CFL=Sf*dx/dmaxvel end !-!初始化!入口: 无；!出口: 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.0 u1=0 v1=0 p1=1.0 rou2=0.125 u2=0 v2=0 p2=0.1 dx=dL/J do 20 i=
13、0,J/2 U(i,0)=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 end!-!边界条件!入口: dx，网格宽度；!出口: U， 已经给定边界。!- subroutine bound(U,dx) implicit double precision (a-h,o-z) common /G_def/ GAMA,PI,J,JJ,dL,T
14、T,Sf dimension U(0:J+1,0:2)!左边界 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 end!-!根据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:J+1,0:2),E(0:J
15、+1,0:2) do 40 i=is,in uu=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+p E(i,2)=(U(i,2)+p)*uu40 continue end!-!一维差分格式求解器!入口: U， 上一时刻U矢量，! Uf、Ef，临时变量，! dx，网格宽度，dt,，时间步长；!出口: U， 计算得到得当前时刻U矢量。!- subroutine MacCormack_1D_Solver(U,Uf,Ef,dx,dt) implicit dou
16、ble 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/dx dnu=0.25 do 60 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
17、*(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,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,J do 64 k=0,2!U(n+1)(i) U(i,k)=0.5*(U(i,k)+Uf(i,k)-0.5
18、*r*(Ef(i,k)-Ef(i-1,k) 64 continue end!-!输出结果, 用数据格式画图!入口: 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+1 rou=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 continue close(1) 81 format(D20.10,D20.10,D20.10,D20.10,D20.10) end-
