利用MacCormack两部差分格式求解一维激波管问题Fortran程序

《利用MacCormack两部差分格式求解一维激波管问题Fortran程序》由会员分享，可在线阅读，更多相关《利用MacCormack两部差分格式求解一维激波管问题Fortran程序（4页珍藏版）》请在装配图网上搜索。

1、!MacCormack1D.for!利用MacCormack两部差分格式求解一维激波管问题!programMacCormack1Dimplicitdoubleprecision(a-h,o-z)parameter(M=1000)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:M+1,0:2),Uf(0:M+1,0:2)dimensionEf(0:M+1,0:2)GAMMA=1.4!气体常数PI=3.1415926J=M!网格数dL=2.0!计算区域TT=0.4!总时间Sf=0.8!时间步长因子callInit(U,dx)T=01dt=CFL(U,

2、dx)T=T+dtwrite(*,*)T=,T,dt=,dtcallMacCormack_1D_Solver(U,Uf,Ef,dx,dt)callbound(U,dx)if(T.lt.TT)goto1callOutput(U,dx)end!计算时间步长!入口：U,当前物理量，dx,网格宽度；!返回:时间步长。!doubleprecisionfunctionCFL(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)dmaxvel=1e-10do10i=1,Juu

3、=U(i,1)/U(i,0)p=(GAMMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)vel=dsqrt(GAMMA*p/U(i,0)+dabs(uu)if(vel.gt.dmaxvel)dmaxvel=vel10continueCFL=Sf*dx/dmaxvelend!初始化!入口:无；!出口：U,已经给定的初始值，dx，网格宽度。!subroutineInit(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)!初始条件rou1=1.0u1

4、=0v1=0p1=1.0rou2=0.125u2=0v2=0p2=0.1dx=dL/Jdo20i=0,J/2U(i,0)=rou1U(i,1)=rou1*u1U(i,2)=p1/(GAMMA-1)+0.5*rou1*u1*u120continuedo21i=J/2+1,J+1U(i,0)=rou2U(i,1)=rou2*u2U(i,2)=p2/(GAMMA-1)+0.5*rou2*u2*u221continueend!边界条件!入口：dx,网格宽度；!出口:U，已经给定边界。!subroutinebound(U,dx)implicitdoubleprecision(a-h,o-z)common

5、/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)!左边界do30k=0,2U(0,k)=U(1,k)30continue!右边界do31k=0,2U(J+1,k)=U(J,k)31continueend!根据U计算E!入口：U,当前U矢量；!出口：E,计算得到的E矢量，!U、E定义见Euler方程组。!subroutineU2E(U,E,is,in)implicitdoubleprecision(a-h,o-z)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2),E(0:J+1,0

6、:2)do40i=is,inuu=U(i,1)/U(i,0)p=(GAMMA-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+pE(i,2)=(U(i,2)+p)*uu40continueend!一维差分格式求解器!入口：U，上一时刻U矢量，!Uf、Ef,临时变量，!dx,网格宽度，dt,，时间步长；!出口：U,计算得到得当前时刻U矢量。!subroutineMacCormack_1D_Solver(U,Uf,Ef,dx,dt)implicitdoubleprecision(a-h,o-z)common/G

7、_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2),Uf(0:J+1,0:2)dimensionEf(0:J+1,0:2)r=dt/dxdnu=0.25do60i=1,Jdo60k=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 continuedo61k=0,2do

8、61i=1,JU(i,k)=Ef(i,k)61 continuecallU2E(U,Ef,0,J+1)do63i=0,Jdo63k=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)callU2E(Uf,Ef,0,J)do64i=1,Jdo64k=0,2!U(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 continueend!输出结果,用数据格式画图!入口：U,当前时刻U矢量，!dx,网格宽度；!出口:无。

9、!subroutineOutput(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)open(1,file=MacCormack_1Dresult.txt,status=unknown)do80i=0,J+1rou=U(i,0)uu=U(i,1)/roup=(GAMMA-1)*(U(i,2)-0.5*U(i,0)*uu*uu)write(1,81)i*dx,rou,uu,p,U(i,2)80continueclose(1)81format(D20.10,D20.10,D20.10,D20.10,D20.10)end4