求助：SIMPLE 方法解Lid-driven cavity flow问题

luckyxu

各位大侠：
我是CFD领域的新兵。现在想用SIMPLE方法解二维Lid-driven cavity flow问题。用有限差分法。程序是在sunmoon大侠二维库塔流的程序基础上改的。流程是按照Anderson: Computational Fluid Dynamics, The Basics with Applications中第九章的第四节中的流程。但是怎么也算不出正确的速度场，为此寝食不安。希望大家能够帮忙看一下程序，找出问题症结所在，本人不胜感激。众人拾柴火焰高，何况本论坛高手如林，请不吝赐教。以下是程序：
PROGRAM driven_cavity
    IMPLICIT NONE
    REAL,DIMENSION(51,51) :: stream,p,px,x(51),y(51),d0,&
    pp,ppm,ppb,u(52,51),v(53,52),ux(52,51),vx(53,52),rux(52,51),rvx(53,52)
    INTEGER :: is_beg,is_end,js_beg,js_end,nstep,cp,i,j,nstep1,nstep2
    REAL :: delx,dely,delt,rho,mu,& a0,b0,c0,fct1,fct2,fct3,delx2,dely2,vb,vbb,&     ruu53,ruu33,ruu43,ruu42,uu1,uu2,ax,bx,ub,ubb,rvv45,rvv43,rvv54,rvv34,vv1,vv2,&
ruv44,ruv42,rvu54,rvu34,height,width,Error0,Error1,Error2,Error00,Error3,&
Error_stream,temp
    DATA is_beg,is_end,js_beg,js_end,rho,mu,delt,fct1,fct2,fct3,height,width&
/1,51,1,51,998.2,1.005e-3,0.0001,0.8,0.8,0.5,0.1,0.1/
        delx=width/REAL(is_end-1)
        dely=height/REAL(js_end-1)
delx2=delx**2
dely2=dely**2
!assign initial values
                ux=0.0
                u=0.0
                u(is_beg:is_end+1,js_end)=2.0e-2
                rux=rho*ux
     vx=0.0
                v=0.0
     rvx=rho*vx
                px=0.0
                p=0.0
                pp=0.0
                nstep=0
     ppm=0.0
     stream=0.0
Error1=1.
DO WHILE (Error1>9.5e-7)
            nstep=nstep+1
!
!     Pressure correction (SIMPLE) algorithm
!     step 1. solve ru* and rv*
!     a. for ru and u, pp.439(9.58)
         DO i=is_beg,is_end-1
            DO j=js_beg+1,js_end-1
                vb=0.5*(v(i+1,j+1)+v(i+2,j+1))
                vbb=0.5*(v(i+1,j)+v(i+2,j))
                ruu53=rho*u(i+2,j)*u(i+2,j)
                ruu33=rho*u(i,j)*u(i,j)
                ruv44=rho*u(i+1,j+1)*vb
                ruv42=rho*u(i+1,j-1)*vbb
                uu1=(u(i+2,j)-2.*u(i+1,j)+u(i,j))/delx2
                uu2=(u(i+1,j+1)-2.*u(i+1,j)+u(i+1,j-1))/dely2
                ax=-((ruu53-ruu33)/(2.*delx)+(ruv44-ruv42)/(2.*dely))+mu*(uu1+uu2)
                rux(i+1,j)=rux(i+1,j)+ax*delt-delt/delx*(px(i+1,j)-px(i,j))
                ux(i+1,j)=rux(i+1,j)/rho
            END DO
         END DO
!Boudary values
              ux(is_beg,js_beg+1:js_end-1)=0.0
              ux(is_end+1,js_beg+1:js_end-1)=0.0
         ux(is_beg:is_end+1,js_beg)=0.0
              ux(is_beg:is_end+1,js_end)=2.0e-2
!
!b. for rv and v, pp.440(9.59)
            DO i=is_beg,is_end
              DO j=js_beg,js_end-1
                ub=0.5*(u(i+1,j)+u(i+1,j+1))
                ubb=0.5*(u(i,j)+u(i,j+1))
                rvv45=rho*v(i+1,j+2)*v(i+1,j+2)
                rvv43=rho*v(i+1,j)*v(i+1,j)
                rvu54=rho*v(i+2,j+1)*ub
                rvu34=rho*v(i,j+1)*ubb
                vv1=(v(i+2,j+1)-2.*v(i+1,j+1)+v(i,j+1))/delx2
                vv2=(v(i+1,j+2)-2.*v(i+1,j+1)+v(i+1,j))/dely2
                bx=-((rvv45-rvv43)/(2.*dely)+(rvu54-rvu34)/(2.*delx))+mu*(vv1+vv2)
                rvx(i+1,j+1)=rvx(i+1,j+1)+bx*delt-delt/dely*(px(i,j+1)-px(i,j))
                vx(i+1,j+1)=rvx(i+1,j+1)/rho
              END DO
            END DO
!Boundary values
              vx(is_end+2,js_beg+1:js_end)=0.0
              vx(is_beg,js_beg+1:js_end)=0.0
              vx(is_beg:is_end+2,js_beg)=0.0
              vx(is_beg:is_end+2,js_end+1)=0.0

!
!     step 2. solve p'; from pressure correction formula
!             by relaxation method (inner iteration)
!
!calculate pressure corrector,  pp.441(9.60)
          a0=2.0*(delt/delx2+delt/dely2)
          b0=-delt/delx2
          c0=-delt/dely2
          DO i=is_beg+1,is_end-1
         DO j=js_beg+1,js_end-1
                pp(i,j)=0.0
                d0(i,j)=(rux(i+1,j)-rux(i,j))/delx+(rvx(i+1,j+1)-rvx(i+1,j))/dely
             END DO
        END DO
