利用MacCormack兩部差分格式求解一維激波管問題Fortran程序

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!氣體常數(shù)PI=3.1415926J=M!網(wǎng)格數(shù)dL=2.0!計算區(qū)域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,網(wǎng)格寬度；!返回:時間步長。!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,Ju

3、u=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,已經(jīng)給定的初始值，dx，網(wǎng)格寬度。!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.0u

4、1=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,網(wǎng)格寬度；!出口:U，已經(jīng)給定邊界。!subroutinebound(U,dx)implicitdoubleprecision(a-h,o-z)commo

5、n/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!根據(jù)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,

6、0: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,網(wǎng)格寬度，dt,，時間步長；!出口：U,計算得到得當前時刻U矢量。!subroutineMacCormack_1D_Solver(U,Uf,Ef,dx,dt)implicitdoubleprecision(a-h,o-z)common/

7、G_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!開關(guān)函數(shù)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)continuedo61k=0,2do61

8、i=1,JU(i,k)=Ef(i,k)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)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)continueend!輸出結(jié)果,用數(shù)據(jù)格式畫圖!入口：U,當前時刻U矢量，!dx,網(wǎng)格寬度；!出口:無。!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)/r