资源描述：
《利用MacCormack两部差分格式求解一维激波管问题Fortran程序.doc》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。
1、!MacCormack1D.for!------------------------------------------------------------------------------------------------------------!利用MacCormack两部差分格式求解一维激波管问题!!-----------------------------------------------------------------------------------------------
2、------------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.J=M!网格数dL=2.0!计算区域TT=0.4!总时间Sf=0.8!时间步长因子callInit(U,dx)T=01d
3、t=CFL(U,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,网格宽度；!返回:时间步长。!---------------------------
4、----------------------------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=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))+da
5、bs(uu)if(vel.gt.dmaxvel)dmaxvel=vel10continueCFL=Sf*dx/dmaxvelend!-------------------------------------------------------!初始化!入口:无；!出口:U,已经给定的初始值，dx，网格宽度。!-------------------------------------------------------subroutineInit(U,dx)implicitdoubleprecisi
6、on(a-h,o-z)common/G_def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)!初始条件rou1=1.0u1=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)
7、=p2/(GAMMA-1)+0.5*rou2*u2*u221continueend!-------------------------------------------------------!边界条件!入口:dx，网格宽度；!出口:U，已经给定边界。!-------------------------------------------------------subroutinebound(U,dx)implicitdoubleprecision(a-h,o-z)common/G_def/GA
8、MMA,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方程组。!------------------