! calculate error of d0
Error1=0.0
DO i=is_beg+1,is_end-1
                DO j=js_beg+1,js_end-1
                  Error0=ABS(ux(i,j)-u(i,j))!/ABS(u(i,j))
                  Error1=MAX(Error1,Error0)
      END DO
           END DO
!calculate p'; using relaxation method
           Error2=1.0
   nstep1=0
DO WHILE (Error2>1.0e-7)
   nstep1=nstep1+1
   Error2=0.0
              DO i=is_beg+1,is_end-1
                DO j=js_beg+1,js_end-1
  ppm(i,j)=-1./a0*(b0*pp(i+1,j)+b0*ppm(i-1,j)+c0*pp(i,j+1)+&
                           c0*ppm(i,j-1)+d0(i,j))
                  !ppm(i,j)=pp(i,j)+fct1*(ppb(i,j)-pp(i,j))
                  Error00=ABS(ppm(i,j)-pp(i,j))/ABS(pp(i,j))
  Error2=MAX(Error2,Error00)
  pp(i,j)=ppm(i,j)
                END DO
              END DO
           END DO
!
!
!     step 3. calculate p at n+1 step
!   
DO i=is_beg+1,is_end-1
         DO j=js_beg+1,js_end-1
p(i,j)=px(i,j)+fct2*ppm(i,j)
           END DO
          END DO
!update boundary values
           p(is_beg:is_end,js_beg)=p(is_beg:is_end,js_beg+1)
           p(is_beg:is_end,js_end)=p(is_beg:is_end,js_end-1)
           p(is_beg,js_beg+1:js_end-1)=p(is_beg+1,js_beg+1:js_end-1)
           p(is_end,js_beg+1:js_end-1)=p(is_end-1,js_beg+1:js_end-1)
!     step 4. calculate u and v at n+1 step,from pp.239(6.100,6,101)
!
             DO i=is_beg,is_end-1
                DO j=js_beg,js_end
                 u(i+1,j)=ux(i+1,j)-delt/delx/rho*(pp(i+1,j)-pp(i,j))
                END DO
             END DO
DO i=is_beg,is_end
                DO j=js_beg,js_end-1
              v(i+1,j+1)=vx(i+1,j+1)-delt/dely/rho*(pp(i,j+1)-pp(i,j))
                END DO
             END DO
!
!     update boundary conditions
!
           u(is_beg:is_end+1,js_end)=2.0e-2
            u(is_beg:is_end+1,js_beg)=0.
            v(is_beg:is_end+2,js_end+1)=0.
            v(is_beg:is_end+2,js_beg)=0.
!
            u(is_beg,js_beg:js_end)=0.0
            u(is_end+1,js_beg:js_end)=0.0
            v(is_beg,js_beg:js_end+1)=0.0
            v(is_end+2,js_beg:js_end+1)=0.0

!     step 5. calculate stream function
      
      Error3=1.0
  nstep2=0
DO WHILE (Error3>1.0e-2)
      nstep2=nstep2+1
  Error3=0.0
  DO i=is_beg+1,is_end
   DO j=js_beg+1,js_end
    temp=stream(i,j)
stream(i,j)=0.5*(stream(i-1,j)+stream(i,j-1))+0.125*dely*&
            (u(i,j)+u(i+1,j)+u(i,j-1)+u(i+1,j-1))-0.125*delx*&(v(i,j)+v(i+1,j)+v(i,j+1)+v(i+1,j+1))
Error_stream=ABS(stream(i,j)-temp)/ABS(temp)
Error3=MAX(Error_stream,Error3)
   END DO
  END DO
  stream(is_beg,js_beg:js_end)=stream(is_beg+1,js_beg:js_end)
      stream(is_beg:is_end,js_beg)=stream(is_beg:is_end,js_beg+1)
     END DO

!
!     step 6. update pressure field
            px=p
      WRITE (*,*) nstep,Error1,nstep1,Error2,nstep2,Error3
      END DO
       OPEN(25,file="C:\Research\LSM\data.dat")
   WRITE (25,*) ';VARIABLES="X","Y",""';
       WRITE (25,*) ';ZONE I=51, J=51 , F=POINT';
       DO i=is_beg,is_end
       DO j=js_beg,js_end
       x(i)=REAL(i)*delx
       y(j)=REAL(j)*dely
       WRITE(25,"(2f10.3,f36.30)") x(i),y(j),p(i,j)
       END DO
       END DO
   OPEN(26,file="C:\Research\LSM\data_u1.dat")
       DO i=is_beg,is_end
       DO j=js_beg,js_end
       x(i)=REAL(i)*delx
       y(j)=REAL(j)*dely
       WRITE(26,"(2f10.3,2f36.30)") x(i),y(j),u(i,j),v(i,j)
       END DO
       END DO
       OPEN(27,file="C:\Research\LSM\stream.dat")
       WRITE (27,*) ';VARIABLES="X","Y","STREAM"';
       WRITE (27,*) ';ZONE I=51, J=51 , F=POINT';
   DO i=is_beg,is_end
       DO j=js_beg,js_end
       x(i)=REAL(i)*delx
       y(j)=REAL(j)*dely
       WRITE(27,"(2f10.3,f36.30)") x(i),y(j),stream(i,j)
   END DO
   END DO
      END PROGRAM driven_cavity
注：
vb: v_bar
vbb:v_bar_bar
ux:u*
vx:v*
px:p*
pp:p';