2次元キャビティ流れNO.1(FORTRAN)プログラム
2003/09/02 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元キャビティ流れを流れ関数-渦度法により差分法で計算する。
2. プログラム
C=================================================================C
C 2-Dimensional Cavity Flow by Stream function and Vorticity C
C with FDM and SOR C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/02 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
PARAMETER(NDX=300, NDY=300)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PHI(0:NDX,0:NDY), OMG(0:NDX,0:NDY)
DIMENSION U(0:NDX,0:NDY), V(0:NDX,0:NDY)
CHARACTER*11 FLOW1
DATA FLOW1/'FLOW1-?.txt'/
C Initial Data
ID = 0
CALL XCLOCK(T1,3)
WRITE(6,*) 'Type In NX,Re'
READ(5,*) NX,Re
WRITE(6,*) 'DT,NT(Total NO.),MT(Output Interval)'
READ(5,*) DT,NT,MT
if(NX.gt.NDX) NX = NDX
NY = NX
WRITE(6,100) NX,NY,Re,DT,NT,MT
C Initial Constant
T = 0.0
DH = NX
H = 1.0/DH
H2 = H**2
D2 = DH**2
C Set Initial Condition
do j=0,NY
do i=0,NX
OMG(i,j) = 0.0
PHI(i,j) = 0.0
end do
end do
C SOR Parameter(ALP)
ALP = 1.0 + DLOG(NX*1.0D0)/DLOG(NDX*1.2D0)
EPS = 1.0e-4
C Main Loop
CALL XCLOCK(T1,5)
do k=1,NT
T = T + DT
C Set Boundary Condition
do i=0,NX
OMG(i,0) = -2.0*PHI(i,1)*D2
OMG(i,NY) = -2.0*(PHI(i,NY-1) + H)*D2
end do
do j=1,NY-1
OMG(0,j) = -2.0*PHI(1,j)*D2
OMG(NX,j) = -2.0*PHI(NX-1,j)*D2
end do
C Compute Omega
do j=1,NY-1
do i=1,NX-1
OMG(i,j) = OMG(i,j) + DT*D2*( (
1 - (PHI(i,j+1)-PHI(i,j-1))*(OMG(i+1,j)-OMG(i-1,j))
2 + (PHI(i+1,j)-PHI(i-1,j))*(OMG(i,j+1)-OMG(i,j-1))
3 ) / 4.0 + (OMG(i,j-1)+OMG(i-1,j)-4.0*OMG(i,j)
4 + OMG(i+1,j)+OMG(i,j+1)) / Re )
end do
end do
C Compute Phi
CALL PHISOR(PHI,OMG,NX,NY,NDX,H2,ALP,EPS)
C Output U,V,Phi
if(MOD(k,MT).eq.0) then
FLOW1(7:7) = CHAR(ID+ICHAR('0'))
OPEN(UNIT=1,FILE=FLOW1)
WRITE(1,100) NX,NY,Re,DT,NT,MT
WRITE(1,150) T,k
WRITE(6,150) T,k
CALL OUTUVP(PHI,U,V,NX,NY,NDX)
CLOSE(UNIT=1)
if(ID.lt.9) ID = ID + 1
end if
end do
CALL XCLOCK(T2,5)
WRITE(6,200) NX,NT,T2-T1
C
STOP
100 FORMAT('# Cavity Flow Numerical Analysis',/'# NX,NY=',2I5,
1' Re=',F6.0,' DT=',F6.2,' NT,MT=',2I5)
150 FORMAT('#Time=',F7.3,' Step=',I5)
200 FORMAT(1H ,' NX,NT,Time(s)=',2I5,F9.2)
END
C=================================================================C
SUBROUTINE PHISOR(PHI,OMG,NX,NY,NDX,H2,ALP,EPS)
C=================================================================C
C Solve Ax=b by SOR with 2 dimensional FDM C
C Given Omega ( Acceleration factor ) C
C-----------------------------------------------------------------C
C PHI(0:NDX,0:NY) R*8, I/O, Phi Value C
C OMG(0:NDX,0:NY) R*8, In, Omega Value C
C NX I*4, In, Grid Numbers on X-axis C
C NY I*4, In, Grid Numbers on Y-axis C
C NDX I*4, In, First Array Size of PHI,OMG C
C H2 R*8, In, H2=H**2 C
C ALP R*8, In, SOR Acceleration factor C
C EPS R*8, In, if ||r||/||b|| <= EPS --> return C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/02 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PHI(0:NDX,0:NY), OMG(0:NDX,0:NY)
C Get 2-Norm B=D2*OMG
BN = 0.0
do j=1,NY-1
do i=1,NX-1
BN = BN + (H2*OMG(i,j))**2
end do
end do
C Main Loop
do k=1,NX*NY
RN = 0.0
do j=1,NY-1
do i=1,NX-1
R = (H2*OMG(i,j) + PHI(i,j-1) + PHI(i-1,j) + PHI(i+1,j)
1 + PHI(i,j+1) )/4.0 - PHI(i,j)
PHI(i,j) = PHI(i,j) + ALP*R
RN = RN + (R*4.0)**2
end do
end do
C if(ERR <= EPS) return
ERR = DSQRT(RN/BN)
if(ERR.le.EPS) goto 100
end do
100 CONTINUE
C
RETURN
END
C=================================================================C
SUBROUTINE OUTUVP(PHI,U,V,NX,NY,NDX)
C=================================================================C
C Compute U,V and Output U,V,Phi C
C U=d(PHI)/dy, V=-d(PHI)/dx C
C-----------------------------------------------------------------C
C PHI(0:NDX,0:NY) R*8, In, Phi Value C
C U(0:NDX,0:NY) R*8, Out, U Value C
C V(0:NDX,0:NY) R*8, Out, V Value C
C NX I*4, In, Grid Numbers on X-axis C
C NY I*4, In, Grid Numbers on Y-axis C
C NDX I*4, In, First Array Size of PHI,OMG C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/02 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION PHI(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
C Initial
XV = NX/2.0D0
YV = NY/2.0D0
C Compute U,V
do j=1,NY-1
do i=1,NX-1
U(i,j) = (PHI(i,j+1) - PHI(i,j-1))*YV
V(i,j) = (PHI(i-1,j) - PHI(i+1,j))*XV
end do
end do
C Boundary on Y=0,1
do i=0,NX
U(i,NY) = 1.0
V(i,NY) = 0.0
U(i,0) = 0.0
V(i,0) = 0.0
end do
C Boundary on X=0,1
do j=1,NY-1
U(0,j) = 0.0
V(0,j) = 0.0
U(NX,j) = 0.0
V(NX,j) = 0.0
end do
C Output
DX = 1.0D0/NX
DY = 1.0D0/NY
do j=0,NY
Y = j*DY
do i=0,NX
X = i*DX
WRITE(1,100) X,Y,PHI(i,j),U(i,j),V(i,j)
end do
end do
C
RETURN
100 FORMAT(1H ,5F9.5)
END
C=================================================================C
SUBROUTINE XCLOCK(CPU,ID)
C=================================================================C
C CPU time Subroutine C
C CPU R*8 Out, CPU Time C
C ID I*4 In, Dummy ( Same Hitachi FORTRAN XCLOCK ) C
C-----------------------------------------------------------------C
C Written by Y. Ushiro ( Waseda and Hitachi ) 2002.10.29 C
C=================================================================C
REAL*8 CPU
INTEGER*2 I1, I2, I3, I4
C
IF(ID.GE.1) THEN
CALL GETTIM(I1,I2,I3,I4)
CPU = ( I1*60.0 + I2 )*60.0 + I3 + I4*0.01
END IF
C
RETURN
END