资源描述:
《利用maccormack两部差分格式求解一维激波管问题fortran程序》由会员上传分享,免费在线阅读,更多相关内容在行业资料-天天文库。
1、!MacCormack1D.for!------------------------------------------------------------------------------------------------------------!利用MacCormack两部差分格式求解一维激波管问题!!-----------------------------------------------------------------------------------------------------------programMacCormack1Dimpli
2、citdoubleprecision(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,dx)T=T+dtwrite(*,*)'T=',T,'dt=',dtcallMacCormack_1D_Solve
3、r(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_
4、def/GAMMA,PI,J,JJ,dL,TT,SfdimensionU(0:J+1,0:2)5dmaxvel=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))+dabs(uu)if(vel.gt.dmaxvel)dmaxvel=vel10continueCFL=Sf*dx/dmaxvelend!-------------------------------------------------------!初始化!入口:无;!出口:U
5、,已经给定的初始值,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=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/(
6、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*u221continueend5!-------------------------------------------------------!边界条件!入口:dx,网格宽度;!出口:U,已经给定边界。!-------------------------------------------------------subroutinebound(U,dx)imp
7、licitdoubleprecision(a-h,o-z)common/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方程组。!---------