2次元角柱周りの流れ解析(FORTRAN)プログラム
2003/09/17 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------
1. 概要
2次元角柱周りの流れを速度-圧力法により差分法で計算する。
2. プログラム
C=================================================================C
C 2-Dimensional Flow around Square a Pillar by velocity C
C and pressure with FDM and SOR C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/14 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
PARAMETER(NDX=301, NDY=301)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(-1:NDX,-1:NDY), V(-1:NDX,-1:NDY)
DIMENSION A(-1:NDX,-1:NDY,5), B(-1:NDX,-1:NDY)
DIMENSION P(-1:NDX,-1:NDY), PHI(-1:NDX,-1:NDY)
CHARACTER*12 FLOW3
DATA FLOW3/'FLOW3-??.txt'/
COMMON /MESH/NX,NX1,NX2,NY,NY1,NY2
C Initial Data
ID = 0
CALL XCLOCK(T1,3)
XM = 2.0
YM = 1.0
WRITE(6,*) 'Type In NX,NY'
READ(5,*) NX,NY
WRITE(6,*) 'Type in Re,DT,NT(Total NO.),MT(Output Interval)'
READ(5,*) Re,dt,NT,MT
if(NX.ge.NDX) NX = NDX - 1
if(NY.ge.NDY) NY = NDY - 1
NX1 = NX*2/10
NX2 = NX*3/10
NY1 = NY*2/5
NY2 = NY*3/5
WRITE(6,100) XM,NX1,NX2,NX
WRITE(6,110) YM,NY1,NY2,NY
WRITE(6,120) Re,dt,NT,MT
C Initial Constant
T = 0.0
hx = XM/NX
hy = YM/NY
dx = 1.0/(2.0*hx)
dy = 1.0/(2.0*hy)
ddx = 1.0/hx**2
ddy = 1.0/hy**2
EPS = 1.0e-4
OMG = 1.0 + DLOG((NX+NY)*1.0D0)/DLOG((NDX+NDY)*1.2D0)
C Set Initial U,V
ITER = NX*NY
CALL PHISOL(PHI,A,B,U,V,hx,hy,OMG,NDX,NDY,EPS,ITER)
C Set Initial P
do j=-1,NY+1
do i=-1,NX+1
P(i,j) = 0.0
end do
end do
C Main Loop
CALL XCLOCK(T1,5)
NITER = ITER
do k=1,NT
T = T + DT
C Compute P
ITER = NX*NY
CALL PSOL(P,A,B,U,V,hx,hy,dt,Re,OMG,NDX,NDY,EPS,ITER)
NITER = NITER + ITER
C Compute U,V
OVER = 0.0
C Lower-Y Area
do j=1,NY1-1
do i=1,NX
include 'UV.h'
end do
end do
C Middle-Y Area
do j=NY1,NY2
do i=1,NX1-1
include 'UV.h'
end do
do i=NX2+1,NX
include 'UV.h'
end do
end do
C Upper-Y Area
do j=NY2+1,NY-1
do i=1,NX
include 'UV.h'
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 Output U,V,Phi
if(MOD(k,MT).eq.0) then
NI1 = ID/10
NI2 = ID - NI1*10
FLOW3(7:7) = CHAR(NI1+ICHAR('0'))
FLOW3(8:8) = CHAR(NI2+ICHAR('0'))
OPEN(UNIT=1,FILE=FLOW3)
WRITE(1,100) XM,NX1,NX2,NX
WRITE(1,110) YM,NY1,NY2,NY
WRITE(1,120) Re,dt,NT,MT
WRITE(1,150) T,k,NITER
WRITE(6,150) T,k,NITER
NITER = 0
CALL OUTUVP(P,U,V,PHI,hx,hy,NDX,NDY)
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('# Flow around Square a Pillar Numerical Analysis',/,
1 '# XM=',F8.1,' NX1,NX2,NX=',3I4)
110 FORMAT('# YM=',F8.1,' NY1,NY2,NY=',3I4)
120 FORMAT('# Re=',F8.1,' DT=',F7.3,' 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 PSOL(P,A,B,U,V,hx,hy,dt,Re,OMG,NDX,NDY,EPS,ITER)
C=================================================================C
C Solve Ax=b by SOR with 2 dimensional FDM C
C Solve P by U,V with -div(P)=f(U,V) C
C-----------------------------------------------------------------C
C P(-1:NDX,-1:NDY) R*8, I/O, P vector C
C A(-1:NDX,-1:NDY,5) R*8, Wk, Matrix A for -div(PHI)=0 C
C B(-1:NDX,-1:NDY) R*8, Wk, Right Vector for -div(PHI)=0 C
C U(-1:NDX,-1:NDY) R*8, In, U vector C
C V(-1:NDX,-1:NDY) R*8, In, V vector C
C hx R*8, In. Delta x C
C hy R*8, In, Delta y C
C dt R*8, In, Delta t C
C Re R*8, In, Re Number C
C OMG R*8, In, Acceleration parameter for SOR C
C NDX I*4, In, First Array Size C
C NDY I*4, In, Second Array Sixe 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/14 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(-1:NDX,-1:NDY), V(-1:NDX,-1:NDY)
DIMENSION A(-1:NDX,-1:NDY,3), B(-1:NDX,-1:NDY)
DIMENSION P(-1:NDX,-1:NDY), W(-1:NDX,-1:NDY,3)
COMMON /MESH/NX,NX1,NX2,NY,NY1,NY2
C Set Value
DYX = hy/hx
DXY = hx/hy
DHX = 1.0/(2.0*hx)
DHY = 1.0/(2.0*hy)
C Set Outer Boundary U,V,P
do i=-1,NX+1
P(i,-1) = P(i,1)
P(i,NY+1) = P(i,NY-1)
U(i,-1) = U(i,0)
V(i,-1) = V(i,0)
U(i,NY+1) = U(i,NY)
V(i,NY+1) = V(i,NY)
end do
do j=-1,NY+1
C P=0 on L3
P(NX,j) = 0.0
P(-1,j) = P(1,j)
P(NX+1,j) = P(NX-1,j)
U(-1,j) = U(0,j)
V(-1,j) = V(0,j)
U(NX+1,j) = U(NX,j)
V(NX+1,j) = V(NX,j)
end do
C Set A and B
BN = 0.0
do j=0,NY
do i=0,NX-1
C Inner
A1 = -DXY
A2 = -DYX
A3 = 2.0*(DXY + DYX)
A4 = -DYX
A5 = -DXY
C on B1
if(i.eq.NX1 .and. j.ge.NY1 .and. j.le.NY2) then
A3 = A3 - DYX
A4 = 0.0
UP0 = 0.0
VP0 = 0.0
else
UP0 = U(i+1,j)
VP0 = V(i+1,j)
end if
C on B2
if(j.eq.NY1 .and. i.ge.NX1 .and. i.le.NX2) then
A3 = A3 - DXY
A5 = 0.0
U0P = 0.0
V0P = 0.0
else
U0P = U(i,j+1)
V0P = V(i,j+1)
end if
C on B3
if(i.eq.NX2 .and. j.ge.NY1 .and. j.le.NY2) then
A2 = 0.0
A3 = A3 - DYX
UM0 = 0.0
VM0 = 0.0
else
UM0 = U(i-1,j)
VM0 = V(i-1,j)
end if
C on B4
if(j.eq.NY2 .and. i.ge.NX1 .and. i.le.NX2) then
A1 = 0.0
A3 = A3 - DXY
U0M = 0.0
V0M = 0.0
else
U0M = U(i,j-1)
V0M = V(i,j-1)
end if
C Set A,B
A(i,j,1) = A1
A(i,j,2) = A2
A(i,j,3) = A3
A(i,j,4) = A4
A(i,j,5) = A5
B(i,j) = ((UP0-UM0)**2*DYX + (V0P-V0M)**2*DXY )/4.0
1 + (VP0-VM0)*(U0P-U0M)/2.0
2 - ((UP0-UM0)*hy + (V0P-V0M)*hx)/(2.0*dt)
C in Square a Pillar
if( (i.gt.NX1.and.i.lt.NX2) .and.
1 (j.gt.NY1.and.j.lt.NY2) ) then
A(i,j,1) = 0.0
A(i,j,2) = 0.0
A(i,j,3) = 1.0
A(i,j,4) = 0.0
A(i,j,5) = 0.0
B(i,j) = P(NX1,NY1)
end if
BN = BN + (B(i,j)/A(i,j,3))**2
end do
end do
C Solve Ax=b
do k=1,ITER
RN = 0.0
do j=0,NY
do i=0,NX-1
R = ( B(i,j) - A(i,j,1)*P(i,j-1) - A(i,j,2)*P(i-1,j)
1 - A(i,j,4)*P(i+1,j) - A(i,j,5)*P(i,j+1) )/A(i,j,3)
2 - P(i,j)
P(i,j) = P(i,j) + OMG*R
RN = RN + (R/A(i,j,3))**2
end do
end do
C Check Conversion
ERR = sqrt(RN/BN)
if(ERR.le.EPS) go to 100
end do
100 CONTINUE
ITER = k
C
RETURN
END
C=================================================================C
SUBROUTINE PHISOL(PHI,A,B,U,V,hx,hy,OMG,NDX,NDY,EPS,ITER)
C=================================================================C
C Solve Ax=b by SOR with 2 dimensional FDM C
C Solve PHI and Compute U,V C
C -div(PHI)=0 and d(PHI)/dx=U, d(PHI)/dy=V C
C-----------------------------------------------------------------C
C PHI(-1:NDX,-1:NDY) R*8, I/O, PHI vector C
C A(-1:NDX,-1:NDY,5) R*8, Wk, Matrix A for -div(PHI)=0 C
C B(-1:NDX,-1:NDY) R*8, Wk, Right Vector for -div(PHI)=0 C
C U(-1:NDX,-1:NDY) R*8, In, U vector C
C V(-1:NDX,-1:NDY) R*8, In, V vector C
C hx R*8, In. Delta x C
C hy R*8, In, Delta y C
C OMG R*8, In, Acceleration parameter for SOR C
C NDX I*4, In, First Array Size C
C NDY I*4, In, Second Array Sixe 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/14 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(-1:NDX,-1:NDY), V(-1:NDX,-1:NDY)
DIMENSION A(-1:NDX,-1:NDY,5), B(-1:NDX,-1:NDY)
DIMENSION PHI(-1:NDX,-1:NDY)
COMMON /MESH/NX,NX1,NX2,NY,NY1,NY2
C Set Value
DYX = hy/hx
DXY = hx/hy
DHX = 1.0/(2.0*hx)
DHY = 1.0/(2.0*hy)
C Set Initial PHI=X
do j=0,NY
do i=0,NX+1
PHI(i,j) = i*hx
end do
end do
C Set A and B
BN = 0.0
do j=1,NY-1
do i=1,NX
C Inner
A1 = -DXY
A2 = -DYX
A3 = 2.0*(DXY + DYX)
A4 = -DYX
A5 = -DXY
BV = 0.0
C on L3
if(i.eq.NX) then
A3 = A3 - DYX
A4 = 0.0
BV = BV + hy
end if
C on B1
if(i.eq.NX1 .and. j.ge.NY1 .and. j.le.NY2) then
A3 = A3 - DYX
A4 = 0.0
end if
C on B2
if(j.eq.NY1 .and. i.ge.NX1 .and. i.le.NX2) then
A3 = A3 - DXY
A5 = 0.0
end if
C on B3
if(i.eq.NX2 .and. j.ge.NY1 .and. j.le.NY2) then
A2 = 0.0
A3 = A3 - DYX
end if
C on B4
if(j.eq.NY2 .and. i.ge.NX1 .and. i.le.NX2) then
A1 = 0.0
A3 = A3 - DXY
end if
C in Square a Pillar
if( (i.gt.NX1.and.i.lt.NX2) .and.
1 (j.gt.NY1.and.j.lt.NY2) ) then
A1 = 0.0
A2 = 0.0
A3 = 1.0
A4 = 0.0
A5 = 0.0
BV = 0.0
end if
C Set A,B
A(i,j,1) = A1
A(i,j,2) = A2
A(i,j,3) = A3
A(i,j,4) = A4
A(i,j,5) = A5
B(i,j) = BV
BN = BN + (BV/A3)**2
end do
end do
C Solve Ax=b
do k=1,ITER
RN = 0.0
do j=1,NY-1
do i=1,NX
R = ( B(i,j) - A(i,j,1)*PHI(i,j-1) - A(i,j,2)*PHI(i-1,j)
1 - A(i,j,4)*PHI(i+1,j) - A(i,j,5)*PHI(i,j+1) )/A(i,j,3)
2 - PHI(i,j)
PHI(i,j) = PHI(i,j) + OMG*R
RN = RN + R**2
end do
end do
C Check Conversion
ERR = sqrt(RN/BN)
if(ERR.le.EPS) go to 100
end do
100 CONTINUE
ITER = k
C Compute U,V
do j=1,NY-1
do i=1,NX
U(i,j) = (PHI(i+1,j) - PHI(i-1,j))*DHX
V(i,j) = (PHI(i,j+1) - PHI(i,j-1))*DHY
end do
end do
C in B
do j=NY1,NY2
do i=NX1,NX2
U(i,j) = 0.0
V(i,j) = 0.0
end do
end do
C on L1 Boundary
do j=1,NY-1
U(0,j) = 1.0
V(0,j) = 0.0
end do
C on L2,L4 Boundary
do i=0,NX
U(i,0) = 1.0
V(i,0) = 0.0
U(i,NY) = 1.0
V(i,NY) = 0.0
end do
C
RETURN
END
C=================================================================C
SUBROUTINE OUTUVP(P,U,V,OMG,hx,hy,NDX,NDY)
C=================================================================C
C Compute PHI and Output PHI,U,V C
C OMG(i,j) = dV/dx - dU/dy C
C-----------------------------------------------------------------C
C P(0:NDX,0:NY) R*8, In, 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 OMG(0:NDX,0:NDY) R*8, Wk, Vorticity Vector C
C hx R*8, In, Delta x C
C hy R*8, In, Delta y C
C NDX I*4, In, First Array Size of P,U,V C
C NDY I*4, In, Second Array Size of P,U,V C
C-----------------------------------------------------------------C
C Written by Yasunori Ushiro, 2003/09/14 C
C ( Hitachi Ltd. and Waseda University ) C
C=================================================================C
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION P(-1:NDX,-1:NDY), OMG(-1:NDX,-1:NDY)
DIMENSION U(-1:NDX,-1:NDY), V(-1:NDX,-1:NDY)
COMMON /MESH/NX,NX1,NX2,NY,NY1,NY2
C Set Omega
do j=0,NY
do i=0,NX
OMG(i,j) = (V(i+1,j)-V(i-1,j))/hx - (U(i,j+1)-U(i,j-1))/hy
end do
end do
C Output
do j=0,NY
Y = j*hy
do i=0,NX
X = i*hx
WRITE(1,100) X,Y,OMG(i,j),U(i,j),V(i,j)
end do
end do
C
RETURN
100 FORMAT(1H ,5F10.4)
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
3. include file（ファイル名:UV.h)
U(i,j) = U(i,j) - dt*( (U(i+1,j)**2 - U(i-1,j)**2)*dx
1 + (U(i,j+1)*V(i,j+1) - U(i,j-1)*V(i,j-1))*dy
2 + (P(i+1,j) - P(i-1,j))*dx
3 + ( (2.0*U(i,j) - U(i-1,j) - U(i+1,j))*ddx
4 + (2.0*U(i,j) - U(i,j-1) - U(i,j+1))*ddy )/Re )
V(i,j) = V(i,j) - dt*( (V(i,j+1)**2 - V(i,j-1)**2)*dy
1 + (U(i+1,j)*V(i+1,j) - U(i-1,j)*V(i-1,j))*dx
2 + (P(i,j+1) - P(i,j-1))*dy
3 + ( (2.0*V(i,j) - V(i-1,j) - V(i+1,j))*ddx
4 + (2.0*V(i,j) - V(i,j-1) - V(i,j+1))*ddy )/Re )
OVER = OVER + abs(U(i,j)) + abs(V(i,j))