C ************************************************* C ** THE UNSTEADY VISCOUS FLOW COMPUTATION ** C ** ** C ** VORTICITY TRANSPORT EQUATION AND ** C ** POISSON EQUATION FOR STREAM FUNCTION ** C ** ** C ** FORWARD TIME CENTERED SPACE METHOD ** C ** ** C ** CARTESIAN COORDINATE ** C ** ** C ** FOR PC9801 VERSION ** C ** ** C ** 1999/12/17 fluid reseach ** C ************************************************* C C C ********************** C *** MAIN ROUTINE *** C ********************** C CALL INIT C CALL CLEAR C CALL CALFLW C CALL OUTPUT C C STOP END C C ************************ C ** TIME MARCHING LOOP ** C ************************ C SUBROUTINE CALFLW IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /TIMEC/ TIME COMMON /CONST1/ RE,DT C LOOP=0 C 50 LOOP=LOOP+1 C C CALL STREAM C CALL BOUNDV C CALL VORT C CALL JUDGE(ERR) C TIME=TIME+DT C IF(MOD(LOOP,100).EQ.0) THEN C WRITE(*,600) TIME,ERR 600 FORMAT(' TIME & RESIDUAL=',F10.4,2X,E9.4) C END IF C IF(TIME.GT.500.D0) GOTO 60 IF(ERR.LT.1.D-6) GOTO 90 C GOTO 50 C 60 WRITE(*,610) 610 FORMAT(' NOT CONVERGED!') C C 90 RETURN END C C *********************************** C ** INITIAL CONDITION SETTING ** C *********************************** C SUBROUTINE INIT C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /ITERAT/ ITER1 COMMON /CONV1/ OMEGA1,EPS1 COMMON /CONST1/ RE,DT COMMON /TIMEC/ TIME COMMON /DELTA/ DX,DY,DXI,DYI,DX2I,DY2I C C ***** NUMBER OF FIELD DIVIDE ***** C NX =21 NY =21 NXM1=NX-1 NYM1=NY-1 C C ***** DIVIDE LENGTH ***** C DX =1.0D0/DFLOAT(NXM1) DY =1.0D0/DFLOAT(NYM1) DXI =1.0D0/DX DYI =1.0D0/DY DX2I=1.0D0/DX/DX DY2I=1.0D0/DY/DY C C **** ITERATION CONSTANT FOR POISSON EQ.***** C ITER1 =100 OMEGA1=1.8D0 EPS1 =1.0D-6 C C *** REYNOLDS NUMBER *** C RE=100.0D0 C DT=0.001D0 C C RETURN END C C ********************************** C ** SUBROUTINE ARRAY INITIALIZED ** C ********************************** C SUBROUTINE CLEAR C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 C DO 100 J=1,NY DO 100 I=1,NX C PSI(I,J) =0.0D0 100 ZETA(I,J)=0.0D0 C C C READ(12,600) NX,NY,RE,TIME C 600 FORMAT(2I5,F10.1,F10.4) C C READ(12,610)((PSI(I,J),ZETA(I,J),I=1,NX),J=1,NY) C 610 FORMAT(2F15.7) C C RETURN END C C ****************************** C ** RESULTS OUTPUT TO FILE ** C ****************************** C SUBROUTINE OUTPUT C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /DELTA/ DX,DY,DXI,DYI,DX2I,DY2I COMMON /CONST1/ RE,DT COMMON /TIMEC/ TIME DIMENSION X(MX),UU(MY),VV(MX) C WRITE(3,600) NX,NY,RE,TIME 600 FORMAT(2I5,F10.1,F10.4) C WRITE(3,610)((PSI(I,J),ZETA(I,J),I=1,NX),J=1,NY) 610 FORMAT(2F15.7) C IH=(NX+1)/2 C DO 110 J=2,NYM1 C UU(J)= 0.5D0*(PSI(IH,J+1)-PSI(IH,J-1))*DYI 110 X(J)=DBLE(J-1)/DBLE(NY-1) C JH=(NY+1)/2 C DO 120 I=2,NXM1 C 120 VV(I)=-0.5D0*(PSI(I+1,JH)-PSI(I-1,JH))*DYI C X(NY)=1.D0 UU(NY)=1.D0 C WRITE(10,650) (X(J),UU(J),VV(J),J=1,NY) 650 FORMAT(3F10.5) C C RETURN END C C ******************************* C ** JUDGEMENT OF STEADY STATE ** C ******************************* C SUBROUTINE JUDGE(ERRS) C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /OLDZ/ ZO(MX,MY) C ERRS=0.D0 C DO 100 J=2,NYM1 DO 100 I=2,NXM1 C SQZ=(ZETA(I,J)-ZO(I,J))**2 100 ERRS=ERRS+SQZ C ERRS=DSQRT(ERRS/DFLOAT((NX-2)*(NY-2))) C C RETURN END C C ********************************** C ** VORTICITY BOUNDARY CONDITION ** C ********************************** C SUBROUTINE BOUNDV C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /DELTA/ DX,DY,DXI,DYI,DX2I,DY2I C C DO 100 I=1,NX C ZETA(I,1 )=-2.0D0* PSI(I,2 )*DY2I ZETA(I,NY)=-2.0D0*(PSI(I,NYM1)*DYI+1.D0)*DYI C 100 CONTINUE C DO 110 J=1,NY C ZETA(1 ,J)=-2.0D0*PSI(2 ,J)*DX2I ZETA(NX,J)=-2.0D0*PSI(NXM1,J)*DX2I C 110 CONTINUE C C RETURN END C C ******************************* C ** CALCULATE STREAM FUNCTION ** C ******************************* C SUBROUTINE STREAM C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /ITERAT/ ITER1 COMMON /CONV1/ OMEGA1,EPS1 COMMON /DELTA/ DX,DY,DXI,DYI,DX2I,DY2I C C ***** SUCCESSIVE OVER RELAXATION LOOP ***** C COEF=1.D0/(2.D0*(DX2I+DY2I)) C DO 300 LOOP=1,ITER1 C ERR=0.0D0 C C ***** CAL. INTERNAL POINT ***** C DO 100 J=2,NYM1 DO 100 I=2,NXM1 C PT=(PSI(I+1,J)+PSI(I-1,J))*DX2I & +(PSI(I,J+1)+PSI(I,J-1))*DY2I+ZETA(I,J) C DP=PT*COEF-PSI(I,J) C PSI(I,J)=PSI(I,J)+OMEGA1*DP C ERR=ERR+DP*DP C 100 CONTINUE C ERR=DSQRT(ERR/DFLOAT((NX-2)*(NY-2))) IF(ERR.LT.EPS1) GOTO 99 C 300 CONTINUE C 99 CONTINUE C CC WRITE(*,600) LOOP,ERR CC 600 FORMAT(' LOOP & ERROR=',I5,E12.5) C C RETURN END C C ************************* C ** CALCULATE VORTICITY ** C ************************* C SUBROUTINE VORT C PARAMETER (MX=257,MY=257) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON /FLOW/ PSI(MX,MY),ZETA(MX,MY) COMMON /DIVIDE/ NX,NY,NXM1,NYM1 COMMON /DELTA/ DX,DY,DXI,DYI,DX2I,DY2I COMMON /CONST1/ RE,DT COMMON /OLDZ/ ZO(MX,MY) DIMENSION UU(MX,MY),VV(MX,MY) C DO 100 J=1,NY DO 100 I=1,NX C 100 ZO(I,J)=ZETA(I,J) C DO 110 J=2,NYM1 DO 110 I=2,NXM1 C UU(I,J)= 0.5D0*(PSI(I,J+1)-PSI(I,J-1))*DYI 110 VV(I,J)=-0.5D0*(PSI(I+1,J)-PSI(I-1,J))*DXI C C DO 120 J=2,NYM1 DO 120 I=2,NXM1 C C ***** CONVECTIVE TERM ***** C CONV=(UU(I+1,J)*ZO(I+1,J)-UU(I-1,J)*ZO(I-1,J))*DXI*.5D0 & +(VV(I,J+1)*ZO(I,J+1)-VV(I,J-1)*ZO(I,J-1))*DYI*.5D0 C C ***** VISCOUS TERM ***** C VISC=(ZO(I+1,J)-2.D0*ZO(I,J)+ZO(I-1,J))*DX2I & +(ZO(I,J+1)-2.D0*ZO(I,J)+ZO(I,J-1))*DY2I C ZETA(I,J)=ZO(I,J)+DT*(-CONV+VISC/RE) C 120 CONTINUE C C RETURN END