C ***** File = cavity32.f ******************************************* C * Laminar Flow in a Square Cavity * C * * 2nd Order Central F.D.M. in Space * C * * 2nd Order Adams-Bashforth Scheme for Convection * C * * 2nd Order Crank-Nicholson Scheme for Diffusion * C * * Choice for Interpolation Method in Momentum Flux * C ****************************************************************** PROGRAM CAVITY32 PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*28 FLUX,CFLX(3) CHARACTER*20 FRMK,CFRM(2) CHARACTER*10 FILENAME COMMON /STEP/ISTEP,TSTEP /RE/RE /DT/DT /BT/BT /IADV/IADV COMMON /ITRU/ITRU /ERRU/ERRU /ITRP/ITRP /ERRP/ERRP /ETK/ETK COMMON /DIVMAX/DIVMAX /DXMIN/DXMIN /E/ENERGY /FILE/FILENAME COMMON /JALG/JALG /FLUX/FLUX /IFLX/IFLX COMMON /KATA/KATA /FRMK/FRMK DATA CFRM/'(1) DIVERGENCE FORM ','(2) GRADIENT FORM '/ DATA CFLX/'(1) SIMPLE AVERAGE ' &,'(2) LAGRANGIAN INTERPOLATION','(3) PROPER INTERPOLATION '/ WRITE(6,801) CFRM 801 FORMAT('FORM OF CONVECTION TERMS = ?',2(/3X,A20) ) READ(5,*) KATA FRMK=CFRM(KATA) WRITE(6,802) CFLX 802 FORMAT('INTERPLATION METHOD = ?',3(/3X,A28) ) READ(5,*) IFLX FLUX=CFLX(IFLX) WRITE(6,*) 'REYNOLDS NUMBER = ?' READ(5,*) RE WRITE(6,*) 'MESH DISTRIBUTION PARAMETER = ?' WRITE(6,*) ' JALG= 0 : UNIFORM MESH' WRITE(6,*) ' JALG=500-999 : TANH DISTRIBUTION' WRITE(6,*) ' JALG= 1000 : COS DISTRIBUTION' READ(5,*) JALG CALL SBMESH WRITE(6,900) DXMIN,RE*DXMIN**2/2.D0 900 FORMAT(/'LIMIT FOR TIME-INCREMENT DUE TO CONVECTION = ',F12.7 & /'LIMIT FOR TIME-INCREMENT DUE TO DIFFUSION = ',F12.7 & /'TIME-INCREMENT = ?') READ(5,*) DT BT=1.D0/DT WRITE(6,*) 'TIME FOR STOP = ?' READ(5,*) TSTOP TMIN=DSQRT(RE) WRITE(6,*) 'TIMESTEP FOR SKIP OUTPUT = ?' READ(5,*) ISKIP WRITE(6,*) 'NAME OF OUTPUT FILE = ?????????? (10 CHARACTERS)' READ(5,*) FILENAME CALL SBINIT OPEN(20,file=FILENAME//'.prn',status='new',form='formatted') WRITE(20,1100) WRITE(6,1100) ENERGB=0.D0 100 CONTINUE CALL SBWALL IF(KATA.EQ.1) THEN CALL SBFLUX ELSE CALL SBCONV ENDIF CALL SBPRED CALL SBSORU CALL SBRHSP CALL SBSORP CALL SBCORR ISTEP=ISTEP+1 TSTEP=DT*DBLE(ISTEP) IF(MOD(ISTEP,ISKIP).EQ.0) THEN CALL SBENGY ETK=(ENERGY-ENERGB)/(DBLE(ISKIP)*DT)/ENERGY CALL SBDIVU WRITE(20,1000) ISTEP,TSTEP,ENERGY,ETK,ITRU,ERRU,ITRP,ERRP,DIVMAX WRITE(6,1000) ISTEP,TSTEP,ENERGY,ETK,ITRU,ERRU,ITRP,ERRP,DIVMAX ENERGB=ENERGY ENDIF 1000 FORMAT(I7,F8.3,F14.11,E10.2,I5,E10.3,I5,2E10.3) 1100 FORMAT(3X,4HSTEP,7X,1HT,8X,6HENERGY,10H (DE/DT)/E &,5H ITRU,6X,4HERRU,5H ITRP,6X,4HERRP,4X,6HDIVMAX) IF(TSTEP.GE.TMIN.AND.DABS(ETK).LT.1.D-9) GO TO 200 IF(TSTEP.GE.TSTOP) GO TO 200 GO TO 100 200 CONTINUE CLOSE(20) CALL SBPRNT STOP END C ---------------------------------------------------------------------- C Mesh & FDM Parameters C ---------------------------------------------------------------------- SUBROUTINE SBMESH PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION BXP(NX),BYP(NY) DIMENSION BXUNEU(0:NX),BYVNEU(0:NY) DIMENSION BXUDIR(0:NX),BYVDIR(0:NY) COMMON /XYU/XU(0:NX),YV(0:NY) /XYP/XP(NX),YP(NY) COMMON /DXU/DXU(0:NX),DYV(0:NY) /DXP/DXP(NX),DYP(NY) COMMON /BXFP/BMXFP(NX),BPXFP(NX) /BYFP/BMYFP(NY),BPYFP(NY) COMMON /B0FP/B0FP(NX,NY) COMMON /CXFU/CMXFU(NX-1),CPXFU(NX-1) /CYFU/CMYFU(NY),CPYFU(NY) COMMON /CXFV/CMXFV(NX),CPXFV(NX) /CYFV/CMYFV(NY-1),CPYFV(NY-1) COMMON /C0/C0FU(NX-1,NY),C0FV(NX,NY-1) COMMON /DT/DT /RE/RE /DXMIN/DXMIN /JALG/JALG REB2=1.D0/RE/2.D0 XU(0)=0. XU(NX)=1.D0 YV(0)=0. YV(NY)=1.D0 C +++ UNIFORM +++ IF(JALG.EQ.0) THEN DO 111 I=1,NX-1 XU(I)=1.D0/DBLE(NX)*DBLE(I) 111 CONTINUE DO 211 J=1,NY-1 YV(J)=1.D0/DBLE(NY)*DBLE(J) 211 CONTINUE C +++ COS +++ ELSE PAI=DACOS(-1.D0) IF(JALG.EQ.1000) THEN DO 112 I=1,NX-1 XU(I)=(1.D0-DCOS(PAI*DBLE(I)/DBLE(NX)))/2.D0 112 CONTINUE DO 212 J=1,NY-1 YV(J)=(1.D0-DCOS(PAI*DBLE(J)/DBLE(NY)))/2.D0 212 CONTINUE C +++ TANH +++ ELSE ALG=DBLE(JALG)/1.D3 AT=DLOG((1.D0+ALG)/(1.D0-ALG))/2.D0 DO 110 I=1,NX-1 ETA=AT*(-1.D0+2.D0*DBLE(I)/DBLE(NX)) XU(I)=(DTANH(ETA)/ALG+1.D0)/2.D0 110 CONTINUE DO 210 J=1,NY-1 ETA=AT*(-1.D0+2.D0*DBLE(J)/DBLE(NY)) YV(J)=(DTANH(ETA)/ALG+1.D0)/2.D0 210 CONTINUE ENDIF ENDIF C WRITE(6,1110) XU C 1110 FORMAT(8F9.5) DO 120 I=1,NX XP(I)=(XU(I-1)+XU(I))/2.D0 DXP(I)=-XU(I-1)+XU(I) BXP(I)=1.D0/DXP(I) 120 CONTINUE DXU(0)=DXP(1) BXUNEU(0)=0.D0 BXUDIR(0)=1.D0/DXU(0) DO 130 I=1,NX-1 DXU(I)=-XP(I)+XP(I+1) BXUNEU(I)=1.D0/DXU(I) BXUDIR(I)=1.D0/DXU(I) 130 CONTINUE DXU(NX)=DXP(NX) BXUNEU(NX)=0.D0 BXUDIR(NX)=1.D0/DXU(NX) DO 140 I=1,NX BMXFP(I)=BXP(I)*BXUNEU(I-1) BPXFP(I)=BXP(I)*BXUNEU(I) CMXFV(I)=REB2*BXP(I)*BXUDIR(I-1) CPXFV(I)=REB2*BXP(I)*BXUDIR(I) 140 CONTINUE DO 150 I=1,NX-1 CMXFU(I)=REB2*BXUDIR(I)*BXP(I) CPXFU(I)=REB2*BXUDIR(I)*BXP(I+1) 150 CONTINUE DO 220 J=1,NY YP(J)=(YV(J-1)+YV(J))/2.D0 DYP(J)=-YV(J-1)+YV(J) BYP(J)=1.D0/DYP(J) 220 CONTINUE DYV(0)=DYP(1) BYVNEU(0)=0.D0 BYVDIR(0)=1.D0/DYV(0) DO 230 J=1,NY-1 DYV(J)=-YP(J)+YP(J+1) BYVNEU(J)=1.D0/DYV(J) BYVDIR(J)=1.D0/DYV(J) 230 CONTINUE DYV(NY)=DYP(NY) BYVNEU(NY)=0.D0 BYVDIR(NY)=1.D0/DYV(NX) DO 240 J=1,NY BMYFP(J)=BYP(J)*BYVNEU(J-1) BPYFP(J)=BYP(J)*BYVNEU(J) CMYFU(J)=REB2*BYP(J)*BYVDIR(J-1) CPYFU(J)=REB2*BYP(J)*BYVDIR(J) 240 CONTINUE DO 250 J=1,NY-1 CMYFV(J)=REB2*BYVDIR(J)*BYP(J) CPYFV(J)=REB2*BYVDIR(J)*BYP(J+1) 250 CONTINUE DO 300 J=1,NY DO 300 I=1,NX B0FP(I,J)=-(BMXFP(I)+BPXFP(I)+BMYFP(J)+BPYFP(J)) 300 CONTINUE DO 310 J=1,NY DO 310 I=1,NX-1 C0FU(I,J)=-(CMXFU(I)+CPXFU(I)+CMYFU(J)+CPYFU(J)) 310 CONTINUE DO 320 J=1,NY-1 DO 320 I=1,NX C0FV(I,J)=-(CMXFV(I)+CPXFV(I)+CMYFV(J)+CPYFV(J)) 320 CONTINUE DXMIN=DMIN1(DXP(1),DYP(1)) DXMAX=DMAX1(DXP(NX/2),DYP(NY/2)) WRITE(6,400) DXMIN,DXMAX 400 FORMAT('DXMIN=',F10.6,5X,'DXMAX=',F10.6) RETURN END C ---------------------------------------------------------------------- C Data Initializing C ---------------------------------------------------------------------- SUBROUTINE SBINIT PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /STEP/ISTEP,TSTEP /E/ENERGY COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /P/P(NX,NY) ISTEP=0 TSTEP=0.D0 DO 10 J=1,NY DO 10 I=0,NX 10 U(I,J)=0.D0 DO 20 J=0,NY DO 20 I=1,NX 20 V(I,J)=0.D0 DO 30 J=1,NY DO 30 I=1,NX 30 P(I,J)=0.D0 ENERGY=0.D0 RETURN END C ---------------------------------------------------------------------- C Boundary Condition C ---------------------------------------------------------------------- SUBROUTINE SBWALL PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) DO 10 I=0,NX U(I,0)=-U(I,1) U(I,NY+1)=2.D0-U(I,NY) 10 CONTINUE DO 20 J=0,NY V(0,J)=-V(1,J) V(NX+1,J)=-V(NX,J) 20 CONTINUE RETURN END C ---------------------------------------------------------------------- C Momentum Flux (DIVERGENCE FORM) C ---------------------------------------------------------------------- SUBROUTINE SBFLUX PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) COMMON /UP/UP(NX,NY) /VP/VP(NX,NY) COMMON /UR/UR(NX-1,0:NY) /VR/VR(0:NX,NY-1) /IFLX/IFLX COMMON /UVF/UF(NX-1,NY),VF(NX,NY-1) COMMON /DXU/DXU(0:NX),DYV(0:NY) /DXP/DXP(NX),DYP(NY) DO 10 J=1,NY DO 10 I=1,NX UC=(U(I-1,J)+U(I,J))/2.D0 VC=(V(I,J-1)+V(I,J))/2.D0 UP(I,J)=-UC**2 VP(I,J)=-VC**2 10 CONTINUE C SHEAR FLUX AT THE WALL DO 20 J=1,NY-1 VR( 0,J)=0.D0 VR(NX,J)=0.D0 20 CONTINUE DO 22 I=1,NX-1 UR(I, 0)=0. UR(I,NY)=0. 22 CONTINUE IF(IFLX.EQ.3) GO TO 300 IF(IFLX.EQ.1) THEN DO 30 J=1,NY-1 DO 30 I=1,NX-1 UC=(U(I,J)+U(I,J+1))/2.D0 VC=(V(I,J)+V(I+1,J))/2.D0 UR(I,J)=-UC*VC 30 CONTINUE ELSE DO 40 J=1,NY-1 DO 40 I=1,NX-1 UC=(DYP(J+1)*U(I,J)+DYP(J)*U(I,J+1))/(2.D0*DYV(J)) VC=(DXP(I+1)*V(I,J)+DXP(I)*V(I+1,J))/(2.D0*DXU(I)) UR(I,J)=-UC*VC 40 CONTINUE ENDIF DO 50 J=1,NY-1 DO 50 I=1,NX-1 VR(I,J)=UR(I,J) 50 CONTINUE GO TO 400 300 CONTINUE DO 60 J=1,NY-1 DO 60 I=1,NX-1 UC=(U(I,J)+U(I,J+1))/2.D0 VC=(V(I,J)+V(I+1,J))/2.D0 UJ=(DYP(J)*U(I,J)+DYP(J+1)*U(I,J+1))/(2.D0*DYV(J)) VJ=(DXP(I)*V(I,J)+DXP(I+1)*V(I+1,J))/(2.D0*DXU(I)) UR(I,J)=-UC*VJ VR(I,J)=-UJ*VC 60 CONTINUE 400 CONTINUE DO 100 J=1,NY DO 100 I=1,NX-1 UF(I,J)=(-UP(I,J)+UP(I+1,J))/DXU(I)+(-UR(I,J-1)+UR(I,J))/DYP(J) 100 CONTINUE DO 110 J=1,NY-1 DO 110 I=1,NX VF(I,J)=(-VR(I-1,J)+VR(I,J))/DXP(I)+(-VP(I,J)+VP(I,J+1))/DYV(J) 110 CONTINUE RETURN END C ---------------------------------------------------------------------- C Momentum Flux (GRADIENT FORM) C ---------------------------------------------------------------------- SUBROUTINE SBCONV PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) COMMON /UP/UP(NX,NY) /VP/VP(NX,NY) COMMON /UR/UR(NX-1,0:NY) /VR/VR(0:NX,NY-1) /IFLX/IFLX COMMON /UVF/UF(NX-1,NY),VF(NX,NY-1) COMMON /DXU/DXU(0:NX),DYV(0:NY) /DXP/DXP(NX),DYP(NY) C CONVECTION TERMS AT THE WALL DO 80 J=1,NY-1 VR( 0,J)=0.D0 VR(NX,J)=0.D0 80 CONTINUE DO 82 I=1,NX-1 UR(I, 0)=0. UR(I,NY)=0. 82 CONTINUE IF(IFLX.EQ.3) GO TO 300 IF(IFLX.EQ.1) THEN DO 10 J=1,NY DO 10 I=1,NX UP(I,J)=-(U(I-1,J)+U(I,J))/2.D0*(-U(I-1,J)+U(I,J))/DXP(I) VP(I,J)=-(V(I,J-1)+V(I,J))/2.D0*(-V(I,J-1)+V(I,J))/DYP(J) 10 CONTINUE DO 12 J=1,NY-1 DO 12 I=1,NX-1 UR(I,J)=-(V(I,J)+V(I+1,J))/2.D0*(-U(I,J)+U(I,J+1))/DYV(J) VR(I,J)=-(U(I,J)+U(I,J+1))/2.D0*(-V(I,J)+V(I+1,J))/DXU(I) 12 CONTINUE DO 14 J=1,NY DO 14 I=1,NX-1 UF(I,J)=(UP(I,J)+UP(I+1,J))/2.D0+(UR(I,J-1)+UR(I,J))/2.D0 14 CONTINUE DO 16 J=1,NY-1 DO 16 I=1,NX VF(I,J)=(VR(I-1,J)+VR(I,J))/2.D0+(VP(I,J)+VP(I,J+1))/2.D0 16 CONTINUE ELSE DO 20 J=1,NY DO 20 I=1,NX UP(I,J)=-(U(I-1,J)+U(I,J))/2.D0*(-U(I-1,J)+U(I,J))/DXP(I) VP(I,J)=-(V(I,J-1)+V(I,J))/2.D0*(-V(I,J-1)+V(I,J))/DYP(J) 20 CONTINUE DO 22 J=1,NY-1 DO 22 I=1,NX-1 UC=(DYP(J+1)*U(I,J)+DYP(J)*U(I,J+1))/(2.D0*DYV(J)) VC=(DXP(I+1)*V(I,J)+DXP(I)*V(I+1,J))/(2.D0*DXU(I)) UR(I,J)=-VC*(-U(I,J)+U(I,J+1))/DYV(J) VR(I,J)=-UC*(-V(I,J)+V(I+1,J))/DXU(I) 22 CONTINUE DO 24 J=1,NY DO 24 I=1,NX-1 UF(I,J)=(DXP(I+1)*UP(I,J)+DXP(I)*UP(I+1,J))/(2.D0*DXU(I)) &+(UR(I,J-1)+UR(I,J))/2.D0 24 CONTINUE DO 26 J=1,NY-1 DO 26 I=1,NX VF(I,J)=(VR(I-1,J)+VR(I,J))/2.D0 &+(DYP(J+1)*VP(I,J)+DYP(J)*VP(I,J+1))/(2.D0*DYV(J)) 26 CONTINUE ENDIF GO TO 400 300 CONTINUE DO 30 J=1,NY DO 30 I=1,NX UP(I,J)=-(U(I-1,J)+U(I,J))/2.D0*(-U(I-1,J)+U(I,J)) VP(I,J)=-(V(I,J-1)+V(I,J))/2.D0*(-V(I,J-1)+V(I,J)) 30 CONTINUE DO 32 J=1,NY-1 DO 32 I=1,NX-1 UC=(DYP(J)*U(I,J)+DYP(J+1)*U(I,J+1))/(2.D0*DYV(J)) VC=(DXP(I)*V(I,J)+DXP(I+1)*V(I+1,J))/(2.D0*DXU(I)) UR(I,J)=-VC*(-U(I,J)+U(I,J+1)) VR(I,J)=-UC*(-V(I,J)+V(I+1,J)) 32 CONTINUE DO 34 J=1,NY DO 34 I=1,NX-1 UF(I,J)=(UP(I,J)+UP(I+1,J))/(2.D0*DXU(I)) & +(UR(I,J-1)+UR(I,J))/(2.D0*DYP(J)) 34 CONTINUE DO 36 J=1,NY-1 DO 36 I=1,NX VF(I,J)=(VR(I-1,J)+VR(I,J))/(2.D0*DXP(I)) & +(VP(I,J)+VP(I,J+1))/(2.D0*DYV(J)) 36 CONTINUE 400 CONTINUE RETURN END C ---------------------------------------------------------------------- C R.H.S. of Elliptic Equation for Velocity C ---------------------------------------------------------------------- SUBROUTINE SBPRED PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /P/P(NX,NY) COMMON /UP/QU(NX,NY) /VP/QV(NX,NY) COMMON /UVF/UF(NX-1,NY),VF(NX,NY-1) /UVB/UB(NX-1,NY),VB(NX,NY-1) COMMON /DXU/DXU(0:NX),DYV(0:NY) COMMON /CXFU/CMXFU(NX-1),CPXFU(NX-1) /CYFU/CMYFU(NY),CPYFU(NY) COMMON /CXFV/CMXFV(NX),CPXFV(NX) /CYFV/CMYFV(NY-1),CPYFV(NY-1) COMMON /C0/C0FU(NX-1,NY),C0FV(NX,NY-1) COMMON /STEP/ISTEP,TSTEP /BT/BT IF(ISTEP.EQ.0) THEN DO 10 J=1,NY DO 10 I=1,NX-1 UB(I,J)=UF(I,J) 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX VB(I,J)=VF(I,J) 20 CONTINUE ENDIF DO 50 J=1,NY DO 50 I=1,NX-1 QU(I,J)=-(-P(I,J)+P(I+1,J))/DXU(I)+(3.D0*UF(I,J)-UB(I,J))/2.D0 & +CMYFU(J)*U(I,J-1)+CMXFU(I)*U(I-1,J) & +(BT+C0FU(I,J))*U(I,J) & +CPXFU(I)*U(I+1,J)+CPYFU(J)*U(I,J+1) 50 CONTINUE DO 60 J=1,NY-1 DO 60 I=1,NX QV(I,J)=-(-P(I,J)+P(I,J+1))/DYV(J)+(3.D0*VF(I,J)-VB(I,J))/2.D0 & +CMYFV(J)*V(I,J-1)+CMXFV(I)*V(I-1,J) & +(BT+C0FV(I,J))*V(I,J) & +CPXFV(I)*V(I+1,J)+CPYFV(J)*V(I,J+1) 60 CONTINUE IF(ISTEP.GT.0) THEN DO 30 J=1,NY DO 30 I=1,NX-1 UB(I,J)=UF(I,J) 30 CONTINUE DO 40 J=1,NY-1 DO 40 I=1,NX VB(I,J)=VF(I,J) 40 CONTINUE ENDIF RETURN END C ---------------------------------------------------------------------- C Solving Elliptic Equation for Velocity C ---------------------------------------------------------------------- SUBROUTINE SBSORU PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /P/P(NX,NY) COMMON /UP/QU(NX,NY) /VP/QV(NX,NY) COMMON /CXFU/CMXFU(NX-1),CPXFU(NX-1) /CYFU/CMYFU(NY),CPYFU(NY) COMMON /CXFV/CMXFV(NX),CPXFV(NX) /CYFV/CMYFV(NY-1),CPYFV(NY-1) COMMON /C0/C0FU(NX-1,NY),C0FV(NX,NY-1) COMMON /BT/BT /ITRU/ITRU /ERRU/ERRU ITRUE=NX SUMQU=0.D0 SUMQV=0.D0 DO 10 J=1,NY DO 10 I=1,NX-1 SUMQU=SUMQU+QU(I,J)**2 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX SUMQV=SUMQV+QV(I,J)**2 20 CONTINUE ITRU=0 100 ITRU=ITRU+1 SUMRU=0.D0 SUMRV=0.D0 DO 50 J=1,NY DO 50 I=1,NX-1 RESIU=-CMYFU(J)*U(I,J-1)-CMXFU(I)*U(I-1,J) & +(BT-C0FU(I,J))*U(I,J) & -CPXFU(I)*U(I+1,J)-CPYFU(J)*U(I,J+1)-QU(I,J) U(I,J)=U(I,J)-1.5*RESIU/(BT-C0FU(I,J)) SUMRU=SUMRU+RESIU**2 50 CONTINUE DO 60 J=1,NY-1 DO 60 I=1,NX RESIV=-CMYFV(J)*V(I,J-1)-CMXFV(I)*V(I-1,J) & +(BT-C0FV(I,J))*V(I,J) & -CPXFV(I)*V(I+1,J)-CPYFV(J)*V(I,J+1)-QV(I,J) V(I,J)=V(I,J)-1.5*RESIV/(BT-C0FV(I,J)) SUMRV=SUMRV+RESIV**2 60 CONTINUE ERRU=DSQRT((SUMRU+SUMRV)/(SUMQU+SUMQV)) IF(ITRU.GE.ITRUE) GO TO 200 IF(ERRU.LT.1.D-8) GO TO 200 GO TO 100 200 CONTINUE RETURN END C ---------------------------------------------------------------------- C R.H.S. of Poisson Equation for Pressure C ---------------------------------------------------------------------- SUBROUTINE SBRHSP PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /DIV/DIV(NX,NY) COMMON /UP/Q(NX,NY) /SUMQ/SUMQ /DXP/DXP(NX),DYP(NY) /DT/DT DO 10 J=1,NY DO 10 I=1,NX DIV(I,J)=(-U(I-1,J)+U(I,J))/DXP(I)+(-V(I,J-1)+V(I,J))/DYP(J) 10 CONTINUE SUMQ=0.D0 DO 30 J=1,NY DO 30 I=1,NX Q(I,J)=DIV(I,J)/DT SUMQ=SUMQ+Q(I,J)**2 30 CONTINUE RETURN END C ---------------------------------------------------------------------- C Solving Poisson Equation for Pressure C ---------------------------------------------------------------------- SUBROUTINE SBSORP PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /UP/Q(NX,NY) /S/S(0:NX+1,0:NY+1) COMMON /SUMQ/SUMQ /SUMR/SUMR /ERRP/ERRP /ITRP/ITRP COMMON /BXFP/BMXFP(NX),BPXFP(NX) /BYFP/BMYFP(NY),BPYFP(NY) COMMON /B0FP/B0FP(NX,NY) ITRPE=2*NX ITRP=0 DO 10 J=0,NY+1 DO 10 I=0,NX+1 S(I,J)=0.D0 10 CONTINUE 100 ITRP=ITRP+1 SUMR=0.D0 DO 30 J=1,NY DO 30 I=1,NX RESI=BMYFP(J)*S(I,J-1)+BMXFP(I)*S(I-1,J)+B0FP(I,J)*S(I,J) & +BPXFP(I)*S(I+1,J)+BPYFP(J)*S(I,J+1)-Q(I,J) S(I,J)=S(I,J)-1.5*RESI/B0FP(I,J) SUMR=SUMR+RESI**2 30 CONTINUE ERRP=DSQRT(SUMR/SUMQ) IF(ITRP.GE.ITRPE) GO TO 200 IF(ERRP.LT.1.D-8) GO TO 200 GO TO 100 200 CONTINUE RETURN END C ---------------------------------------------------------------------- C Correction Step C ---------------------------------------------------------------------- SUBROUTINE SBCORR PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /P/P(NX,NY) COMMON /S/S(0:NX+1,0:NY+1) /DT/DT /DXU/DXU(0:NX),DYV(0:NY) COMMON /GP/GPX(0:NX),GPY(NX,0:NY) DO 10 J=1,NY DO 10 I=1,NX-1 U(I,J)=U(I,J)-DT*(-S(I,J)+S(I+1,J))/DXU(I) 10 CONTINUE DO 20 J=1,NY-1 DO 20 I=1,NX V(I,J)=V(I,J)-DT*(-S(I,J)+S(I,J+1))/DYV(J) 20 CONTINUE DO 30 J=1,NY DO 30 I=1,NX P(I,J)=P(I,J)+S(I,J) 30 CONTINUE RETURN END C ---------------------------------------------------------------------- C Kinetic Energy C ---------------------------------------------------------------------- SUBROUTINE SBENGY PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) COMMON /E/ENERGY /DXP/DXP(NX),DYP(NY) SUME=0.D0 DO 10 J=1,NY DO 10 I=1,NX US=(U(I-1,J)+U(I,J))/2.D0 VS=(V(I,J-1)+V(I,J))/2.D0 SUME=SUME+DXP(I)*DYP(J)*(US**2+VS**2) 10 CONTINUE ENERGY=SUME/2.D0 RETURN END C ---------------------------------------------------------------------- C Divergence check C ---------------------------------------------------------------------- SUBROUTINE SBDIVU PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) COMMON /DIVMAX/DIVMAX /DXP/DXP(NX),DYP(NY) /DIV/DIV(NX,NY) DIVMAX=0.D0 DO 90 J=1,NY DO 90 I=1,NX DIV(I,J)=(-U(I-1,J)+U(I,J))/DXP(I)+(-V(I,J-1)+V(I,J))/DYP(J) DIVMAX=DMAX1(DABS(DIV(I,J)),DIVMAX) 90 CONTINUE RETURN END C ---------------------------------------------------------------------- C Output C ---------------------------------------------------------------------- SUBROUTINE SBPRNT PARAMETER (NX=32,NY=32) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*10 FILENAME CHARACTER*28 FLUX CHARACTER*20 FRMK COMMON /U/U(0:NX,0:NY+1) /V/V(0:NX+1,0:NY) /P/P(NX,NY) COMMON /XYP/XP(NX),YP(NY) /DIV/DIV(NX,NY) /DXP/DXP(NX),DYP(NY) COMMON /STEP/ISTEP,TSTEP /RE/RE /DT/DT /FILE/FILENAME COMMON /DIVMAX/DIVMAX /E/ENERGY /ETK/ETK /FLUX/FLUX COMMON /FRMK/FRMK OPEN(10,file=FILENAME//'.dat',status='new',form='formatted') WRITE(10,900) FRMK,FLUX,INT(RE),NX,NY &,DXP(1),DXP(NX/2),DXP(NX/2)/DXP(1) &,DT,TSTEP,DIVMAX,ENERGY,ETK 900 FORMAT('PROGRAM:CAVITY4 (2ND-AB/CN IMPLICIT METHOD)' &/'FORM FOR CONVECTION = ',A20 &/'INTERPOLATION METHOD = ',A28 &/'REYNOLDS NUMBER: RE =',I8 &/'NUMBER OF GRID POINTS: NX =',I4,', NY =',I4 &/'GRID SPACING: DXMIN=',F9.6,' DXMAX=',F9.6,' RATIO=',F7.3 &/'TIME-INCREMENT: DT =',E12.3,', TIME: TSTEP =',F10.5 &/'DIVERGENCE ERROR (MAX): DIVMAX =',E12.3 &/'KINETIC ENERGY: ENERGY =',F15.10,' (DE/DT)/E=',E12.3) WRITE(10,1100) DO 10 L=1,NX I=L J=NY+1-L WRITE(10,1000) & J,YP(J),U(NX/2,J),(P(NX/2,J)+P(NX/2+1,J))/2.D0, & I,XP(I),V(I,NY/2),(P(I,NY/2)+P(I,NY/2+1))/2.D0 10 CONTINUE CLOSE(10) 1000 FORMAT(2(I5,3F10.6)) 1100 FORMAT(4X,1HJ,9X,1HY,9X,1HU,9X,1HP, & 4X,1HI,9X,1HX,9X,1HV,9X,1HP) OPEN(30,file=FILENAME//'.d',status='new',form='unformatted') WRITE(30) RE,TSTEP WRITE(30) U WRITE(30) V WRITE(30) P CLOSE(30) RETURN END