基于steger warming 矢通量分裂 的 二双马赫数反射模拟程序请教

lvp490

大家好，
我写了一段基于steger warming 矢通量分裂 的 二维双马赫数反射模拟程序。但是程序运行只能迭代2步，在这里请教各位专家，能否帮我看看程序哪里不对。双马赫数反射是计算流体力学验证算法的基本例子， 大家帮我看看哪里不对？自己是电气工程专业，自学了一阵CFD。再次谢谢大家。得到指点后，本人会有感谢。

        PROGRAM main

      IMPLICIT NONE

      Integer i,j,iter

      INTEGER NTMAX,IIX, IIY, MT,  F_SAVE,MN
      Parameter (NTMAX=1000000,IIX=960,IIY=240,MT=2)
      PARAMETER(F_SAVE=100)

      Real*8 gamma, tout, cfl
      Parameter(gamma=1.4,tout=0.2,cfl=0.45)

      REAL*8 xl,yl
      common /group1/ xl,yl
      Real*8 dx,dy
      common /group3/ dx,dy

      REAL*8 X(0:IIX), Y(0:IIY)
      common /group2/X,Y

      Real*8 D(-MT:IIX+MT, -MT:IIY+MT), U(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 V(-MT:IIX+MT, -MT:IIY+MT), P(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 C(-MT:IIX+MT, -MT:IIY+MT)
     
      Real*8 U1(-MT:IIX+MT, -MT:IIY+MT),U2(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 U3(-MT:IIX+MT, -MT:IIY+MT),U4(-MT:IIX+MT, -MT:IIY+MT)

      Real*8 time, dt
     
      REAL*8 f1(-MT:IIX+MT, -MT:IIY+MT),f2(-MT:IIX+MT, -MT:IIY+MT)
      REAL*8 f3(-MT:IIX+MT, -MT:IIY+MT),f4(-MT:IIX+MT, -MT:IIY+MT)
      REAL*8 g1(-MT:IIX+MT, -MT:IIY+MT),g2(-MT:IIX+MT, -MT:IIY+MT)
      REAL*8 g3(-MT:IIX+MT, -MT:IIY+MT),g4(-MT:IIX+MT, -MT:IIY+MT)

      REAL*8 dsdt1(-MT:IIX+MT, -MT:IIY+MT),dsdt2(-MT:IIX+MT, -MT:IIY+MT)
      REAL*8 dsdt3(-MT:IIX+MT, -MT:IIY+MT),dsdt4(-MT:IIX+MT, -MT:IIY+MT)


      xl=4.0
      yl=1.0

      Call mesher
      Call init(D,U,V,P,C,U1,U2,U3,U4)

      time = 0.0
*
C.....输出初始时刻的数值流场
      CALL SAVE(0, D, U, V, C,P)
      
      WRITE(*,*)'*********************************************'
      WRITE(*,*)'   Time step N        TIME             '
      WRITE(*,*)'*********************************************'
*
      DO  iter = 1, NTMAX
*
C.....应用CFL条件计算时间步长
*
        CALL cflcon (U, V, C,time,dt)
*
        time = time + dt
      
     
       CALL BOUNDARY (D, U, V, P,time)
       

        Do j=0,IIY
            Do i=-1,IIX
C.....求解X方向的数值通量
*

          CALL swx(U1(i,j), U2(i,j), U3(i,j), U4(i,j),
     &             U1(i+1,j), U2(i+1,j), U3(i+1,j), U4(i+1,j),
     &             f1(i,j), f2(i,j), f3(i,j), f4(i,j))
            enddo
        enddo
C.....求解Y方向的数值通量
       Do i=0,IIX
            Do j=-1,IIY
*
           CALL swy(U1(i,j),U2(i,j),U3(i,j),U4(i,j),
     &    U1(i,j+1), U2(i,j+1), U3(i,j+1), U4(i,j+1),
     &     g1(i,j), g2(i,j), g3(i,j), g4(i,j))

            enddo
        enddo


        DO i=0,IIX
            Do j=0,IIY
        
        dsdt1(i,j) = -(f1(i,j)-f1(i-1,j))/dx - (g1(i,j)-g1(i,j-1))/dy
        dsdt2(i,j) = -(f2(i,j)-f1(i-1,j))/dx - (g2(i,j)-g1(i,j-1))/dy
        dsdt3(i,j) = -(f3(i,j)-f1(i-1,j))/dx - (g3(i,j)-g1(i,j-1))/dy
        dsdt4(i,j) = -(f4(i,j)-f1(i-1,j))/dx - (g4(i,j)-g1(i,j-1))/dy

         U1(i,j) = U1(i,j) + dt*dsdt1(i,j)
         U2(i,j) = U2(i,j) + dt*dsdt2(i,j)
         U3(i,j) = U3(i,j) + dt*dsdt3(i,j)
         U4(i,j) = U4(i,j) + dt*dsdt4(i,j)
         D(i,j) = U1(i,j)
         U(i,j) = U2(i,j)/U1(i,j)
         V(i,j) = U3(i,j)/U1(i,j)
         P(i,j)=(gamma-1)*(U4(i,j)-0.5*(U2(i,j)**2+U3(i,j)**2)*U1(i,j))
         C(i,j) = SQRT(gamma*P(i,j)/D(i,j))
        
            Enddo
       Enddo
      

        WRITE(*,*) iter, time

*C.....在指定时间步上输出文件
*
           IF(MOD(iter, F_SAVE).EQ.0) THEN
             CALL SAVE(iter/F_SAVE, D, U, V, C,P)
           ENDIF
C
        IF(ABS(TIME - TOUT).LE.1.0E-06) THEN
*
C.....输出最终结果
*
           CALL SAVE(iter/F_SAVE+1, D, U, V, C,P)
*
           WRITE(*,*)'*********************************************'
           WRITE(*,*)'   Number of time steps = ',iter
           WRITE(*,*)'*********************************************'
*
           EXIT
*
        ENDIF

*
        Enddo

        Stop
      END Program
*----------------------------------------------------------------------------------------

       SUBROUTINE mesher
*
C.....该子程序的作用：生成结构网格
*
        IMPLICIT NONE

      INTEGER  IIX, IIY, MT,i,j
      Parameter (IIX=960,IIY=240,MT=2)
       
      REAL*8 X(0:IIX), Y(0:IIY)
      common /group2/X,Y

      REAL*8 xl,yl
      common /group1/ xl,yl
      Real*8 dx,dy
      common /group3/ dx,dy
      
      dx=xl/IIX
      dy=yl/IIY
       
        DO i =0,IIX
          X(i) = i*dx
        ENDDO

        DO j = 0, IIY
          Y(j) = j*dy
        ENDDO

        END
*----------------------------------------------------------------------------------------
*
      SUBROUTINE init(D,U,V,P,C,U1,U2,U3,U4)
C.....该子程序的作用：给计算区域赋初场
*
        IMPLICIT NONE

      INTEGER IIX, IIY, MT,i,j
      Parameter (IIX=960,IIY=240,MT=2)
      Real gamma
      Parameter(gamma=1.4)

        REAL*8 X(0:IIX), Y(0:IIY)
      common /group2/X,Y
   
      
      Real*8 D(-MT:IIX+MT, -MT:IIY+MT), U(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 V(-MT:IIX+MT, -MT:IIY+MT), P(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 C(-MT:IIX+MT, -MT:IIY+MT)
     
      Real*8 U1(-MT:IIX+MT, -MT:IIY+MT),U2(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 U3(-MT:IIX+MT, -MT:IIY+MT),U4(-MT:IIX+MT, -MT:IIY+MT)
      


        DO i = 0, IIX
          DO j = 0, IIY
C.....赋初值

            IF(Y(j).LT.(X(i)-1./6.)*SQRT(3.0)) THEN
              D(i,j) = gamma
              U(i,j) = 0.0
              V(i,j) = 0.0
              P(i,j) = 1.0
            ELSE
            D(i,j) = 8.0
              U(i,j) = 7.14471
              V(i,j) =-4.125
              P(i,j) = 116.5
            ENDIF

            C(i,j) = SQRT(gamma*P(i,j)/D(i,j))
          U1(i,j)=D(i,j)
          U2(i,j)=D(i,j)*U(i,j)
          U3(i,j)=D(i,j)*V(i,j)
          U4(i,j)=.5*D(i,j)*(U(i,j)**2+V(i,j)**2)+P(i,j)/(gamma-1)
C
          ENDDO
        ENDDO
C
        END
*----------------------------------------------------------------------------------------
     
      SUBROUTINE cflcon (U, V, C,time,dt)
*
C.....该子程序的作用：计算时间步长
*
      IMPLICIT NONE

      INTEGER IIX, IIY, MT,i,j
      Parameter (IIX=960,IIY=240,MT=2)
     
      Real*8 gamma, tout, cfl
      Parameter(gamma=1.4,tout=0.2,cfl=0.45)

      REAL*8  DTL,SPX, SPY
   
     
      Real*8 dx,dy
      common /group3/ dx,dy

     
      Real*8 U(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 V(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 C(-MT:IIX+MT, -MT:IIY+MT)
      

      Real*8 time, dt
     
      
      dt = 1.0E+10
     
      DO i = 0, IIX
        DO j = 0, IIY
*
C.....计算X和Y方向上的特征速度
*
          SPX = C(i,j) + ABS(U(i,j))
          SPY = C(i,j) + ABS(V(i,j))
C
          DTL = MIN(dx/SPX, dy/SPY)
C
          IF(DTL.LT.dt)dt = DTL
           
         
        ENDDO
      ENDDO
C
      dt = cfl*dt

       write(*,*) dt

C.....防止TIME大于TOUT
*
      IF((time + dt).GT.tout)THEN
         dt = tout - time
      ENDIF
C
      END
*----------------------------------------------------------------------------------------
*
      SUBROUTINE BOUNDARY (D, U, V, P,time)
*
C.....该子程序的作用：定义边界条件
*
      IMPLICIT NONE

      INTEGER IIX, IIY, MT,i,j
      Parameter (IIX=960,IIY=240,MT=2)

      Real*8 gamma
      Parameter(gamma=1.4)

      REAL*8 xl,yl
      common /group1/ xl,yl
      Real*8 dx,dy
      common /group3/ dx,dy

      REAL*8 X(0:IIX), Y(0:IIY)
      common /group2/X,Y
   
      
      Real*8 D(-MT:IIX+MT, -MT:IIY+MT), U(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 V(-MT:IIX+MT, -MT:IIY+MT), P(-MT:IIX+MT, -MT:IIY+MT)
      
     

      Real*8 time
     


      REAL SS, XIN, YIN   

C....设置左、右边界条件
*
        i = 1
        DO j = 0, IIY
C
            D(-i,j) = D(i-1,j)
            U(-i,j) = U(i-1,j)
            V(-i,j) = V(i-1,j)
            P(-i,j) = P(i-1,j)
C
                D(IIX+i,j) = D(IIX-i+1,j)
            U(IIX+i,j) = U(IIX-i+1,j)
            V(IIX+i,j) = V(IIX-i+1,j)
            P(IIX+i,j) = P(IIX-i+1,j)
C
          ENDDO

*
C.....设置上边界条件
*
        j = 1
        DO i = 0, IIX
C
            SS = 10.0
            XIN = X(i) + 0.5*dx
            YIN = Y(IIY) + 0.5*dy
C
                IF(YIN-SQRT(3.)*(XIN-1./6.).GT.-SS*time/0.5) THEN
              D(i,IIY+j) = 8.0
              U(i,IIY+j) = 7.14471
              V(i,IIY+j) =-4.125
              P(i,IIY+j) = 116.5
            ELSE
              D(i,IIY+j) = gamma
              U(i,IIY+j) = 0.0
              V(i,IIY+j) = 0.0
              P(i,IIY+j) = 1.0
            ENDIF
C
          ENDDO
       
*
C.....设置下边界条件
*
        j = 1
        DO i = 0, IIX
C
            XIN = X(i) + 0.5*dx
            YIN = Y(0) + 0.5*dy
C
                IF(XIN.LT.1./6.) THEN
              D(i,-j) = 8.0
              U(i,-j) = 7.14471
              V(i,-j) =-4.125
              P(i,-j) = 116.5
            ELSE
              D(i,-j) = D(i,j-1)
              U(i,-j) = U(i,j-1)
              V(i,-j) =-V(i,j-1)
              P(i,-j) = P(i,j-1)
            ENDIF
C
          ENDDO
       
C      
      END
*----------------------------------------------------------------------------------------
      
c...Steger-Warming flux vector splitting along X-direction
      
      subroutine swx(rp,rup,rvp,retp,rm,rum,rvm,retm,f1,f2,f3,f4)

      Real*8 gamma
      Parameter(gamma=1.4)
      
      Real*8 rp,rup,rvp,retp,rm,rum,rvm,retm
      Real*8 rhop,rhom,up,um,pp,pm,vp,vm
      Real*8 lambdaf1p, lambdaf2p, lambdaf3p, lambdaf4p
      Real*8 lambdaf1m, lambdaf2m, lambdaf3m, lambdaf4m
      Real*8 fp1, fp2, fp3, fp4,fm1, fm2, fm3,fm4
      Real*8 f1, f2, f3, f4
     

c...Convert conservative variables to primitive variables.

      rhop = rp
      rhom = rm
      up =   rup/rp
      um =   rum/rm
      vp =   rvp/rp
      vm =   rvm/rm
     
      pp = (gamma-1.)*(retp -
     & .5*(rup*rup+rvp*rvp)/rp)

       pm = (gamma-1.)*(retm -
     & .5*(rum*rum+rvm*rvm)/rm)
      
      cp = sqrt(gamma*pp/rp)
      cm = sqrt(gamma*pm/rm)
      

c...Compute wave speed splitting along X-dirction

      lambdaf1p =0.5*(up+abs(up))
      lambdaf2p=0.5*(up+abs(up))
      lambdaf3p =0.5*((up-cp)+abs(up-cp))
      lambdaf4p =0.5*((up+cp)+abs(up+cp))

      lambdaf1m = 0.5*(um-abs(um))
      lambdaf2m = 0.5*(um-abs(um))
      lambdaf3m = 0.5*((um-cm)-abs(um-cm))
      lambdaf4m = 0.5*((um+cm)-abs(um+cm))

c...Compute flux splitting along X-direction
      fp1 = .5*rhop*(2.*(gamma-1.)*lambdaf1p+lambdaf3p+lambdaf4p)
     &         /gamma
      fp2 = .5*rhop*(2.*(gamma-1.)*up*lambdaf1p+(up
     & -cp)*lambdaf3p+ (up+cp)*lambdaf4p)/gamma
      fp3 = .5*rhop*(2.*(gamma-1.)*vp*lambdaf1p+vp
     & *lambdaf3p+ vp*lambdaf4p)/gamma
      fp4 = .5*rhop*((gamma-1.)*lambdaf1p*(up*up+vp
     & *vp)+ .5*((up-cp)**2+vp**2)*lambdaf3p+
     &  . 5*((up+cp)**2+vp**2)*lambdaf4p+.5*(3-gamma)*
     & (lamdaf3p+lamdaf4p)*cp*cp/(gamma-1))/gamma


      fm1 = .5*rhom*(2.*(gamma-1.)*lambdaf1m+lambdaf3m+lambdaf4m)
     &        /gamma
      fm2 = .5*rhom*(2.*(gamma-1.)*um*lambdaf1m+(um
     & -cm)*lambdaf3m + (um+cm)*lambdaf4m)/gamma
      fm3 = .5*rhom*(2.*(gamma-1.)*vm*lambdaf1m+vm
     & *lambdaf3m + vm*lambdaf4m)/gamma
      fm4 = .5*rhom*((gamma-1.)*lambdaf1m*(um*um
     & +vm*vm)+ .5*((um-cm)**2+vm**2)*lambdaf3m+.5*((um+cm)**2+vm**2)*
     & lambdaf4m +.5*(3-gamma) *(lamdaf3m+lamdaf4m)*cm*cm/(gamma-1))
     & /gamma

c...Compute conserative numerical fluxes along X-direction
      f1 = fp1 + fm1
      f2 = fp2 + fm2
      f3 = fp3 + fm3
      f4 = fp4 + fm4

      return
      End

*----------------------------------------------------------------------------------------

c...Steger-Warming flux vector splitting along Y-direction

      subroutine swy(rp,rup,rvp,retp,rm,rum,rvm,retm,g1,g2,g3,g4)

      real*8 gamma
      parameter(gamma=1.4)
      real*8 rp,rup,rvp,retp,rm,rum,rvm,retm
      real*8 rhop,rhom,up,um,pp,pm,vp,vm
      real*8 lambdag1p, lambdag2p, lambdag3p, lambdag4p
      real*8 lambdag1m, lambdag2m, lambdag3m, lambdag4m
      real*8 gp1, gp2, gp3, gp4,gm1, gm2, gm3,gm4
      real*8 g1, g2, g3, g4
     

c...Convert conservative variables to primitive variables.

      rhop = rp
      rhom = rm
      up =   rup/rp
      um =   rum/rm
      vp =   rvp/rp
      vm =   rvm/rm
     
      pp = (gamma-1.)*(retp -
     & .5*(rup*rup+rvp*rvp)/rp)

       pm = (gamma-1.)*(retm -
     & .5*(rum*rum+rvm*rvm)/rm)
      
      cp = sqrt(gamma*pp/rp)
      cm = sqrt(gamma*pm/rm)



c...Compute wave speed splitting along Y-dirction
      lambdag1p =0.5*(vp+abs(vp))
      lambdag2p =0.5*(vp+abs(vp))
      lambdag3p = 0.5*((vp-cp)+abs(vp-cp))
      lambdag4p = 0.5*((vp+cp)+abs(vp+cp))
   
      lambdag1m = 0.5*(vm-abs(vm))
      lambdag2m = 0.5*(vm-abs(vm))
      lambdag3m = 0.5*((vm-cm)-abs(vm-cm))
      lambdag4m =0.5*((vm+cm)-abs(vm+cm))


c...Compute flux splitting along Y-direction



      gp1 = .5*rhop*(2.*(gamma-1.)*lambdag1p+lambdag3p+lambdag4p)
     &        /gamma
      gp2 = .5*rhop*(2.*(gamma-1.)*up*lambdag1p+up
     & *lambdag3p+ up*lambdag4p)/gamma
      gp3 = .5*rhop*(2.*(gamma-1.)*vp*lambdag1p+
     & (vp-cp)*lambdag3p+ (vp+cp)*lambdaf4p)/gamma
      gp4 = .5*rhop*((gamma-1.)*lambdag1p*(up*up+vp
     & *vp)+ .5*((vp-cp)**2+up**2)*lambdag3p+
     & .5*((vp+cp)**2+up**2)*lambdag4p+
     & .5*(3-gamma)*(lamdag3p+lamdag4p)*cp*cp/(gamma-1))/gamma


      gm1 = .5*rhom*(2.*(gamma-1.)*lambdag1m+lambdag3m+lambdag4m)
     &      /gamma
      gm2 = .5*rhom*(2.*(gamma-1.)*um*lambdag1m+um
     & *lambdaf3m+ um*lambdaf4m)/gamma
      gm3 = .5*rhom*(2.*(gamma-1.)*vm*lambdag1m+
     & (vm-cm)*lambdag3m+ (vm+cm)*lambdag4m)
     &  /gamma
      gm4 = .5*rhom*((gamma-1.)*lambdag1m*(um*um+
     & vm*vm)+ .5*((vm-cm)**2+um**2)
     & *lambdag3m+.5*((vm+cm)**2+um**2)*lambdag4m
     & +.5*(3-gamma)*(lamdag3m+lamdag4m)*cm*cm/(gamma-1))
     & /gamma

c...Compute conserative numerical fluxes along Y-direction
      g1 = gp1 + gm1
      g2 = gp2 + gm2
      g3 = gp3 + gm3
      g4 = gp4 + gm4
      
      return
      End
*----------------------------------------------------------------------------------------
      SUBROUTINE SAVE(INDEX_FILE, D, U, V, C, P)
*
C.....该子程序的作用：输出计算结果
*
      IMPLICIT NONE
      
      INTEGER IIX, IIY, MT,i,j
      Parameter (IIX=960,IIY=240,MT=2)

      INTEGER IX, JY, INDEX_FILE

      Real*8 D(-MT:IIX+MT, -MT:IIY+MT), U(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 V(-MT:IIX+MT, -MT:IIY+MT), P(-MT:IIX+MT, -MT:IIY+MT)
      Real*8 C(-MT:IIX+MT, -MT:IIY+MT)
  
      REAL*8 X(0:IIX), Y(0:IIY)
      common /group2/X,Y
*
        CHARACTER*4  INDEX_NAME
        CHARACTER*8 TITLED,TITLEU,TITLEV,TITLEC,TITLEP
       
        WRITE(INDEX_NAME,'(I4.4)') INDEX_FILE
        TITLED="D"//INDEX_NAME
      TITLEU="U"//INDEX_NAME
      TITLEV="V"//INDEX_NAME
      TITLEC="C"//INDEX_NAME
      TITLEP="P"//INDEX_NAME
        OPEN(UNIT=3,FILE=TITLED//'.dat')
      OPEN(UNIT=4,FILE=TITLEU//'.dat')
      OPEN(UNIT=5,FILE=TITLEV//'.dat')
      OPEN(UNIT=7,FILE=TITLEC//'.dat')
      OPEN(UNIT=6,FILE=TITLEP//'.dat')

*
        OPEN(8, FILE = 'x_1d.dat')
      OPEN(9, FILE = 'y_1d.dat')
*
C.....输出二维结果
*
      Do i=0,IIX
C
      WRITE(3,'(481(F12.6,1X))') (D(i,j),j=0,IIY)
      WRITE(4,'(481(F12.6,1X))') (U(i,j),j=0,IIY)
      WRITE(5,'(481(F12.6,1X))') (V(i,j),j=0,IIY)
      
      WRITE(7,'(481(F12.6,1X))') (C(i,j),j=0,IIY)
      WRITE(6,'(481(F12.6,1X))') (P(i,j),j=0,IIY)
C
        Enddo
     
C
      IX = IIX/2
      JY = IIY/2
*
C.....输出X方向的一维结果
*
      DO 10 i = 0, IIX
         WRITE(8,40) X(i),D(i,JY),U(i,JY),V(i,JY),P(i,JY)
10   CONTINUE       
*
C.....输出Y方向的一维结果
*
      DO 30 j = 0, IIY
         WRITE(9,40) Y(j),D(IX,j),U(IX,j),V(IX,j),P(IX,j)
30   CONTINUE
C
      CLOSE(8)
      CLOSE(9)

40   FORMAT(2(F12.6,2X),5X,3(F12.4,2X))
C
      END
*----------------------------------------------------------------------------------------