戻る

C
C ******************************************************************** C
C C
C F 3 D - F D M C
C C
C Three-Dimensional Fluid Analysis Software C
C C
C by 4th order F.D.M. C
C C
C ( Version 4.0 ) C
C C
C 　Copyright：Yasuhiro MATSUDA 　　　 C
C C
C ******************************************************************** C
C
PROGRAM MAIN
C
IMPLICIT REAL*8( A- H,O-Z )
C
C ***** DATA ( 1 ) *****************************************************
C
*** PARAMETER ( NX = 5, NY = 5, NZ = 5 )
C
PARAMETER ( NX = 11, NY = 11, NZ = 11 )
C
*** PARAMETER ( NX = 21, NY = 21, NZ = 21 )
*** PARAMETER ( NX = 31, NY = 31, NZ = 31 )
C
C **********************************************************************
C
PARAMETER ( NEX= NX-1, NEY= NY-1, NEZ= NZ-1, NNP= NX* NY* NZ )
C
DIMENSION
C
1 U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ), UA( NX,NY,NZ),
2 VA( NX,NY,NZ), WA( NX,NY,NZ), UB( NX,NY,NZ), VB( NX,NY,NZ),
3 WB( NX,NY,NZ), OU( NX,NY,NZ), OV( NX,NY,NZ), OW( NX,NY,NZ),
C
4 P ( NX,NY,NZ), PA( NX,NY,NZ), DP( NX,NY,NZ), OP( NX,NY,NZ),
5 AX( NX,NY,NZ), AY( NX,NY,NZ), AZ( NX,NY,NZ),
6 UN( NNP), VN( NNP), WN( NNP), PN( NNP),
C
7 B ( NX,NY,NZ), HX( NX), HY( NY), HZ( NZ),
8 X ( NX), Y ( NY), Z ( NZ),
9 XL( NEX), YL( NEY), ZL( NEZ)
C
C ----- Attention: Size = Max.( NX,NY,NZ ) + 3 etc. ---
C
DIMENSION FU( NX+3), FV( NX+3), FW( NX+3)
C
COMMON /PAR/ DT, RE, WEI
COMMON /VVP/ VRV( 3000), VRP( 3000)
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
C ***** DATA ( 2 ) *****************************************************
C
*** OPEN ( 11,FILE ='UE4_F3D_FDM_R100_FIG' )
C
OPEN ( 11,FILE ='T_UE10_F3D_FDM_R100_FIG' ) ! Re=100 DT=0.0384
C
*** OPEN ( 11,FILE ='UE20_F3D_FDM_R400_FIG' ) ! Re=400 DT=0.02
*** OPEN ( 11,FILE ='UE30_F3D_FDM_R1000_FIG' ) ! Re=1000 DT=0.01 OK
C
C **********************************************************************
C
CALL GMSH ( XL, YL, ZL, HX, HY, HZ, NEX, NEY, NEZ, MES )
C
555 CONTINUE
C
CALL SINT ( U, V, W, UA, VA, WA, UB, VB, WB, P,
1 PA, DP, OU, OV, OW, OP, B, X, Y, Z,
2 HX, HY, HZ, NX, NY, NZ )
C
LOP = 0
C
C ----------------------------------------------------------
999 LOP = LOP + 1
C
CALL FVLB ( U, V, W, UB, VB, WB, HX, HY, HZ, NX, NY, NZ )
C
C ----- PA ---
C
CALL SORP ( UB, VB, WB, PA, HX, HY, HZ, B, NX, NY, NZ, 1 )
C
C ----- UA, VA, WA ---
C
CALL VLSX ( U, V, W, UA, VA, WA, FU, FV, FW, AX, HX,
1 NX, NY, NZ )
C
CALL VLSY ( U, V, W, UA, VA, WA, FU, FV, FW, AY, HY,
1 NX, NY, NZ )
C
CALL VLSZ ( U, V, W, UA, VA, WA, FU, FV, FW, AZ, HZ,
1 NX, NY, NZ )
C
CALL VLSF ( U, V, W, UA, VA, WA, PA, HX, HY, HZ,
1 NX, NY, NZ )
C
C ----- DP ---
C
CALL SORP ( UA, VA, WA, DP, HX, HY, HZ, B, NX, NY, NZ, 2 )
C
C ----- U, V, W ---
C
CALL PRVL ( U, V, W, UA, VA, WA, OU, OV, OW, P, PA,
1 DP, OP, HX, HY, HZ, NX, NY, NZ )
C
C ----- Stationarity etc. Check ---
C
CALL CHCK ( U, V, W, P, HX, NX, NY, NZ, NNP, *999 )
C
C ----- Ouput for Figures etc.---
C
CALL PLTO ( U, V, W, P, UN, VN, WN, PN, X, Y, Z,
1 NX, NY, NZ, NNP, *555 )
C
CLOSE ( 11)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CHCK ( U, V, W, P, HX, NX, NY, NZ, NNP, * )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
1 P ( NX,NY,NZ), HX( NX)
C
COMMON /CIT/ IT1, IT2
COMMON /PAR/ DT, RE, WEI
COMMON /VVP/ VRV( 3000), VRP( 3000)
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
IF ( MOD( LOP,NPR).EQ.0)
1WRITE(6,2000) LOP, IT1, IT2, EPS*100.D0, PMI, PMX
2000 FORMAT('LOP=',I5,' IT1=',I3,' IT2=',I3,' EPS=',F8.3,'(%)',
1 ' PMI=',F7.3,' PMX=',F7.3)
C
C --------------------------------------------------------
IF ( LOP.EQ.LOPE) THEN
WRITE(6,*) ' '
WRITE(6,*) ' *** Loop End ***'
WRITE(6,2000) LOP, IT1, IT2, EPS*100.D0, PMI, PMX
GO TO 900
END IF
C --------------------------------------------------------
IF ( ITR.GE.( NX-1)* 100 .OR. UMX.GT.1.5D0) THEN
WRITE(6,2200) LOP, IT1, IT2, UMX
2200 FORMAT(' *** Diverged ! * STOP : Loop =',I5,' * IT1=',I3,
1 ' * IT2=',I3,' Umax.=',F11.2 )
STOP
END IF
C --------------------------------------------------------
C
C ----- Stationarity check ---
C
IF ( LOP.EQ.1) RETURN 1
C
IF ( EPS.LE.0.001D0) THEN
C
WRITE(6,*) ' '
WRITE(6,*) ' *** Converged ! *** Loop =', LOP
WRITE(6,*) ' '
WRITE(6,2000) LOP, IT1, IT2, EPS* 100.D0, PMI, PMX
C
C ------------------------------------------------
BMX = 0.D0
DO 100 I = 1, NX
C
BVAR= UMX* DT/ HX( I) ! Max. bx.-value
IF ( BVAR.GE.BMX) BMX = BVAR
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,2100) VRV( LOP), VRP( LOP), BMX
2100 FORMAT(' *** VRV=', F9.4,' * VRP=', F9.4,' * Max.bx=', F9.3)
GO TO 900
C
END IF
C ----------------------------------------------------------
C
RETURN 1
C
900 RETURN
END
C **********************************************************************
C
SUBROUTINE COFS ( HH, NXYZ, L )
C
IMPLICIT REAL*8( A-H,O-Z )
C
DIMENSION HH( NXYZ)
C
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C -------------------------------------------------
C
HA14 = HH( L-2)* ( HH( L-2) + HH( L-1))
1 * ( HH( L-2) + HH( L-1) + HH( L))
2 * ( HH( L-2) + HH( L-1) + HH( L) + HH( L+1))
C
A1 = HH( L-1)* HH( L)* ( HH( L) + HH( L+1))/ HA14
C
A2 = - 2.D0* ( HH( L-1)* ( 2.D0* HH( L) + HH( L+1))
1 - HH( L)* ( HH( L) + HH( L+1)))/ HA14
C
A3 = 6.D0* HH( L)/ HA14
C
A4 = 24.D0* ( HH( L-1)- 2.D0* HH( L)- HH( L+1))/ HA14
C
AG1 = - A1
AG2 = A2/ 2.D0
AG3 = - A3/ 6.D0
AG4 = A4/ 24.D0
C
C ------------------------------------------------
C
HA58 = HH( L-3)*( HH( L-3) + HH( L-2))
1 * ( HH( L-3) + HH( L-2) + HH( L-1))
2 * ( HH( L-3) + HH( L-2) + HH( L-1) + HH( L))
C
A5 = HH( L-2)* HH( L-1)* ( HH( L-1) + HH( L))/ HA58
C
A6 = - 2.D0* ( HH( L-2)* ( 2.D0* HH( L-1) + HH( L))
1 - HH( L-1)* ( HH( L-1) + HH( L)))/ HA58
C
A7 = 6.D0* HH( L-1)/ HA58
C
A8 = 24.D0* ( HH( L-2)-2.D0* HH( L-1)- HH( L))/ HA58
C
C --------------------------------------------------
C
HA912 = HH( L-4)* ( HH( L-4) + HH( L-3))
1 * ( HH( L-4) + HH( L-3) + HH( L-2))
2 * ( HH( L-4) + HH( L-3) + HH( L-2) + HH( L-1))
C
A9 = HH( L-3)* HH( L-2)* ( HH( L-2) + HH( L-1))/ HA912
C
A10 = - 2.D0* ( HH( L-3)* ( 2.D0* HH( L-2) + HH( L-1))
1 - HH( L-2)* ( HH( L-2) + HH( L-1)))
2 / HA912
C
A11 = 6.D0* HH( L-2)/ HA912
C
A12 = 24.D0* ( HH( L-3) - 2.D0* HH( L-2) - HH( L-1))
1 / HA912
C
C ----------------------------------------------------------
C
HA1316 = HH( L-1)* ( HH( L-1) + HH( L))
1 * ( HH( L-1) + HH( L) + HH( L+1))
2 * ( HH( L-1) + HH( L) + HH( L+1) + HH( L+2))
C
A13 = HH( L)* HH( L+1)* ( HH( L+1) + HH( L+2))/ HA1316
C
A14 = - 2.D0* ( HH( L)* ( 2.D0* HH( L+1) + HH( L+2))
1 - HH( L+1)* ( HH( L+1) + HH( L+2)))
2 / HA1316
C
A15 = 6.D0* HH( L+1)/ HA1316
C
A16 = 24.D0* ( HH( L) - 2.D0* HH( L+1) - HH( L+2))/ HA1316
C
C ----------------------------------------------------------
C
HB14 = HH( L-1)* ( HH( L-1) + HH( L))
1 * ( HH( L-1) + HH( L) + HH( L+1))
C
HB141 = ( HH( L-2) + HH( L-1))
1 * ( HH( L-2) + HH( L-1) + HH( L))
2 * ( HH( L-2) + HH( L-1) + HH( L) + HH( L+1))
C
B1 = ( - A1* HB141 - HH( L)* ( HH( L) + HH( L+1)))/ HB14
B2 = ( - A2* HB141 + 2.D0* ( 2.D0* HH( L) + HH( L+1)))/ HB14
B3 = ( - A3* HB141 - 6.D0)/ HB14
B4 = ( - A4* HB141 )/ HB14
C
BG1 = - ( A13 + B1)
BG2 = ( A14 + B2)/ 2.D0
BG3 = - ( A15 + B3)/ 6.D0
BG4 = ( A16 + B4)/ 24.D0
C
C ----------------------------------------------------------
C
HB58 = HH( L-2)* ( HH( L-2) + HH( L-1))
1 * ( HH( L-2) + HH( L-1) + HH( L))
C
HB581 = ( HH( L-3) + HH( L-2))
1 * ( HH( L-3) + HH( L-2) + HH( L-1))
2 * ( HH( L-3) + HH( L-2) + HH( L-1) + HH( L))
C
B5 = ( - A5* HB581 - HH( L-1)* ( HH( L-1) + HH( L)))/ HB58
B6 = ( - A6* HB581 + 2.D0* ( 2.D0* HH( L-1) + HH( L)))/ HB58
B7 = ( - A7* HB581 - 6.D0)/ HB58
B8 = ( - A8* HB581 )/ HB58
C
C ----------------------------------------------------------
C
HB912 = HH( L-3)* ( HH( L-3) + HH( L-2))
1 * ( HH( L-3) + HH( L-2) + HH( L-1))
C
HB9121 = ( HH( L-4) + HH( L-3))
1 * ( HH( L-4) + HH( L-3) + HH( L-2))
2 * ( HH( L-4) + HH( L-3) + HH( L-2) + HH( L-1))
C
B9 = ( - A9* HB9121 - HH(L-2)* ( HH( L-2) + HH( L-1)))/ HB912
C
B10 = ( - A10* HB9121 + 2.D0* ( 2.D0* HH( L-2) + HH( L-1)))
1 / HB912
C
B11 = ( - A11* HB9121 - 6.D0)/ HB912
C
B12 = ( - A12* HB9121 )/ HB912
C
C ----------------------------------------------------------
C
HD14 = HH( L)* HH( L+1)
C
HD141 = ( HH( L-2) + HH( L-1))
1 * ( HH( L-2) + HH( L-1) + HH( L) + HH( L+1))
C
HD142 = HH( L-1)* ( HH( L-1) + HH( L) + HH( L+1))
C
D1 = ( A1* HD141 + B1* HD142 + HH( L) + HH( L+1))/ HD14
D2 = ( A2* HD141 + B2* HD142 - 2.D0)/ HD14
D3 = ( A3* HD141 + B3* HD142 )/ HD14
D4 = ( A4* HD141 + B4* HD142 )/ HD14
C
C ----------------------------------------------------------
C
HD58 = HH( L-1)* HH( L)
C
HD581 = ( HH( L-3) + HH( L-2))
1 * ( HH( L-3) + HH( L-2) + HH( L-1) + HH( L))
C
HD582 = HH( L-2)* ( HH( L-2) + HH( L-1) + HH( L))
C
D5 = ( A5* HD581 + B5* HD582 + HH( L-1) + HH( L))/ HD58
D6 = ( A6* HD581 + B6* HD582 - 2.D0)/ HD58
D7 = ( A7* HD581 + B7* HD582 )/ HD58
D8 = ( A8* HD581 + B8* HD582 )/ HD58
C
C ----------------------------------------------------------
C
HD912 = HH( L-2)* HH( L-1)
C
HD9121 = ( HH( L-4) + HH( L-3))
1 * ( HH( L-4) + HH( L-3) + HH( L-2) + HH( L-1))
C
HD9122 = HH( L-3)* ( HH( L-3) + HH( L-2) + HH( L-1))
C
D9 = ( A9* HD9121 + B9* HD9122 + HH( L-2) + HH( L-1))/ HD912
D10 = ( A10* HD9121 + B10* HD9122 - 2.D0)/ HD912
D11 = ( A11* HD9121 + B11* HD9122 )/ HD912
D12 = ( A12* HD9121 + B12* HD9122 )/ HD912
C
C ----------------------------------------------------------
C
HE14 = HH( L) + HH( L+1)
C
E1 = ( A1* ( HH( L-2) + HH( L-1))
1 + B1* HH( L-1) - D1* HH( L) + 1.D0)/ HE14
C
E2 = ( A2* ( HH( L-2) + HH( L-1))
1 + B2* HH( L-1) - D2* HH( L))/ HE14
C
E3 = ( A3* ( HH( L-2) + HH( L-1))
1 + B3* HH( L-1) - D3* HH( L))/ HE14
C
E4 = ( A4* ( HH( L-2) + HH( L-1))
1 + B4* HH( L-1) - D4* HH( L))/ HE14
C
C ----------------------------------------------------------
C
HE58 = HH( L-1) + HH( L)
C
E5 = ( A5* ( HH( L-3) + HH( L-2))
1 + B5* HH( L-2) - D5* HH( L-1)+ 1.D0)/ HE58
C
E6 = ( A6* ( HH( L-3) + HH( L-2))
1 + B6* HH( L-2) - D6* HH( L-1))/ HE58
C
E7 = ( A7* ( HH( L-3) + HH( L-2))
1 + B7* HH( L-2) - D7* HH( L-1))/ HE58
C
E8 = ( A8* ( HH( L-3) + HH( L-2))
1 + B8* HH( L-2) - D8* HH( L-1))/ HE58
C
DG1 = E5
DG2 = - E6/ 2.D0
DG3 = E7/ 6.D0
DG4 = - E8/ 24.D0
C
C ----------------------------------------------------------
C
HE912 = HH( L-2) + HH( L-1)
C
E9 = ( A9* ( HH( L-4) + HH( L-3))
1 + B9* HH( L-3) - D9* HH( L-2)+ 1.D0)/ HE912
C
E10 = ( A10* ( HH( L-4) + HH( L-3))
1 + B10* HH( L-3) - D10* HH( L-2))/ HE912
C
E11 = ( A11* ( HH( L-4) + HH( L-3))
1 + B11* HH( L-3) - D11* HH( L-2))/ HE912
C
E12 = ( A12* ( HH( L-4) + HH( L-3))
1 + B12* HH( L-3) - D12* HH( L-2))/ HE912
C
CG1 = D5 + E9
CG2 = - ( D6 + E10)/ 2.D0
CG3 = ( D7 + E11)/ 6.D0
CG4 = - ( D8 + E12)/ 24.D0
C
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FCFX ( FU, FV, FW, U, V, W, AX, NX, NY, NZ,
1 M, I, J, K )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION FU( NX+3), FV( NX+3), FW( NX+3),
1 U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
2 AX( NX,NY,NZ)
C
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C ----------------------------
C
FU( M) = ! FU( I-1/2) & FU( I+1/2)
C
1 AX( I,J,K) * ( AG1* U( M-2,J,K) + BG1* U( M-1,J,K)
2 + CG1* U( M,J,K) + DG1* U( M+1,J,K))
C
3 + AX( I,J,K)** 2 * ( AG2* U( M-2,J,K) + BG2* U( M-1,J,K)
4 + CG2* U( M,J,K) + DG2* U( M+1,J,K))
C
5 + AX( I,J,K)** 3 * ( AG3* U( M-2,J,K) + BG3* U( M-1,J,K)
6 + CG3* U( M,J,K) + DG3* U( M+1,J,K))
C
7 + AX( I,J,K)** 4 * ( AG4* U( M-2,J,K) + BG4* U( M-1,J,K)
8 + CG4* U( M,J,K) + DG4* U( M+1,J,K))
C
C ----------------------------
C
FV( M) = ! FV( I-1/2) & FV( I+1/2)
C
1 AX( I,J,K) * ( AG1* V( M-2,J,K) + BG1* V( M-1,J,K)
2 + CG1* V( M,J,K) + DG1* V( M+1,J,K))
C
3 + AX( I,J,K)** 2 * ( AG2* V( M-2,J,K) + BG2* V( M-1,J,K)
4 + CG2* V( M,J,K) + DG2* V( M+1,J,K))
C
5 + AX( I,J,K)** 3 * ( AG3* V( M-2,J,K) + BG3* V( M-1,J,K)
6 + CG3* V( M,J,K) + DG3* V( M+1,J,K))
C
7 + AX( I,J,K)** 4 * ( AG4*V( M-2,J,K) + BG4* V( M-1,J,K)
8 + CG4*V( M,J,K) + DG4* V( M+1,J,K))
C
C ----------------------------
C
FW( M) = ! FW( I-1/2) & FW( I+1/2)
C
1 AX( I,J,K) * ( AG1* W( M-2,J,K) + BG1* W( M-1,J,K)
2 + CG1* W( M,J,K) + DG1* W( M+1,J,K))
C
3 + AX( I,J,K)** 2 * ( AG2* W( M-2,J,K) + BG2* W( M-1,J,K)
4 + CG2* W( M,J,K) + DG2* W( M+1,J,K))
C
5 + AX( I,J,K)** 3 * ( AG3* W( M-2,J,K) + BG3* W( M-1,J,K)
6 + CG3* W( M,J,K) + DG3* W( M+1,J,K))
C
7 + AX( I,J,K)** 4 * ( AG4* W( M-2,J,K) + BG4* W( M-1,J,K)
8 + CG4* W( M,J,K) + DG4* W( M+1,J,K))
C
C ----------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FCFY ( FU, FV, FW, U, V, W, AY, NX, NY, NZ,
1 M, I, J, K )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION FU( NY+3), FV( NY+3), FW( NY+3),
1 U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
2 AY( NX,NY,NZ)
C
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C ----------------------------
C
FU( M) = ! FU( I-1/2) & FU( I+1/2)
C
1 AY( I,J,K) * ( AG1* U( I,M-2,K) + BG1* U( I,M-1,K)
2 + CG1* U( I,M,K) + DG1* U( I,M+1,K))
C
3 + AY( I,J,K)**2 * ( AG2* U( I,M-2,K) + BG2* U( I,M-1,K)
4 + CG2* U( I,M,K) + DG2* U( I,M+1,K))
C
5 + AY( I,J,K)**3 * ( AG3* U( I,M-2,K) + BG3* U( I,M-1,K)
6 + CG3* U( I,M,K) + DG3* U( I,M+1,K))
C
7 + AY( I,J,K)**4 * ( AG4* U( I,M-2,K) + BG4* U( I,M-1,K)
8 + CG4* U( I,M,K) + DG4* U( I,M+1,K))
C
C ----------------------------
C
FV( M) = ! FV( I-1/2) & FV( I+1/2)
C
1 AY( I,J,K) * ( AG1* V( I,M-2,K) + BG1* V( I,M-1,K)
2 + CG1* V( I,M,K) + DG1* V( I,M+1,K))
C
3 + AY( I,J,K)**2 * ( AG2* V( I,M-2,K) + BG2* V( I,M-1,K)
4 + CG2* V( I,M,K) + DG2* V( I,M+1,K))
C
5 + AY( I,J,K)**3 * ( AG3* V( I,M-2,K) + BG3* V( I,M-1,K)
6 + CG3* V( I,M,K) + DG3* V( I,M+1,K))
C
7 + AY( I,J,K)**4 * ( AG4* V( I,M-2,K) + BG4* V( I,M-1,K)
8 + CG4* V( I,M,K) + DG4* V( I,M+1,K))
C
C ----------------------------
C
FW( M) = ! FW( I-1/2) & FW( I+1/2)
C
1 AY( I,J,K) * ( AG1* W( I,M-2,K) + BG1* W( I,M-1,K)
2 + CG1* W( I,M,K) + DG1* W( I,M+1,K))
C
3 + AY( I,J,K)**2 * ( AG2* W( I,M-2,K) + BG2* W( I,M-1,K)
4 + CG2* W( I,M,K) + DG2* W( I,M+1,K))
C
5 + AY( I,J,K)**3 * ( AG3* W( I,M-2,K) + BG3* W( I,M-1,K)
6 + CG3* W( I,M,K) + DG3* W( I,M+1,K))
C
7 + AY( I,J,K)**4 * ( AG4* W( I,M-2,K) + BG4* W( I,M-1,K)
8 + CG4* W( I,M,K) + DG4* W( I,M+1,K))
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FCFZ ( FU, FV, FW, U, V, W, AZ, NX, NY, NZ,
1 M, I, J, K )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION FU( NZ+3), FV( NZ+3), FW( NZ+3),
1 U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
2 AZ( NX,NY,NZ)
C
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C --------------------------------------------------------------------
C
FU( M) = ! FU( I-1/2) & FU( I+1/2)
C
1 AZ( I,J,K) *( AG1* U( I,J,M-2) + BG1* U( I,J,M-1)
2 + CG1* U( I,J,M ) + DG1* U( I,J,M+1))
C
3 + AZ( I,J,K)** 2 *( AG2* U( I,J,M-2) + BG2* U( I,J,M-1)
4 + CG2* U( I,J,M ) + DG2* U( I,J,M+1))
C
5 + AZ( I,J,K)** 3 *( AG3* U( I,J,M-2) + BG3* U( I,J,M-1)
6 + CG3* U( I,J,M ) + DG3* U( I,J,M+1))
C
7 + AZ( I,J,K)** 4 *( AG4* U( I,J,M-2) + BG4* U( I,J,M-1)
8 + CG4* U( I,J,M ) + DG4* U( I,J,M+1))
C
C --------------------------------------------------------------------
C
FV( M) = ! FV( I-1/2) & FV( I+1/2)
C
1 AZ( I,J,K) *( AG1* V( I,J,M-2) + BG1* V( I,J,M-1)
2 + CG1* V( I,J,M ) + DG1* V( I,J,M+1))
C
3 + AZ( I,J,K)** 2 *( AG2* V( I,J,M-2) + BG2* V( I,J,M-1)
4 + CG2* V( I,J,M ) + DG2* V( I,J,M+1))
C
5 + AZ( I,J,K)** 3 *( AG3* V( I,J,M-2) + BG3* V( I,J,M-1)
6 + CG3* V( I,J,M ) + DG3* V( I,J,M+1))
C
7 + AZ( I,J,K)** 4 *( AG4* V( I,J,M-2) + BG4* V( I,J,M-1)
8 + CG4* V( I,J,M ) + DG4* V( I,J,M+1))
C
C --------------------------------------------------------------------
C
FW( M) = ! FW( I-1/2) & FW( I+1/2)
C
1 AZ( I,J,K) *( AG1* W( I,J,M-2) + BG1* W( I,J,M-1)
2 + CG1* W( I,J,M ) + DG1* W( I,J,M+1))
C
3 + AZ( I,J,K)** 2 *( AG2* W( I,J,M-2) + BG2* W( I,J,M-1)
4 + CG2* W( I,J,M ) + DG2* W( I,J,M+1))
C
5 + AZ( I,J,K)** 3 *( AG3* W( I,J,M-2) + BG3* W( I,J,M-1)
6 + CG3* W( I,J,M ) + DG3* W( I,J,M+1))
C
7 + AZ( I,J,K)** 4 *( AG4* W( I,J,M-2) + BG4* W( I,J,M-1)
8 + CG4* W( I,J,M ) + DG4* W( I,J,M+1))
C
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE FVLB ( U, V, W, UB, VB, WB, HX, HY, HZ, NX,NY,NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), UB( NX,NY,NZ), V ( NX,NY,NZ),
1 VB( NX,NY,NZ), W ( NX,NY,NZ), WB( NX,NY,NZ),
2 HX( NX), HY( NY), HZ( NZ )
C
COMMON /PAR/ DT, RE, WEI
C
C ----------------------------------------------------------
DO 100 K = 2, NZ-1
C
HZZ = HZ( K-1) + HZ( K)
C ----------------------------------------------------------
DO 100 J = 2, NY-1
C
HYY = HY( J-1) + HY( J)
C ----------------------------------------------------------
DO 100 I = 2, NX-1
C
HXX = HX( I-1) + HX( I)
C
E1 = ( 2.D0 / HX( I-1)/ HXX )
E2 = ( - 2.D0 / HX( I-1)/ HX( I))
E3 = ( 2.D0 / HX( I) / HXX )
C
E4 = ( 2.D0 / HY( J-1)/ HYY )
E5 = ( - 2.D0 / HY( J-1)/ HY( J))
E6 = ( 2.D0 / HY( J) / HYY )
C
E7 = ( 2.D0 / HZ( K-1)/ HZZ )
E8 = ( - 2.D0 / HZ( K-1)/ HZ( K))
E9 = ( 2.D0 / HZ( K) / HZZ )
C
C ------------------------------------------------------
F1 = ( - HX( I)/ HX( I-1) / HXX )
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I) / HXX )
C
F4 = ( - HY( J)/ HY( J-1) / HYY )
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J) / HYY )
C
F7 = ( - HZ( K)/ HZ( K-1) / HZZ )
F8 = (( HZ( K) - HZ( K-1))/ HZ( K-1)/ HZ( K))
F9 = ( HZ( K-1)/ HZ( K) / HZZ )
C
C --------------------------------------
UB( I,J,K) = U( I,J,K)
C
1 + DT* (( E1* U( I-1,J,K) + E2* U( I,J,K) + E3* U( I+1,J,K)
2 + E4* U( I,J-1,K) + E5* U( I,J,K) + E6* U( I,J+1,K)
3 + E7* U( I,J,K-1) + E8* U( I,J,K) + E9* U( I,J,K+1))
4 / RE
C
5 - U( I,J,K)*( F1* U( I-1,J,K) + F2* U( I,J,K) + F3* U( I+1,J,K))
6 - V( I,J,K)*( F4* U( I,J-1,K) + F5* U( I,J,K) + F6* U( I,J+1,K))
7 - W( I,J,K)*( F7* U( I,J,K-1) + F8* U( I,J,K) + F9* U( I,J,K+1)))
C
C --------------------------------------
VB( I,J,K) = V( I,J,K)
C
1 + DT* (( E1* V( I-1,J,K) + E2* V( I,J,K) + E3* V( I+1,J,K)
2 + E4* V( I,J-1,K) + E5* V( I,J,K) + E6* V( I,J+1,K)
3 + E7* V( I,J,K-1) + E8* V( I,J,K) + E9* V( I,J,K+1))
4 / RE
C
5 - U( I,J,K)*( F1* V( I-1,J,K) + F2* V( I,J,K) + F3* V( I+1,J,K))
6 - V( I,J,K)*( F4* V( I,J-1,K) + F5* V( I,J,K) + F6* V( I,J+1,K))
7 - W( I,J,K)*( F7* V( I,J,K-1) + F8* V( I,J,K) + F9* V( I,J,K+1)))
C
C --------------------------------------
WB( I,J,K) = W( I,J,K)
C
1 + DT* (( E1* W( I-1,J,K) + E2* W( I,J,K) + E3* W( I+1,J,K)
2 + E4* W( I,J-1,K) + E5* W( I,J,K) + E6* W( I,J+1,K)
3 + E7* W( I,J,K-1) + E8* W( I,J,K) + E9* W( I,J,K+1))
4 / RE
C
5 - U( I,J,K)*( F1* W( I-1,J,K) + F2* W( I,J,K) + F3* W( I+1,J,K))
6 - V( I,J,K)*( F4* W( I,J-1,K) + F5* W( I,J,K) + F6* W( I,J+1,K))
7 - W( I,J,K)*( F7* W( I,J,K-1) + F8* W( I,J,K) + F9* W( I,J,K+1)))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE GMSH ( XL, YL, ZL, HX, HY, HZ, LX, LY, LZ, MES )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XL( LX), YL( LY), ZL( LZ), HX( LX+1), HY( LY+1),
1 HZ( LZ+1)
C
C ------------------------------------------------
C
CALL TTLE
C
C ------------------------------------------------
C
C MES = 0 : Meshes of equal length ----- EM
C
C MES = 1 : JSME-Paper Version --------- UEM
C
C ------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,*) ' ** MES = 1 ( Meshes of unequal length ) or 0 ( Meshes
1of equal length ) ?'
** READ(5,*) MES
MES = 1
WRITE(6,*) ' ** MES =', MES
WRITE(6,*) ' '
C
C --------------------------------------
IF ( MES.EQ.0) GO TO 333
IF ( MES.EQ.1) GO TO 250
C --------------------------------------
C
333 WRITE(6,*) '* Meshes of equal length '
C
C ------------------------------------------------
DO 100 I = 1, LX
C
XL( I) = DFLOAT( I-1)/ DFLOAT( LX-1)
C
100 CONTINUE
C ------------------------------------------------
DO 200 J = 1, LY
C
YL( J) = DFLOAT( J-1)/ DFLOAT( LY-1)
C
200 CONTINUE
C ------------------------------------------------
DO 300 K = 1, LZ
C
ZL( K) = DFLOAT( K-1)/ DFLOAT( LZ-1)
C
300 CONTINUE
C ------------------------------------------------
GO TO 990
C
C ----------------------------------------------------------
250 WRITE(6,*) '*** Meshes of unequal length ( UEM ) ***'
C
C ---- JSME-Paper Ref.
C
AA = 2.D0
DX = 2.D0/ DFLOAT( LX-1)
DY = 2.D0/ DFLOAT( LY-1)
DZ = 2.D0/ DFLOAT( LZ-1)
C
C ----- Important ----
C
MDX = IDNINT( DFLOAT( LX+1)/ 2.D0)
MX1 = MDX - 1
C
MDY = IDNINT( DFLOAT( LY+1)/ 2.D0)
MY1 = MDY - 1
C
MDZ = IDNINT( DFLOAT( LZ+1)/ 2.D0)
MZ1 = MDZ - 1
C ---------------------
C
XL( MDX) = 0.5D0
YL( MDY) = 0.5D0
ZL( MDZ) = 0.5D0
C
C ----------------------------------------------------
DO 400 I = 1, MX1
C
XL( I) = ( DEXP( AA* DFLOAT( I-1)* DX) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
XL( LX+1-I ) = 1.D0 - XL( I)
C
400 CONTINUE
C ----------------------------------------------------
DO 500 J = 1, MY1
C
YL( J) = ( DEXP( AA* DFLOAT( J-1)* DY) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
YL( LY+1-J) = 1.D0 - YL( J)
C
500 CONTINUE
C ----------------------------------------------------
DO 600 K = 1, MZ1
C
ZL( K) = ( DEXP( AA* DFLOAT( K-1)* DZ) - 1.D0)
1 /( 2.D0* ( DEXP( AA) - 1.D0))
ZL( LZ+1-K) = 1.D0 - ZL( K)
C
600 CONTINUE
C -----------------------------------------------------------
C
990 CONTINUE
C -----------------------------------------------------------
DO 700 I = 1, LX-1
C
HX( I) = DABS( XL( I+1) - XL( I))
C
700 CONTINUE
C ------------------------------------------------
DO 800 J = 1, LY-1
C
HY( J) = DABS( YL( J+1) - YL( J))
C
800 CONTINUE
C ------------------------------------------------
DO 900 J = 1, LZ-1
C
HZ( J) = DABS( ZL( J+1) - ZL( J))
C
900 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTO ( U, V, W, P, UN, VN, WN, PN,
1 X, Y, Z, NX, NY, NZ, NNP, * )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U( NX,NY,NZ), V( NX,NY,NZ), W( NX,NY,NZ), P( NX,NY,NZ),
1 UN( NNP), VN( NNP), WN( NNP), PN( NNP),
2 X ( NX), Y ( NY), Z ( NZ )
C
COMMON /PAR/ DT, RE, WEI
COMMON /VVP/ VRV( 3000), VRP( 3000)
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
C -------------
J = 1
C ------------------------------------------------
DO 100 IZ = 1, NZ
DO 100 IY = 1, NY
DO 100 IX = 1, NX
C
UN( J) = U( IX, IY, IZ)
VN( J) = V( IX, IY, IZ)
WN( J) = W( IX, IY, IZ)
PN( J) = P( IX, IY, IZ)
J = J + 1
C
100 CONTINUE
C ------------------------------------------------
C
NEX = NX-1
C
WRITE(11,*) NX, NY, NZ
WRITE(11,*) NEX, NNP, LOP
WRITE(6,2000) NEX, NNP, RE, DT, LOP
2000 FORMAT(/,' *** NEX =', I3,' ( NNP =',I6,') * Re =',F7.1,' DT =',
1 F7.4,' LOP =',I4)
C
WRITE(11,*) ( X ( I), I = 1, NX )
WRITE(11,*) ( Y ( I), I = 1, NY )
WRITE(11,*) ( Z ( I), I = 1, NZ )
WRITE(11,*) ( UN( I), I = 1, NNP )
WRITE(11,*) ( VN( I), I = 1, NNP )
WRITE(11,*) ( WN( I), I = 1, NNP )
WRITE(11,*) ( PN( I), I = 1, NNP )
WRITE(11,*) ( VRV( I), I = 1, LOP )
WRITE(11,*) ( VRP( I), I = 1, LOP )
C
IF ( LOP.LT.10) THEN
NLAST = LOP
ELSE
NLAST = ( LOP/10)* 10
END IF
NLAST = LOP
WRITE(6,*) ' '
WRITE(6,*) ' *** Computation End ***'
C
NYN = 0
IF ( NYN.EQ.1) RETURN 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRVL ( U, V, W, UA, VA, WA, OU, OV, OW, P,
1 PA, DP, OP, HX, HY, HZ, NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
1 UA( NX,NY,NZ), VA( NX,NY,NZ), WA( NX,NY,NZ),
2 P ( NX,NY,NZ), PA( NX,NY,NZ), DP( NX,NY,NZ),
3 OP( NX,NY,NZ),
4 OU( NX,NY,NZ), OV( NX,NY,NZ), OW( NX,NY,NZ),
5 HX( NX), HY( NY), HZ( NZ )
C
COMMON /PAR/ DT,RE,WEI
C
C ------------------------------------------------
DO 100 K = 2, NZ-1
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C
OP( I,J,K) = P( I,J,K)
C
100 CONTINUE
C
C ----- Pressure ---------------------------------
C
DO 200 K = 1, NZ
DO 200 J = 1, NY
DO 200 I = 1, NX
C
P( I,J,K) = PA( I,J,K) + DP( I,J,K)
C
200 CONTINUE
C ------------------------------------------------
C
C ----- Expression (1.4) ---
C
C ------------------------------------------------
DO 300 K = 2, NZ-1
C
HZZ = HZ( K-1) + HZ( K)
C
C ------------------------------------------------
DO 300 J = 2, NY-1
C
HYY = HY( J-1) + HY( J)
C
C ------------------------------------------------
DO 300 I = 2, NX-1
C
HXX = HX( I-1) + HX( I)
C
OU( I,J,K) = U( I,J,K)
OV( I,J,K) = V( I,J,K)
OW( I,J,K) = W( I,J,K)
C
F1 = ( - HX( I)/ HX( I-1) / HXX )
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I) / HXX )
C
F4 = ( - HY( J)/ HY( J-1) / HYY )
F5 = (( HY( J) - HY( J-1))/ HY( J-1)/ HY( J))
F6 = ( HY( J-1)/ HY( J) / HYY )
C
F7 = ( - HZ( K)/ HZ( K-1) / HZZ )
F8 = (( HZ( K) - HZ( K-1))/ HZ( K-1)/ HZ( K))
F9 = ( HZ( K-1)/ HZ( K) / HZZ )
C
C --------------------------------------
C
U( I,J,K) = UA( I,J,K)
1 - DT* ( F1* DP( I-1,J,K) + F2* DP( I,J,K) + F3* DP( I+1,J,K))
C --------------------------------------
C
V( I,J,K) = VA( I,J,K)
1 - DT* ( F4* DP( I,J-1,K) + F5* DP( I,J,K) + F6* DP( I,J+1,K))
C
C --------------------------------------
C
W( I,J,K) = WA( I,J,K)
1 - DT* ( F7* DP( I,J,K-1) + F8* DP( I,J,K) + F9* DP( I,J,K+1))
C --------------------------------------
C
300 CONTINUE
C -----------------------------------------------------------
C
CALL VARC ( U, V, W, P, OU, OV, OW, OP, NX, NY, NZ )
C
C -----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SINT ( U, V, W, UA, VA, WA, UB, VB, WB,
1 P, PA, DP, OU, OV, OW, OP, B, X,
2 Y, Z, HX, HY, HZ, NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
1 UA( NX,NY,NZ), VA( NX,NY,NZ), WA( NX,NY,NZ),
2 UB( NX,NY,NZ), VB( NX,NY,NZ), WB( NX,NY,NZ),
3 OU( NX,NY,NZ), OV( NX,NY,NZ), OW( NX,NY,NZ),
C
4 P ( NX,NY,NZ), PA( NX,NY,NZ), DP( NX,NY,NZ),
5 OP( NX,NY,NZ), B ( NX,NY,NZ),
6 X ( NX), Y ( NY), Z ( NZ),
7 HX( NX), HY( NY), HZ( NZ )
C
COMMON /PAR/ DT, RE, WEI
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
WRITE (6,2000)
2000 FORMAT(/' * Reynolds number = ? ')
** READ ( 5,*) RE
RE = 100.D0
C ================
C
WRITE (6,2100)
2100 FORMAT(/' * Time step = ? ')
** READ ( 5,*) DT
DT = 0.0384D0
C ==================
C
C ----- Weighting Parameter ---
C
WRITE (6,2200)
2200 FORMAT(/' * WEI = ?')
C ***** Attention ***
C
WEI = 1.D0
C ===============
C
WRITE (6,2300)
2300 FORMAT(/' * Timing of print output = ? ')
** READ ( 5,*) NPR
NPR = 10
C ==============
C
WRITE (6,2400)
2400 FORMAT(/' * How many Loops = ? ')
** READ ( 5,*) LOPE
LOPE = 200
WRITE(6,2500) LOPE
2500 FORMAT(' ** LOPE =',I5,/)
C
C ----- Initialization ---
C
NEX = NX-1
C
MNX = ( NX+1)/ 2
MNY = ( NY+1)/ 2
MNZ = ( NZ+1)/ 2
C
C --------------------------------------
DO 100 I = 1, MNX-1
C
X( I) = ( DEXP( 2.D0* DFLOAT( I-1)* 2.D0/ DFLOAT( NEX))- 1.D0 )
1 /( 2.D0* ( DEXP( 2.D0) - 1.D0))
C
100 CONTINUE
X( MNX) = 0.5D0
C --------------------------------------
DO 200 I = MNX+1, NEX
C
X( I) = 1.D0 - X( NX+1-I)
C
200 CONTINUE
X( NX) = 1.D0
C --------------------------------------
DO 300 I = 1, NX
C
Y( I) = X( I)
Z( I) = X( I)
C
300 CONTINUE
C --------------------------------------
DO 400 I = 1, NX-1
C
HX( I) = X( I+1) - X( I)
HY( I) = Y( I+1) - Y( I)
HZ( I) = Z( I+1) - Z( I)
C
400 CONTINUE
C --------------------------------------
HX( NX) = HX( NX-1)
HY( NY) = HY( NY-1)
HZ( NZ) = HZ( NZ-1)
C --------------------------------------
C
C ----- Initial value ---
C
C --------------------------------------
DO 500 K = 1, NZ
DO 500 J = 1, NY
DO 500 I = 1, NX
C
U ( I,J,K) = 0.D0
UA( I,J,K) = 0.D0
UB( I,J,K) = 0.D0
OU( I,J,K) = 0.D0
C
V ( I,J,K) = 0.D0
VA( I,J,K) = 0.D0
VB( I,J,K) = 0.D0
OV( I,J,K) = 0.D0
C
W ( I,J,K) = 0.D0
WA( I,J,K) = 0.D0
WB( I,J,K) = 0.D0
OW( I,J,K) = 0.D0
C
B ( I,J,K) = 0.D0
C
P ( I,J,K) = 0.D0
PA( I,J,K) = 0.D0
DP( I,J,K) = 0.D0
OP( I,J,K) = 0.D0
C
500 CONTINUE
C --------------------------------------
C
C ----- Moving wall boundary ---
C
C --------------------------------------
DO 600 K = 1, NZ
DO 600 I = 1, NX
C
U ( I,NY,K) = 1.D0
UA( I,NY,K) = 1.D0
UB( I,NY,K) = 1.D0
C
600 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SORP ( SU, SV, SW, SP, HX, HY, HZ, B, NX, NY, NZ, KY )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION SU( NX,NY,NZ), SV( NX,NY,NZ), SW( NX,NY,NZ),
1 SP( NX,NY,NZ), B ( NX,NY,NZ),
2 HX( NX), HY( NY), HZ( NZ)
C
DIMENSION EP( 6), ! EP: Left hand side EU: Right --
1 EU( 6)
C
COMMON /CIT/ IT1, IT2
COMMON /PAR/ DT, RE, WEI
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
C ----- To solve Expression (1.7) or Expression (1.17) ---
C
C --- SP(I,J,K): P* or Del_P ---
C
C AP = 1.D0, 1.4D0, 1.5D0, 1.6D0
C
AP = 1.5D0 ! Attention ( Acceleration Par.)
C ====================================================
C
ITR = 0
ORM = 0.D0
C
C ----- Boundary conditions ----------------------
C
DO 100 I = 2, NX-1
C
SP( I, 1, 1 ) = SP( I, 2, 2 )
SP( I, 1, NZ ) = SP( I, 2, NZ-1 )
SP( I, NY, NZ ) = SP( I, NY-1, NZ-1 )
SP( I, NY, 1 ) = SP( I, NY-1, 2 )
C
100 CONTINUE
C ------------------------------------------------
DO 200 J = 2, NY-1
C
SP( 1, J, 1 ) = SP( 2, J, 2)
SP( NX, J, 1 ) = SP( NX-1, J, 2)
SP( NX, J, NZ ) = SP( NX-1, J, NZ-1)
SP( 1, J, NZ ) = SP( 2, J, NZ-1)
C
200 CONTINUE
C ------------------------------------------------
DO 300 K = 2, NZ-1
C
SP( 1, 1, K) = SP( 2, 2, K)
SP( 1, NY, K) = SP( 2, NY-1, K)
SP( NX, NY, K) = SP( NX-1, NY-1, K)
SP( NX, 1, K) = SP( NX-1, 2, K)
C
300 CONTINUE
C -------------------------------------------
C
SP( 1, 1, 1) = SP( 2, 2, 2)
SP( NX, 1, 1) = SP( NX-1,2, 2)
SP( NX, 1, NZ) = SP( NX-1,2, NZ-1)
SP( 1, 1, NZ) = SP( 2, 2, NZ-1)
C
SP( 1, NY, 1) = SP( 2, NY-1, 2)
SP( NX, NY, 1) = SP( NX-1,NY-1, 2)
SP( NX, NY, NZ)= SP( NX-1,NY-1, NZ-1)
SP( 1, NY, NZ)= SP( 2, NY-1, NZ-1)
C
999 ITR = ITR + 1
ITN = 0
C
C ----------------------------------------------------------
DO 400 K = 1, NZ
DO 400 J = 1, NY
DO 400 I = 1, NX
C ----------------------------------------------------------
C
IF ( I.EQ.1) THEN
EU( 1) = SU( I,J,K)
EP( 1) = SP( I,J,K)
ELSE
EU( 1) = SU( I-1,J,K)
EP( 1) = SP( I-1,J,K)
END IF
C
IF ( I.EQ.NX ) THEN
EU( 2) = SU( I,J,K)
EP( 2) = SP( I,J,K)
ELSE
EU( 2) = SU( I+1,J,K)
EP( 2) = SP( I+1,J,K)
END IF
C
IF ( J.EQ.1) THEN
EU( 3) = SV( I,J,K)
EP( 3) = SP( I,J,K)
ELSE
EU( 3) = SV( I,J-1,K)
EP( 3) = SP( I,J-1,K)
END IF
C
IF ( J.EQ.NY ) THEN
EU( 4) = SV( I,J,K)
EP( 4) = SP( I,J,K)
ELSE
EU( 4) = SV( I,J+1,K)
EP( 4) = SP( I,J+1,K)
END IF
C
IF ( K.EQ.1) THEN
EU( 5) = SW( I,J,K)
EP( 5) = SP( I,J,K)
ELSE
EU( 5) = SW( I,J,K-1)
EP( 5) = SP( I,J,K-1)
END IF
C
IF ( K.EQ.NZ ) THEN
EU( 6) = SW( I,J,K)
EP( 6) = SP( I,J,K)
ELSE
EU( 6) = SW( I,J,K+1)
EP( 6) = SP( I,J,K+1)
END IF
C
C -------------------------------
IF ( I.EQ.1) THEN
II = I+1
ELSE
II = I
END IF
C --------------------------------
IF ( J.EQ.1) THEN
JJ = J+1
ELSE
JJ = J
END IF
C --------------------------------
IF ( K.EQ.1) THEN
KK = K+1
ELSE
KK = K
END IF
C --------------------------------
C
E1 = ( 2.D0 / HX( II-1)/( HX( II-1) + HX( I)))
E2 = ( -2.D0 / HX( II-1)/ HX( I))
E3 = ( 2.D0 / HX( I) /( HX( II-1) + HX( I)))
C
E4 = ( 2.D0 / HY( JJ-1)/( HY( JJ-1) + HY( J)))
E5 = ( -2.D0 / HY( JJ-1)/ HY( J))
E6 = ( 2.D0 / HY( J) /( HY( JJ-1) + HY( J)))
C
E7 = ( 2.D0 / HZ( KK-1)/( HZ( KK-1) + HZ( K)))
E8 = ( -2.D0 / HZ( KK-1)/ HZ( K))
E9 = ( 2.D0 / HZ( K) /( HZ( KK-1) + HZ( K)))
C
G1 = ( - HX( I)/ HX( II-1)/( HX( II-1) + HX( I)))
G2 = (( HX( I)- HX( II-1))/ HX( II-1)/ HX( I))
G3 = ( HX( II-1)/ HX( I)/( HX( II-1) + HX( I)))
C
G4 = ( - HY( J)/ HY( JJ-1)/( HY( JJ-1) + HY( J)))
G5 = (( HY( J)- HY( JJ-1))/ HY( JJ-1)/ HY( J))
G6 = ( HY( JJ-1)/ HY( J)/( HY( JJ-1) + HY( J)))
C
G7 = ( - HZ( K)/ HZ( KK-1)/( HZ( KK-1) + HZ( K)))
G8 = (( HZ( K)- HZ( KK-1))/ HZ( KK-1)/ HZ( K))
G9 = ( HZ( KK-1)/ HZ( K)/( HZ( KK-1) + HZ( K)))
C
B( I,J,K)= ( G1* EU( 1) + G2* SU( I,J,K) + G3* EU( 2)
1 + G4* EU( 3) + G5* SV( I,J,K) + G6* EU( 4)
2 + G7* EU( 5) + G8* SW( I,J,K) + G9* EU( 6))/ DT
C
SSP = ( B( I,J,K) - E1* EP( 1) - E3* EP( 2) - E4* EP( 3)
1 - E6* EP( 4) - E7* EP( 5) - E9* EP( 6))
2 /( E2 + E5 + E8 )
C
ORM = DABS( SSP - SP( I,J,K))
C
SP( I,J,K) = SP( I,J,K) + AP* ( SSP - SP( I,J,K))
C ======================================================
C
IF ( ORM.GT.5.D-4) ITN = 1
C
400 CONTINUE
C ----------------------------------------------------------
C
IF ( ITN.EQ.0) GO TO 333
IF ( ITR.LT.( NX-1)* 100) GO TO 999
C
WRITE(6,2000) LOP, ITR, ORM
2000 FORMAT(' *** Diverged & Stop *** LOP=',I4,' * ITR =',I5,
1 '* ORM =',F10.5 )
STOP
C
333 IF ( KY.EQ.1) IT1 = ITR
IF ( KY.EQ.2) IT2 = ITR
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,*) '*****************************************************'
WRITE(6,*) '* *'
WRITE(6,*) '* F 3 D - F D M *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( V.4.0 ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* 3D Fluid Flow Analysis Software *'
WRITE(6,*) '* *'
WRITE(6,*) '* ( Fourth-ordered F.D.M. ) *'
WRITE(6,*) '* *'
WRITE(6,*) '* Copyright 2011 : Yasuhiro MATSUDA *'
WRITE(6,*) '* *'
WRITE(6,*) '*****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VARC ( U, V, W, P, OU, OV, OW, OP, NX,NY,NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
1 P ( NX,NY,NZ), OU( NX,NY,NZ), OV( NX,NY,NZ),
2 OW( NX,NY,NZ), OP( NX,NY,NZ)
C
COMMON /VVP/ VRV( 3000), VRP( 3000)
COMMON /ABC/ PMI, PMX, UMX, EPS, LOP, ITR, LOPE, NPR
C
EPS = 0.D0
C
C --------------------------------------
DO 100 K = 2, NZ-1
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C
IF ( DABS( U( I,J,K)).GT.0.01)
1 EPS = DMAX1( EPS, DABS(( OU( I,J,K) - U( I,J,K))/ U( I,J,K)))
C
IF ( DABS( V( I,J,K)).GT.0.01)
1 EPS = DMAX1( EPS, DABS(( OV( I,J,K) - V( I,J,K))/ V( I,J,K)))
C
IF ( DABS( W( I,J,K)).GT.0.01)
1 EPS = DMAX1( EPS, DABS(( OW( I,J,K) - W( I,J,K))/ W( I,J,K)))
C
100 CONTINUE
C --------------------------------------
C
EPSP = 0.D0
C --------------------------------------
DO 200 K = 2, NZ-1
DO 200 J = 2, NY-1
DO 200 I = 2, NX-1
C
IF ( DABS( P( I,J,K)).GT.0.01D0)
1 EPSP = DMAX1( EPSP, DABS(( OP( I,J,K)-P( I,J,K))/ P( I,J,K)))
C
200 CONTINUE
C --------------------------------------
C
C ----- (%) ---
C
VRV( LOP) = EPS * 100.D0
VRP( LOP) = EPSP* 100.D0
C --------------------------------------
C
PMI = 1000.D0
PMX = -1000.D0
UMX = -1000.D0
C
C --------------------------------------
DO 300 K = 1, NZ
DO 300 J = 1, NY
DO 300 I = 1, NX
C
PMI = DMIN1( PMI, P( I,J,K))
PMX = DMAX1( PMX, P( I,J,K))
UMX = DMAX1( UMX, U( I,J,K))
C
300 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************

C
SUBROUTINE VLSF ( U, V, W, UA, VA, WA, PA, HX, HY, HZ,
1 NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), UA( NX,NY,NZ), V ( NX,NY,NZ),
1 VA( NX,NY,NZ), W ( NX,NY,NZ), WA( NX,NY,NZ),
2 PA( NX,NY,NZ),
3 HX( NX), HY( NY), HZ( NZ)
C
COMMON /PAR/ DT, RE, WEI
C
C ----------------------------------------------------------
DO 100 K = 2, NZ-1
C
HZD = HZ( K-1) + HZ( K)
C
E7 = 2.D0 / HZ( K-1)/ HZD
E8 = - 2.D0 / HZ( K-1)/ HZ( K)
E9 = 2.D0 / HZ( K) / HZD
C
G7 = - HZ( K)/( HZ( K-1)* HZD )
G8 = ( HZ( K) - HZ( K-1))/( HZ( K-1)* HZ( K))
G9 = HZ( K-1)/( HZ( K)* HZD )
C
C ----------------------------------------------------------
DO 100 J = 2, NY-1
C
HYD = HY( J-1) + HY( J)
C
E4 = 2.D0/ HY( J-1)/ HYD
E5 = - 2.D0/ HY( J-1)/ HY( J)
E6 = 2.D0/ HY( J) / HYD
C
G4 = - HY( J)/( HY( J-1)* HYD )
G5 = ( HY( J) - HY( J-1))/( HY( J-1)* HY( J))
G6 = HY( J-1)/( HY( J)* HYD )
C
C ----------------------------------------------------------
DO 100 I = 2, NX-1
C
HXD = HX( I-1) + HX( I)
C
E1 = 2.D0/ HX( I-1)/ HXD
E2 = - 2.D0/ HX( I-1)/ HX( I)
E3 = 2.D0/ HX( I) / HXD
C
G1 = - HX( I)/( HX( I-1)* HXD )
G2 = ( HX( I) - HX( I-1))/( HX( I-1)* HX( I))
G3 = HX( I-1)/( HX( I)* HXD )
C
C ----- Correction coefficient of Diffusivity for 4th FDM ----
C
BX = U( I,J,K)* DT/ HX( I)
BY = V( I,J,K)* DT/ HY( J)
BZ = W( I,J,K)* DT/ HZ( K)
C
RX = DT/( RE* HX( I)**2 )
RY = DT/( RE* HY( J)**2 )
RZ = DT/( RE* HZ( K)**2 )
C
F = 1.D0 + ( BX* BY + BX* BZ + BY* BZ)/( RX + RY + RZ)
C ============================================================
C
C ----- Latter half of Expr.(8), (9) and (10) ---
C
UA( I,J,K) = UA( I,J,K)
C
1 +( F* DT/RE)*( E1* U( I-1,J,K) + E2* U( I,J,K) + E3* U( I+1,J,K)
2 + E4* U( I,J-1,K) + E5* U( I,J,K) + E6* U( I,J+1,K)
3 + E7* U( I,J,K-1) + E8* U( I,J,K) + E9* U( I,J,K+1))
C
4 - DT* ( G1* PA( I-1,J,K) + G2* PA( I,J,K) + G3* PA( I+1,J,K))
C
C --------------------------------------
C
VA( I,J,K) = VA( I,J,K)
C
1 + ( F*DT/RE)*( E1* V( I-1,J,K) + E2* V( I,J,K) + E3* V( I+1,J,K)
2 + E4* V( I,J-1,K) + E5* V( I,J,K) + E6* V( I,J+1,K)
3 + E7* V( I,J,K-1) + E8* V( I,J,K) + E9* V( I,J,K+1))
C
4 - DT* ( G4* PA(I,J-1,K) + G5* PA( I,J,K)+ G6* PA( I,J+1,K))
C
C --------------------------------------
C
WA( I,J,K) = WA( I,J,K)
C
1 + ( F*DT/RE)*( E1* W( I-1,J,K) + E2* W( I,J,K) + E3* W( I+1,J,K)
2 + E4* W( I,J-1,K) + E5* W( I,J,K) + E6* W( I,J+1,K)
3 + E7* W( I,J,K-1) + E8* W( I,J,K) + E9* W( I,J,K+1))
C
4 - DT* ( G7* PA( I,J,K-1) + G8* PA( I,J,K) + G9* PA( I,J,K+1))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSX ( U, V, W, UA, VA, WA, FU, FV, FW,
1 AX, HX, NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), V ( NX,NY,NZ), W ( NX,NY,NZ),
1 UA( NX,NY,NZ), VA( NX,NY,NZ), WA( NX,NY,NZ),
2 FU( NX+3), FV( NX+3), FW( NX+3),
3 AX( NX,NY,NZ), HX( NX )
C
COMMON /PAR/ DT, RE, WEI
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C --- du/dt: -u*( du/dx) dv/dt: -u*( dv/dx) dw/dt: -u*( dw/dx) ---
C
C ----------------------------------------------------------
DO 100 K = 2, NZ-1
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C ----------------------------------------------------------
C
F1 = ( - HX( I)/ HX( I-1) /( HX( I-1) + HX( I)))
F2 = (( HX( I) - HX( I-1))/ HX( I-1)/ HX( I))
F3 = ( HX( I-1)/ HX( I) /( HX( I-1) + HX( I)))
C
AX( I,J,K) = U( I,J,K)* DT
C ===============================
C
IF ( I.LE.4.OR.I.EQ.NX-1) GO TO 50
GO TO 60
C
50 CONTINUE
C
UA( I,J,K) = U( I,J,K) - AX( I,J,K)
1 * ( F1* U( I-1,J,K) + F2* U( I,J,K) + F3* U( I+1,J,K))
C
VA( I,J,K) = V( I,J,K) - AX( I,J,K)
1 * ( F1* V( I-1,J,K) + F2* V( I,J,K) + F3* V( I+1,J,K))
C
WA( I,J,K) = W( I,J,K) - AX( I,J,K)
1 * ( F1* W( I-1,J,K) + F2* W( I,J,K) + F3* W( I+1,J,K))
C
GO TO 100
C
60 CONTINUE
C
CALL COFS ( HX, NX, I )
C
CALL FCFX ( FU, FV, FW, U, V, W, AX, NX, NY, NZ, I, I, J, K )
C
CALL FCFX ( FU, FV, FW, U, V, W, AX, NX, NY, NZ, I+1, I, J, K )
C
C ----- Former half of Expr.(1.8), (1.9) and (1.10) (x-Dir.) ---
C
UA( I,J,K) = WEI* U( I,J,K) + ( 1.D0 - WEI )* U( I-1,J,K)
C ---------------------------------------------
1 + WEI* ( FU( I) - FU( I+1))
** 1 + ( 1.D0 - WEI )* ( FU( I+2) - FU( I+3))
C
VA( I,J,K) = WEI* V( I,J,K) + ( 1.D0 - WEI )* V( I,J-1,K)
C --------------------------------------------
1 + WEI* ( FV( I) - FV( I+1))
** 1 + ( 1.D0 - WEI )* ( FV( I+2) - FV( I+3))
C
WA( I,J,K) = WEI* W( I,J,K) + ( 1.D0 - WEI )* W( I,J,K-1)
C --------------------------------------------
1 + WEI* ( FW( I) - FW( I+1))
** 1 + ( 1.D0 - WEI )* ( FW( I+2) - FW( I+3))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSY ( U, V, W, UA, VA, WA, FU, FV, FW,
1 AY, HY, NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), UA( NX,NY,NZ), V ( NX,NY,NZ),
1 VA( NX,NY,NZ), W ( NX,NY,NZ), WA( NX,NY,NZ),
2 HY( NY), FU( NY + 3), FV( NY + 3),
3 FW( NY + 3), AY( NX,NY,NZ )
C
COMMON /PAR/ DT, RE, WEI
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C du/dt : -v* ( du/dy) dv/dt : -v* ( dv/dy) dw/dt : -v* ( dw/dy)
C
C ----------------------------------------------------------
DO 100 K = 2, NZ-1
DO 100 J = 2, NY-1
C ----------------------------------------------------------
C
F1 = ( - HY( J)/ HY( J-1) /( HY( J-1) + HY( J)))
F2 = (( HY( J)- HY( J-1))/ HY( J-1)/ HY( J))
F3 = ( HY( J-1)/ HY( J) /( HY( J-1) + HY( J)))
C
C ----------------------------------------------------------
DO 100 I = 2, NX-1
C
AY( I,J,K) = V( I,J,K)* DT
C ===============================
C
IF ( J.LE.4.OR.J.EQ.NY-1) GO TO 50
GO TO 60
C
50 CONTINUE
C
UA( I,J,K) = UA( I,J,K) - AY( I,J,K)
1 * ( F1* U( I, J-1, K) + F2* U( I,J,K) + F3* U( I, J+1, K))
C
VA( I,J,K) = VA( I,J,K) - AY( I,J,K)
1 * ( F1* V( I, J-1, K) + F2* V( I,J,K) + F3* V( I, J+1, K))
C
WA( I,J,K) = WA( I,J,K) - AY( I,J,K)
1 * ( F1* W( I, J-1, K) + F2* W( I,J,K) + F3* W( I, J+1, K))
GO TO 100
C
60 CONTINUE
C
CALL COFS ( HY, NY, J )
C
CALL FCFY ( FU, FV, FW, U, V, W, AY, NX, NY, NZ, J, I, J, K )
C
CALL FCFY ( FU, FV, FW, U, V, W, AY, NX, NY, NZ, J+1, I, J, K )
C
C ----- Former half of Expr.(1.8),(1.9),(1.10) (y-Dir.) ---
C
UA( I,J,K) = UA( I,J,K) + WEI* ( FU( J) - FU( J+1))
** 1 + ( 1.D0 - WEI )* ( FU( J+2) - FU( J+3))
C
VA( I,J,K) = VA( I,J,K) + WEI* ( FV( J) - FV( J+1))
** 1 + ( 1.D0 - WEI )* ( FV( J+2) - FV( J+3))
C
WA( I,J,K) = WA( I,J,K) + WEI* ( FW( J) - FW( J+1))
** 1 + ( 1.D0 - WEI )* ( FW( J+2) - FW( J+3))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VLSZ ( U, V, W, UA, VA, WA, FU, FV, FW, AZ,
1 HZ, NX, NY, NZ )
C
IMPLICIT REAL*8( A- H,O-Z )
C
DIMENSION U ( NX,NY,NZ), UA( NX,NY,NZ), V ( NX,NY,NZ),
1 VA( NX,NY,NZ), W ( NX,NY,NZ), WA( NX,NY,NZ),
2 HZ( NZ), FU( NZ+3), FV( NZ+3),
3 FW( NZ+3), AZ( NX,NY,NZ )
C
COMMON /PAR/ DT, RE, WEI
COMMON /CCF/ AG1, AG2, AG3, AG4, BG1, BG2, BG3, BG4,
1 CG1, CG2, CG3, CG4, DG1, DG2, DG3, DG4
C
C --- du/dt = -w*( du/dz) dv/dt = -w*( dv/dz) dw/dt = -w*( dw/dz) ---
C
C ----------------------------------------------------------
DO 100 K = 2, NZ-1
C
F1 = ( - HZ( K) / HZ( K-1) /( HZ( K-1) + HZ( K)))
F2 = ( ( HZ( K) - HZ( K-1))/ HZ( K-1)/ HZ( K))
F3 = ( HZ( K-1)/ HZ( K) /( HZ( K-1) + HZ( K)))
C
C ----------------------------------------------------------
DO 100 J = 2, NY-1
DO 100 I = 2, NX-1
C ----------------------------------------------------------
C
AZ( I,J,K) = W( I,J,K)* DT
C ===============================
C
IF ( K.LE.4.OR.K.EQ.NZ-1) GO TO 50
GO TO 60
C
50 CONTINUE
UA( I,J,K) = UA( I,J,K) - AZ( I,J,K)
1 * ( F1* U( I,J,K-1) + F2* U( I,J,K) + F3* U( I,J,K+1))
C
VA( I,J,K) = VA( I,J,K) - AZ( I,J,K)
1 * ( F1* V( I,J,K-1) + F2* V( I,J,K) + F3* V( I,J,K+1))
C
WA( I,J,K) = WA( I,J,K) - AZ( I,J,K)
1 * ( F1* W( I,J,K-1) + F2* W( I,J,K) + F3* W( I,J,K+1))
GO TO 100
C
60 CONTINUE
C
CALL COFS ( HZ, NZ, K )
C
CALL FCFZ ( FU, FV, FW, U, V, W, AZ, NX,NY,NZ, K, I, J, K )
C
CALL FCFZ ( FU, FV, FW, U, V, W, AZ, NX,NY,NZ, K+1, I, J, K )
C
C ----- Former half of Expr.(1.8), (1.9) and (1.10) (z-Dir.) ---
C
UA( I,J,K) = UA( I,J,K) + WEI* ( FU( K) - FU( K+1))
** 1 + ( 1.D0 - WEI )* ( FU( K+2) - FU( K+3))
C
VA( I,J,K) = VA( I,J,K) + WEI* ( FV( K) - FV( K+1))
** 1 + ( 1.D0 - WEI )* ( FV( K+2) - FV( K+3))
C
WA( I,J,K) = WA( I,J,K) + WEI* ( FW( K) - FW( K+1))
** 1 + ( 1.D0 - WEI )* ( FW( K+2) - FW( K+3))
C
100 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************