2次元キャビティ流れNO.2(FORTRAN)プログラム
2003/09/10 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元キャビティ流れを速度-圧力法により差分法で計算する。
2. プログラム
C=================================================================C
C 2-Dimensional Cavity Flow by velocity and pressure C
C with FDM and SOR C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/10 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
PARAMETER(NDX=300, NDY=300)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(0:NDX,0:NDY), V(0:NDX,0:NDY), P(0:NDX,0:NDY)
DIMENSION B(0:NDX,0:NDY)
CHARACTER*12 FLOW2
DATA FLOW2/'FLOW2-??.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
D5 = DH/2.0
H = 1.0/DH
H2 = H**2
D2 = DH**2
C Set Initial Condition
do j=0,NY
do i=0,NX
U(i,j) = 0.0
V(i,j) = 0.0
P(i,j) = 0.0
end do
end do
C U=1.0 on L1
do i=0,NX
U(i,NY) = 1.0
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)
NITER = 0
do k=1,NT
T = T + DT
C Compute U,V
OVER = 0.0
do j=1,NY-1
do i=1,NX-1
U(i,j) = U(i,j) + DT*( ( U(i-1,j)**2 - U(i+1,j)**2
1 + U(i,j-1)*V(i,j-1) - U(i,j+1)*V(i,j+1)
2 + P(i-1,j) - P(i+1,j) )*D5 - ( 4.0*U(i,j) - U(i,j-1)
3 - U(i-1,j) - U(i+1,j) - U(i,j+1) )*D2/Re )
V(i,j) = V(i,j) + DT*( ( V(i,j-1)**2 - V(i,j+1)**2
1 + U(i-1,j)*V(i-1,j) - U(i+1,j)*V(i+1,j)
2 + P(i,j-1) - P(i,j+1) )*D5 - ( 4.0*V(i,j) - V(i,j-1)
3 - V(i-1,j) - V(i+1,j) - V(i,j+1) )*D2/Re )
OVER = OVER + abs(U(i,j)) + abs(V(i,j))
end do
end do
C Check Divergent
if(OVER.gt.1.0D10) then
WRITE(6,*) '** Stop for over flow computation **'
WRITE(6,*) ' You have to give smaller DT'
STOP
end if
C Compute P
CALL PSOR(P,U,V,B,NX,NY,NDX,DT,H,ALP,EPS,ITER)
NITER = NITER + ITER
C Output U,V,Phi
if(MOD(k,MT).eq.0) then
NI1 = ID/10
NI2 = ID - NI1*10
FLOW2(7:7) = CHAR(NI1+ICHAR('0'))
FLOW2(8:8) = CHAR(NI2+ICHAR('0'))
OPEN(UNIT=1,FILE=FLOW2)
WRITE(1,100) NX,NY,Re,DT,NT,MT
WRITE(1,150) T,k,NITER
WRITE(6,150) T,k,NITER
NITER = 0
CALL OUTUVP(P,U,V,NX,NY,NDX)
CLOSE(UNIT=1)
if(ID.lt.99) 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,' SOR Loops=',I6)
200 FORMAT(1H ,' NX,NT,Time(s)=',2I5,F9.2)
END
C=================================================================C
SUBROUTINE PSOR(P,U,V,B,NX,NY,NDX,DT,H,ALP,EPS,ITER)
C=================================================================C
C Solve Ax=b by SOR with 2 dimensional FDM C
C Get P with U,V ( Acceleration factor ) C
C-----------------------------------------------------------------C
C P(0:NDX,0:NY) R*8, I/O, P vector C
C U(0:NDX,0:NY) R*8, In, U vector C
C V(0:NDX,0:NY) R*8, In, V vector C
C B(0:NDX,0:NDY) R*8, Wk, Work 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 DT R*8, In, Delta T C
C H R*8, In, H=1.0/NX C
C ALP R*8, In, SOR Acceleration factor C
C EPS R*8, In, if ||r||/||b|| <= EPS --> return C
C ITER I*4, Out, Number of Iteration C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/10 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
DIMENSION B(0:NDX,0:NY)
C Set B and BN (B norm)
BN = 0.0
do j=0,NY
do i=0,NX
C Boundary on L1
if(j.ne.NY) then
U0P = U(i,j+1)
V0P = V(i,j+1)
else
U0P = 1.0
V0P = 0.0
end if
C Boundary on L2
if(i.ne.0) then
UM0 = U(i-1,j)
VM0 = V(i-1,j)
else
UM0 = U(1,j)
VM0 = -V(1,j)
end if
C Boundary on L3
if(j.ne.0) then
U0M = U(i,j-1)
V0M = V(i,j-1)
else
U0M = -U(i,1)
V0M = V(i,1)
end if
C Boundary on L4
if(i.ne.NX) then
UP0 = U(i+1,j)
VP0 = V(i+1,j)
else
UP0 = U(NX-1,j)
VP0 = -V(NX-1,j)
end if
C Computation
B(i,j) = ( (UP0 - UM0)**2 + (V0P - V0M)**2 )/4.0
1 + (VP0 - VM0)*(U0P - U0M)/2.0
2 - (UP0 - UM0 + V0P - V0M)*H/(2.0*DT)
BN = BN + B(i,j)**2
end do
end do
C Main Loop
do k=1,NX*NY
RN = 0.0
do j=0,NY
JP = j + 1
JM = j - 1
C Boundary on L1,L3
if(j.eq.NY) JP = NY - 1
if(j.eq.0) JM = 1
do i=0,NX
IP = i + 1
IM = i - 1
C Boundary on L2,L4
if(i.eq.NX) IP = NX - 1
if(i.eq.0) IM = 1
C Computation
R = (B(i,j)+P(i,JM)+P(IM,j)+P(IP,j)+P(i,JP))/4.0 - P(i,j)
P(i,j) = P(i,j) + ALP*R
C Reset on P1
if(i.eq.0 .and. j.eq.0) then
P(i,j) = 0.0
R = 0.0
end if
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
ITER = k
C
RETURN
END
C=================================================================C
SUBROUTINE OUTUVP(P,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 P(0:NDX,0:NY) R*8, In, P vector C
C U(0:NDX,0:NY) R*8, Out, U vector C
C V(0:NDX,0:NY) R*8, Out, V vector 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/10 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(0:NDX,0:NY), U(0:NDX,0:NY), V(0:NDX,0:NY)
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,P(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