戻る

C ******************************************************************** C
C C
PROGRAM MAIN
C C
C *** F 2 D - M A C *** C
C C
C ( Two-Dimensional Viscous Fluid Analysis Software ) C
C C
C ( Finite Difference Method ) C
C C
C Continuous computation ( To change DT ） C
C C
C ( Data : 10 X 10 / 20 X 20 / 40 X 40 ) C
C C
C =====> C
C (0,0)+------+-----> X C
C | |(W,0) C
C | | C
C (0,H)+------+(W,H) C
C | C
C V C
C Y C
C C
C ( V.4.2 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /RTM/ RE, DT, T, LOP, N
C
C ----- DATA for 10 X 10, 20 X 20, 40 X 40 ( Max.: 50 X 50 ) ---
C
PARAMETER ( IX= 14, IY= 14) ! Re= 100 DT = 0.1
C
** PARAMETER ( IX= 24, IY= 24) ! Re= 400 DT = 0.05
C
** PARAMETER ( IX= 44, IY= 44) ! Re=1000 DT = 0.021
C
C --------------------------------------------------------------------
C
DIMENSION U ( IX,IY,2), V ( IX,IY,2), PR ( IX,IY),
1 D ( IX,IY), R ( IX,IY), UT ( IX,IY),
2 VT ( IX,IY), OVL( IX,IY), CON( IX,IY),
C
3 TST( IX,IY), CST( IX,IY),
4 XX (( IX-3)* ( IY-3)), YY(( IX-3)* ( IY-3)),
5 NDE((IX-4)* (IY-4), 4)
C
C ----- Attention : Max. of LOP < = 2000 ( Sub."CHCK" ) ---
C
DIMENSION RMX( 2000)
C
C ----- Attention : AAA( IX* IY* 2) in "VCST" ---
C
OPEN ( 3, FILE='R100_10_UVVAL.DATA') ! DT = 0.1
OPEN ( 7, FILE='R100_10_F2D_MAC.DATA')
OPEN ( 8, FILE='R100_10_VARP.DATA')
C
** OPEN ( 3, FILE='R100_40_UVVAL.DATA')
** OPEN ( 7, FILE='R100_40_F2D_MAC.DATA')
** OPEN ( 8, FILE='R100_40_VARP.DATA')
C
** OPEN ( 3, FILE='R400_20_UVVAL.DATA') ! DT = 0.05
** OPEN ( 7, FILE='R400_20_F2D_MAC.DATA')
** OPEN ( 8, FILE='R400_20_VARP.DATA')
C
** OPEN ( 3, FILE='R400_40_UVVAL.DATA') ! DT = 0.035
** OPEN ( 7, FILE='R400_40_F2D_MAC.DATA')
** OPEN ( 8, FILE='R400_40_VARP.DATA')
C
** OPEN ( 3, FILE='R1000_40_UVVAL.DATA') ! DT = 0.021
** OPEN ( 7, FILE='R1000_40_F2D_MAC.DATA')
** OPEN ( 8, FILE='R1000_40_VARP.DATA')
C
CALL PREP ( IX, IY)
C
C --------------------------------------------------------------------
999 CONTINUE
C
C --------------------------------------------------------------------
C
CALL INIT ( U, V, PR, D, R, UT, VT, XX, YY, NDE,
1 IX, IY, L, K )
C
C --------------------------------------------------------------------
T = 0.D0
NTL = 0
C
C ----------------------------------------------------------
DO 100 N = 1, LOP
C
T = T + DT
C
C ----- D & R from U & V -------------------------
C
CALL DRCL ( U, V, K, L, D, IX, IY, R )
C
C ----- Pressure ---------------------------------------------------
C
CALL PCAL ( U, V, L, R, PR, CON, TST, CST, IX, IY )
C
C ----- U & V ------------------------------------------------------
C
CALL UVCL ( U, V, PR, K, L, IX, IY )
C
C ----------------------------------------------------------------
C
CALL CHCK ( RMX, UT, VT, U, V, OVL, K, IX, IY, ICK,
1 NTL )
C
C ----------------------------------------------------------------
IF ( ICK.NE.0) GO TO 990
GO TO ( 10, 20), L
10 L = 2
K = 1
GO TO 100
20 L = 1
K = 2
C
100 CONTINUE
C ----------------------------------------------------------
C
990 CONTINUE
C -------------------------------------------------------------------
C
CALL PLTO ( RMX, UT, VT, PR, XX, YY, NDE, IX, IY,
1 ICK, ICP, NTL )
C
C -------------------------------------------------------------------
C
IF ( ICP.EQ.1) GO TO 999
C -------------------------------------------------------------------
CLOSE ( 3)
CLOSE ( 7)
CLOSE ( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CHCK ( RMX, UT, VT, U, V, OVL, K, IX, IY, ICK, NTL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /UVP/ PCN, MXT, ITO, ITI
COMMON /RMA/ RMD, ISS, JSS, RML
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
C
DIMENSION UT ( IX,IY), VT ( IX,IY), U( IX,IY,2), V( IX,IY,2),
1 OVL( IX,IY)
C
C ----- Attention : Max. of LOP < = 2000 ( Sub."CHCK" ) ---
C
DIMENSION RMX( 2000)
C
ISS = 0
JSS = 0
RMD = 0.D0
ICK = 0
C
C ----- CHCK of Divergence ---
C
IF ( ITO.GT.100) WRITE(6,*) ' ** ITO > = 100 ** STOP *** '
IF ( ITO.GT.100) ICK = 3
IF ( ITO.GT.100) GO TO 900
IF ( N.NE.1) GO TO 333
C
C ------------------------------------------------
DO 100 I = 3, II
DO 100 J = 3, JJ
C ------------------------------------------------
C
AV = ( V( I,J,K) + V( I-1,J,K))/ 2.D0
AU = ( U( I,J,K) + U( I,J-1,K))/ 2.D0
AUV = DSQRT( AV* AV + AU* AU)
OVL( I,J) = AUV
C
100 CONTINUE
C ------------------------------------------------
GO TO 900
C
333 CONTINUE
C
RML = 0.D0
C ------------------------------------------------
DO 200 I = 3, II
DO 200 J = 3, JJ
C ------------------------------------------------
C
UT( I,J) = ( U( I,J,K) + U( I,J-1,K))/ 2.D0
VT( I,J) = ( V( I-1,J,K) + V( I,J,K) )/ 2.D0
ABSV = DSQRT( UT( I,J)* UT( I,J) + VT( I,J)* VT( I,J))
IF ( ABSV.GT.RML) RML = ABSV
C
200 CONTINUE
C ------------------------------------------------
C
C ----- RML : Max. Velocity Value ---
C
RML = RML* 0.01D0
C
C ------------------------------------------------
DO 300 I = 3, II
DO 300 J = 3, JJ
C ------------------------------------------------
C
AV = ( V( I,J,K) + V( I-1,J,K))/ 2.D0
AU = ( U( I,J,K) + U( I,J-1,K))/ 2.D0
AUV = DSQRT( AV* AV + AU* AU)
C
C ----- > = Max.Vel(1 ) X 0.01 ---
C
IF ( AUV - RML) 195, 250, 250
C
250 CONTINUE
DFR = DABS(( AUV - OVL( I,J))/ OVL( I,J))
IF ( DFR.LE.RMD) GO TO 195
RMD = DFR
ISS = I
JSS = J
195 CONTINUE
OVL( I,J) = AUV
C
300 CONTINUE
C ------------------------------------------------
C
C ------ Max. Rel. Vel Var. Ratio = < 0.1 (%) ---
C
IF ( RMD.LE.0.001D0) GO TO 350
C ------------------------------------------------
DO 400 I = 3, II
DO 400 J = 3, JJ
C ------------------------------------------------
C
IF ( DABS( V( I,J,K)) - 1.1D0) 400, 50, 50
50 WRITE(6,2000) I, J
2000 FORMAT(' *** Diverged ! : V.GE. 1.1D0 ','I =',I3,'J =',I3)
C
C ----- Diverged ---
C
ICK = 2
GO TO 900
C
400 CONTINUE
C ------------------------------------------------
GO TO 900
C
C ----- Stationary condition reached ---
C
350 ICK = 1
900 CONTINUE
C
C ----- (%) ----------
C
RMX( N) = RMD* 100.D0
IF ( N.EQ.1) WRITE(6,2100) T
2100 FORMAT(' ** TIME =',F7.3)
NTL = NTL + ITO
C --------------------
C
NPR = 50
IF ( MOD( N,NPR).EQ.0)
1 WRITE(6,2200) N, T, ITO, RMX( N)
2200 FORMAT(' ** LOP =',I4,' TIME =',F7.3,' ITR=', I4,' Max.Var.of Ve
1l.=', F8.3,'(%)')
C
C ----- Stationary condition reached ---
C
IF ( ICK.EQ.1) THEN
WRITE(6,*) ' '
WRITE(6,*) '** Stationary condition reached **'
WRITE(6,2200) N, T, ITO, RMX( N)
END IF
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CSTB ( IX, IY)
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /GRV/ GY, GDY
COMMON /CTR/ OMG, CMG, OMZ
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
COMMON /VSK/ VIC, VSX, VSY, VX, VY, VXY
COMMON /XYT/ DEX, DEY, DXQ, DYQ, DX4, DY4, TXY, FDX, FDY, DYX, DXY
C
DEX = 1.D0/ DFLOAT( IX-4)
DEY = 1.D0/ DFLOAT( IY-4)
C
DXQ = DEX** 2
DYQ = DEY** 2
DX4 = 4.D0* DXQ
DY4 = 4.D0* DYQ
C
GY = 0.D0
GDY = 0.D0
C
VSX = 2.D0* VIC/ DEX
VSY = 2.D0* VIC/ DEY
VY = VIC/ DYQ
VX = VIC/ DXQ
VXY = VIC/( DEX* DEY)
C
C ----- OMG : Relaxation factor ---
C
CMG = 1.D0 - OMG
OMZ = OMG/( 2.D0* ( 1.D0/ DXQ + 1.D0/ DYQ))
C
FDX = 4.D0* DEX
FDY = 4.D0* DEY
C
TXY = 2.D0* DEX* DEY
DXY = DEX/ DEY
DYX = DEY/ DEX
C
IP1 = II + 1
IM1 = II - 1
C
JP1 = JJ + 1
JM1 = JJ - 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE DRCL ( U, V, K, L, D, IX, IY, R )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /GRV/ GY, GDY
COMMON /CTR/ OMG, CMG, OMZ
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
COMMON /VSK/ VIC, VSX, VSY, VX, VY, VXY
COMMON /XYT/ DEX, DEY, DXQ, DYQ, DX4, DY4, TXY, FDX, FDY, DYX, DXY
C
DIMENSION U( IX,IY,2), V( IX,IY,2), D( IX,IY), R( IX,IY)
C
DMAX = 0.D0
C
C ------------------------------------------------
DO 100 J = 3, JJ
DO 100 I = 3, II
C ------------------------------------------------
C
D( I,J) = (( U( I,J,L) - U( I,J-1,L))/ DEX )
1 + (( V( I,J,L) - V( I-1,J,L))/ DEY )
C
IF ( DABS( D( I,J)).LE.DMAX) GO TO 50
DMAX = DABS( D( I,J))
SDMX = D( I,J)
C
C ----- General term ---
C
50 CONTINUE
RV1 = (( U( I,J+1,L) + U( I,J,L))** 2
1 - 2.D0 * ( U( I,J,L) + U( I,J-1,L))** 2
2 + ( U( I,J-1,L) + U( I,J-2,L))** 2 )/ DX4
C
3 + (( V( I+1,J,L) + V( I, J,L))** 2
4 - 2.D0 * ( V( I, J,L) + V( I-1,J,L))** 2
5 + ( V( I-1,J,L) + V( I-2,J,L))** 2 )/ DY4
C
RV2 = ((( U( I+1,J, L) + U( I,J,L))* ( V( I,J+1,L) + V( I,J,L))
1 + ( U( I,J-1,L) + U( I-1,J-1,L))
2 * ( V( I-1,J,L) + V( I-1,J-1,L))
C
3 - ( U( I,J,L) + U( I-1,J,L))* ( V( I-1,J+1,L) + V( I-1,J,L))
4 - ( U( I+1,J-1,L) + U( I,J-1,L))
5 * ( V( I,J,L) + V( I,J-1,L)))/ TXY ) - D( I,J)/ DT
C
R( I,J) = RV1 + RV2
C
100 CONTINUE
C ------------------------------------------------
DO 200 J = 3, JJ
C
C ----- Above Upper Wall ---
C
R( 2,J)= 0.25D0* ( DYX
1 * (( U( 3,J,L) + U( 2,J,L))* ( V( 2,J+1,L) + V( 2,J,L))
2 - ( U( 3,J-1,L) + U( 2,J-1,L))* ( V( 2,J,L) + V( 2,J-1,L)))
C
3 + (( V( 3,J,L) + V( 2,J,L))** 2
4 - ( V( 2,J,L) + V( 1,J,L))** 2 )) - 0.5D0* VSX
5 * ( DYX* ( V( 2,J+1,L) - 2.D0* V( 2,J,L) + V( 2,J-1,L))
C
6 - ( U( 3,J,L) - U( 3,J-1,L) - U( 2,J,L) + U( 2,J-1,L)))
7 - GDY
C *************
C
C ----- Below Lower Wall ---
C
R( II+1,J)= - 0.25D0* (( ( V( II+1,J,L) + V( II, J,L))** 2
1 - ( V( II, J,L) + V( II-1,J,L))** 2 )
2 + DYX* (( U( II+1,J,L) + U( II,J,L))
3 * ( V( II,J+1,L) + V( II,J,L))
C
4 - ( U( II+1,J-1,L) + U( II,J-1,L))
5 * ( V( II, J, L) + V( II,J-1,L))))
C
6 + 0.5D0* VSX*
7 ( DYX* ( V( II,J+1,L) - 2.D0* V( II,J,L) + V( II,J-1,L))
8 - ( U( II+1,J,L) - U( II+1,J-1,L)
9 - U( II, J,L) + U( II, J-1,L))) + GDY
C
200 CONTINUE
C ------------------------------------------------
DO 300 I = 3, II
C
C ----- Before Left Wall ---
C
R( I,2)= 0.25D0* (( ( U( I,3,L) + U( I,2,L))** 2
1 - ( U( I,2,L) + U( I,1,L))** 2 )
C
2 + DXY* (( U( I+1,2,L) + U( I,2,L))
3 * ( V( I,3,L) + V( I,2,L)) - ( U( I-1,2,L) + U( I,2,L))
4 * ( V( I-1,3,L) + V( I-1,2,L))))
C
5 - 0.5D0* VSY
6 * ( DXY* ( U( I+1,2,L) - 2.D0* U( I,2,L) + U( I-1,2,L))
7 - ( V( I,3,L) - V( I-1,3,L) - V( I,2,L) + V( I-1,2,L)))
C
C ----- After Right Wall ---
C
R( I,JJ+1 )= - 0.25D0* (( ( U( I,JJ+1,L) + U( I,JJ,L))** 2
1 - ( U( I,JJ,L) + U( I,JJ-1,L))** 2 )
2 + DXY * (( U( I+1,JJ,L) + U( I,JJ,L))
3 * ( V( I,JJ+1,L) + V( I,JJ,L))
C
4 - ( U( I-1,JJ,L)+U( I,JJ,L))*( V( I-1,JJ+1,L)+V( I-1,JJ,L))))
C
5 + 0.5D0* VSY* ( DXY* ( U( I+1,JJ,L) - 2.D0* U( I,JJ,L)
6 + U( I-1,JJ,L)) - ( V( I,JJ+1,L) - V( I-1,JJ+1,L)
7 - V( I,JJ,L) + V( I-1,JJ,L)))
C
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE INIT ( U, V, PR, D, R, UT, VT, X, Y, NDE,
1 IX, IY, L, K )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY,2), V ( IX,IY,2), PR( IX,IY), D( IX,IY),
1 R( IX,IY), UT( IX,IY), VT( IX,IY)
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
COMMON /NNL/ NNP, NEL
COMMON /CTR/ OMG, CMG, OMZ
COMMON /UVP/ PCN, MXT, ITO, ITI
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
C
C ----- Setting including Boundfary Values ---
C
CALL STUV ( U, V, PR, D, R, UT, VT, IX, IY )
C
MXT = 100
ITO = 200
ITI = 5
C
PCN = 0.0001D0
OMG = 1.8D0
C
WRITE(6,*) ' '
WRITE(6,2000) IX, IY
2000 FORMAT(' *** IX=', I3,' * IY=', I3 )
WRITE(6,*) ' '
WRITE(6,2100)
2100 FORMAT(' * DT = ?')
C
** READ(5,*) DT
DT = 0.1D0
C
WRITE(6,2200) RE, LOP, DT, OMG
2200 FORMAT(' * Re =',F9.3,' LOP =',I4,' DT =',F7.3,
1 ' OMG =',F6.2,/)
C
II = IX - 2
JJ = IY - 2
C ----------------------------------------
C
CALL CSTB ( IX, IY)
C
C ----------------------------------------
C
CALL STNDE ( IX, IY, X, Y, NDE )
C
C ----------------------------------------
L = 1
K = 2
C
RETURN
END
C **********************************************************************
C
SUBROUTINE ITER ( PR, R, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION PR( IX,IY), R( IX,IY)
C
COMMON /CTR/ OMG, CMG, OMZ
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
COMMON /XYT/ DEX, DEY, DXQ, DYQ, DX4, DY4, TXY, FDX, FDY, DYX, DXY
C
C ----- Row above Upper Wall ---
C
C ------------------------------------------------
DO 100 J = 3, JJ
C
PR( 2,J) = PR( 3,J) + R( 2,J)
100 CONTINUE
C ------------------------------------------------
C
C ----- Column before Left Wall ---
C
C ------------------------------------------------
DO 200 I = 3, II
C
PR( I,2) = PR( I,3 ) + R( I,2)
C
C --------------------------------------
DO 250 J = 3, JJ
C
C ----- General Term ---
C
PR( I,J) = CMG* PR( I,J)
1 + OMZ* (( PR( I,J+1) + PR( I,J-1))/ DXQ
2 + ( PR( I+1,J) + PR( I-1,J))/ DYQ + R( I,J))
C
250 CONTINUE
C --------------------------------------
C
C ----- Column after Right Wall ---
C
PR( I,JJ+1) = PR( I,JJ) + R( I,JJ+1)
C
200 CONTINUE
C ------------------------------------------------
C
C ----- Row below Lower Wall ---
C
C ------------------------------------------------
DO 300 J = 3, JJ
C
PR( II+1,J) = PR( II,J) + R( II+1,J)
300 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCAL ( U, V, L, R, PR, CON, TST, CST, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Poisson's Equation Solver ---
C
DIMENSION R( IX,IY), U( IX,IY,2), V( IX,IY,2), PR( IX,IY)
C
COMMON /GRV/ GY, GDY
COMMON /UVP/ PCN, MXT, ITO, ITI
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
C
DIMENSION CON( IX,IY), TST( IX,IY), CST( IX,IY)
C
C ----- PCN = 0.0001D0 MXT = 100 ITO = 200 ITI = 5 ---
C
IL = 2 + ( II-2)/ 2
C
C ------------------------------------------------------------
DO 100 J = 2, JP1
DO 100 I = 2, IP1
C
CST( I,J) = 0.25D0* (( U( I,J,L) + U( I,J-1,L))** 2
1 + ( V( I,J,L) + V( I-1,J,L))** 2 )
C
100 CONTINUE
C ------------------------------------------------------------
C
IF ( ITO.NE.1) GO TO 150
IF ( ITI.EQ.5) GO TO 150
ITI = ITI - 1
150 CONTINUE
C
C ----------------------------------------------------------
DO 200 IOUT = 1, MXT
C
ITO = IOUT
C ------------------------------------------------
DO 220 IIN = 1, ITI
C
C ----------------------------------
C
CALL ITER ( PR, R, IX, IY )
C
C ----------------------------------
C
220 CONTINUE
C ------------------------------------------------
C
C ------------------------------------------------
DO 240 J = 2, JP1
DO 240 I = 2, IP1
C
TST( I,J) = PR( I,J)
C
240 CONTINUE
C ------------------------------------------------
C
CALL ITER ( PR, R, IX, IY)
C
C ----- Test for Convergence ---------------
C
CMAX = 0.D0
C ------------------------------------------------
DO 260 J = 3, JJ
DO 260 I = 3, II
C ------------------------------------------------
C
C ----- Convergence Check ---
C
CON( I,J) = DABS( PR ( I,J) - TST( I,J))/( DABS( PR( I,J))
1 + DABS( TST( I,J)) + CST( I,J))
IF ( CON( I,J).GT.CMAX) CMAX = CON( I,J)
C
260 CONTINUE
C ------------------------------------------------
C
IF ( CMAX.LE.PCN) GO TO 250 ! PCN = 0.0001D0
C
200 CONTINUE
C ----------------------------------------------------------
ITO = MXT
RETURN
C
250 CONTINUE
CHG = 1.D0 - 0.25D0* ( PR( IL, JJ) + PR( IL, JP1 )
1 + PR( IL+1, JJ) + PR( IL+1, JP1 ))
C
C ----- Normalization ---
C
C ----------------------------------------------------------
DO 300 J = 2, JP1
DO 300 I = 2, IP1
C
PR( I,J) = PR( I,J) + CHG
300 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTO ( RMX, UT, VT, PR, XX, YY, NDE, IX, IY,
1 ICK, ICP, NTL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NNL/ NNP, NEL
COMMON /RTM/ RE, DT, T, LOP, N
C
C ----- Attention : Max. of Loop < = 2000 ( Sub."CHCK" ) ---
C
DIMENSION RMX( 2000)
C
DIMENSION UT( IX,IY), VT ( IX,IY), PR( IX,IY), XX( NNP),
1 YY( NNP), NDE( NEL,4)
C
C --------------------------------------------
C
IF ( ICK.EQ.1.OR.ICK.EQ.2) GO TO 990
WRITE(7,2000) IX, IY, RE, DT, N-1
2000 FORMAT(2I3,F7.1,F7.4,I5)
WRITE(7,2100) ( RMX( I), I= 2, N-1)
C
WRITE(8,2200) N-1, DT
WRITE(8,2300) ( RMX( I), I= 2, N-1)
GO TO 998
C
990 CONTINUE
WRITE(7,2000) IX, IY, RE, DT, N
WRITE(7,2100) ( RMX( I), I= 2, N)
2100 FORMAT( 10D11.4)
WRITE(8,2200) N, DT
2200 FORMAT(I5,D8.3)
WRITE(8,2300) ( RMX( I), I= 2, N)
2300 FORMAT(D11.4)
C
998 CONTINUE
WRITE(7,2400) (( UT( I,J), J=1, IY), I=1, IX)
WRITE(7,2400) (( VT( I,J), J=1, IY), I=1, IX)
WRITE(7,2400) (( PR( I,J), J=1, IY), I=1, IX)
WRITE(7,2400) ( XX( I), I=1, NNP)
WRITE(7,2400) ( YY( I), I=1, NNP)
2400 FORMAT(10D11.4)
C
WRITE(7,2500) (( NDE( I,J), J = 1, 4), I = 1, NEL)
2500 FORMAT(30I4)
C
C ----- Comparison with other Data ---
C
IF ( ICK.EQ.0) LPN = N - 1
IF ( ICK.NE.0) LPN = N
C
WRITE(3,2600) IX, IY, RE, DT, LPN
2600 FORMAT(2I3,F7.1,F7.4,I4)
WRITE(3,2400) (( UT( I,J), J = 1, IY), I = 1, IX)
WRITE(3,2400) (( VT( I,J), J = 1, IY), I = 1, IX)
C
UMX = 0.D0
C
C ----------------------------------------------------------
DO 100 I = 1, IX
DO 100 J = 1, IY
C
IF ( DABS( UT(I,J)).GE.UMX ) UMX = DABS( UT(I,J))
C
100 CONTINUE
C ----------------------------------------------------------
C
WRITE(6,*) ' '
WRITE(6,2700) NTL, UMX
2700 FORMAT(' * Tot. Iter. No=', I6,' ** Umax.=',F7.3)
WRITE(6,2800)
2800 FORMAT(/' *** Computation to be continued ? ( 1: Yes 2: No)')
C
** READ (5,*) ICP
ICP = 2
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PREP ( IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON /NNL/ NNP, NEL
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /VSK/ VIC, VSX, VSY, VX, VY, VXY
C
WRITE(6,100)
WRITE(6,101)' ************************************************** '
WRITE(6,*) ' '
WRITE(6,101)' ##### ### ### ## ## ## #### '
WRITE(6,101)' ## # # ## ## ### ### #### ## ## '
WRITE(6,101)' #### # # # ### ## # ## ## ## ## '
WRITE(6,101)' ## # # ## ## ## ###### ## ## '
WRITE(6,101)' ## #### ### ## ## ## ## #### '
WRITE(6,101)' '
WRITE(6,101)' '
WRITE(6,101)' Two Dimensional Viscous Fluid Flow Analysis '
WRITE(6,101)' '
WRITE(6,101)' Software by Finite Difference Method '
WRITE(6,101)' '
WRITE(6,101)' ( V.4.2 ) '
WRITE(6,101)' '
WRITE(6,101)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,101)' '
WRITE(6,101)' ************************************************** '
100 FORMAT(/)
101 FORMAT(A55)
C
NNP = ( IX-3)* ( IY-3)
NEL = ( IX-4)* ( IY-4)
C
WRITE(6,*) ' '
WRITE(6,2000)
2000 FORMAT(' ',' Re = ?')
** READ(5,*) RE
RE = 100.D0
C
VIC = 1.D0/ RE
C -------------------
C
WRITE(6,2100)
2100 FORMAT(' ',' How many Loops ? ')
** READ(5,*) LOP
LOP = 200
C
WRITE(6,*) ' '
WRITE(6,2200) RE, LOP
2200 FORMAT(' * Re =',F9.3,' Loop =',I4)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STNDE ( IXX, IYY, X, Y, NDE )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), NDE( NEL,4)
C
COMMON /NNL/ NNP, NEL
C
IF ( IXX.EQ.IYY) GO TO 1000
STOP
C
1000 N = IXX - 4
DX = 1.D0/ DFLOAT( N)
DY = 1.D0/ DFLOAT( N)
C
NNP = 0
C ------------------------------------------------
DO 100 I = 1, N+1
DO 100 J = 1, N+1
C ------------------------------------------------
C
NNP = NNP + 1
X( NNP) = DX* DFLOAT( N+1-J)
Y( NNP) = DY* DFLOAT( I-1)
C
100 CONTINUE
C ------------------------------------------------
DO 200 J = 1, N
DO 200 I = 1, N
C ------------------------------------------------
C
NUM = ( J-1)* N + I
NDE( NUM,1) = ( J-1)* ( N+1) + I + 1
NDE( NUM,2) = ( J-1)* ( N+1) + I
NDE( NUM,3) = ( J-1)* ( N+1) + I + ( N+1)
NDE( NUM,4) = ( J-1)* ( N+1) + I + ( N+2)
C
200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE STUV ( U, V, PR, D, R, UT, VT, IX, IY )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY,2), V ( IX,IY,2), PR( IX,IY), D( IX,IY),
1 R( IX,IY), UT( IX,IY), VT( IX,IY)
C
INN = IX* IY
IN2 = INN* 2
C
C ----- Attention ---
C
CALL VCST ( U, 0.D0, 1, IN2 )
C
CALL VCST ( V, 0.D0, 1, IX+1 )
C
CALL VCST ( V, 1.D0, IX+2, IX*2-2 )
C
CALL VCST ( V, 0.D0, IX*2-1, IX*(IX+1)+1 )
C
CALL VCST ( V, 1.D0, IX*(IX+1)+2, IX*(IX+1)+2+(IX-4))
C
CALL VCST ( V, 0.D0, IX*IY*2-(IX-1)*(IY-1), IN2 )
C
CALL VCST ( PR, 1.D0, 1, INN )
CALL VCST ( D, 0.D0, 1, INN )
CALL VCST ( R, 0.D0, 1, INN )
C
CALL VCST ( UT, 0.D0, 1, INN )
CALL VCST ( VT, 0.D0, 1, INN )
C
C --- 10 x 10 -------------------------------------------------------
C
* DATA U / 392*0./, V/ 15*0., 11*1., 185*0., 11*1., 170*0./
* DATA PR /196*1./, D /196*0./, R/ 196*0./
C
C --- 20 x 20 -------------------------------------------------------
C
* DATA U / 1152*0./, V/ 25*0., 21*1., 555*0., 21*1., 530*0./
* DATA PR /576*1./, GX/576*0./, D /576*0./, R/ 576*0./
C
C --- 40 X 40 -------------------------------------------------------
C
* CALL VCST ( U, 0.D0, 1, 3872)
C
* CALL VCST ( V, 0.D0, 1, 45)
* CALL VCST ( V, 1.D0, 46, 86)
* CALL VCST ( V, 0.D0, 87, 1981)
* CALL VCST ( V, 1.D0, 1982, 2022)
* CALL VCST ( V, 0.D0, 2023, 3872)
C
* CALL VCST ( PR, 1.D0, 1, 1936)
* CALL VCST ( GX, 0.D0, 1, 1936)
* CALL VCST ( D, 0.D0, 1, 1936)
* CALL VCST ( R, 0.D0, 1, 1936)
C
* CALL VCST ( OVEL, 0.D0, 1, 1936)
C
* CALL VCST ( UT, 0.D0, 1, 1936 )
* CALL VCST ( VT, 0.D0, 1, 1936 )
C
* DATA U / 3872*0./, V / 45*0., 41*1., 1895*0., 41*1., 1850*0./
* DATA PR/ 1936*1./, GX / 1936*0./, D / 1936*0./, R / 1936*0./
C ---------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE UVCL ( U, V, PR, K, L, IX, IY )
C
C ----- U - V Velocities ---
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION U( IX,IY,2), V( IX,IY,2), PR( IX,IY)
C
COMMON /GRV/ GY, GDY
COMMON /RTM/ RE, DT, T, LOP, N
COMMON /VSK/ VIC, VSX, VSY, VX, VY, VXY
COMMON /SIZ/ II, JJ, IP1, JP1, IM1, JM1
COMMON /XYT/ DEX, DEY, DXQ, DYQ, DX4, DY4, TXY, FDX, FDY, DYX, DXY
C
GXA = 0.D0
C
C ----- U, V - Velocity ---
C
C ----------------------------------------------------------
DO 100 I = 3, II
DO 100 J = 2, JJ
C ----------------------------------------------------------
C
C ----- General Terms ---
C
U( I,J,K) = U( I,J,L)
C
1 + DT* (( ( U( I,J,L) + U( I,J-1,L))** 2
2 - ( U( I,J+1,L) + U( I,J,L ))** 2)/ FDX
C
3 + (( U( I,J,L) + U( I-1,J,L))* ( V( I-1,J+1,L) + V( I-1,J,L))
4 - ( U( I+1,J,L) + U( I,J,L))*( V( I,J+1,L) + V( I,J,L)))/ FDY
C
5 + GXA + ( PR( I,J) - PR( I,J+1 ))/ DEX
6 + VY * ( U( I+1,J,L) - 2.D0* U( I,J,L) + U( I-1,J,L))
7 - VXY* ( V( I,J+1,L) - V( I-1,J+1,L)- V( I,J,L)+ V( I-1,J,L)))
C
V( I,J,K) = V( I,J,L)
C
1 + DT* ((( U( I+1,J-1,L) + U( I,J-1,L))*( V( I,J,L) + V( I,J-1,L))
2 - ( U( I+1,J,L) + U( I,J,L))* ( V( I,J+1,L) + V( I,J,L)))/ FDX
C
5 + (( V( I,J,L) + V( I-1,J,L))** 2
6 - ( V( I+1,J,L) + V( I, J,L))** 2 )/ FDY
7 + GY + ( PR( I,J) - PR( I+1,J))/ DEY
C *******
C
8 + VX * ( V( I,J+1,L) - 2.D0* V( I,J,L) + V( I,J-1,L))
9 - VXY* ( U( I+1,J,L) - U( I,J,L)- U( I+1,J-1,L)+ U( I,J-1,L)))
C
100 CONTINUE
C ----------------------------------------------------------
C
C ----- Lower No-Slip ---
C
C --------------------------------------
DO 120 J = 3, JJ
C
V( II,J,K) = 0.D0
120 CONTINUE
C --------------------------------------
C
C ----- Right No-Slip Wall -------------
C
DO 140 I = 3, II
C
U( I,JJ,K) = 0.D0
140 CONTINUE
C --------------------------------------
C
C ----- Boundary values ---
C
C --------------------------------------
DO 160 J = 3, JM1
C
C ----- Above Upper, No-Slip Wall ---
C
U( 2,J,K) = - U( 3,J,K)
V( 1,J,K) = - V( 3,J,K)
C
C ----- Below, Lower, No-Slip Wall ---
C
U( II+1,J,K) = - U( II,J,K)
V( II+1,J,K) = - V( II-1,J,K)
C
160 CONTINUE
C --------------------------------------
C
V( 1,JJ,K) = - V( 3,JJ,K)
V( II+1,JJ,K) = - V( II-1,JJ,K)
C
C --------------------------------------
DO 200 I = 3, IM1
C
C ----- Left, Moving, No-Slip Wall ---
C
U( I,1,K)= - U( I,3,K)
C
C ----- Attention : Velocity direction : Right ---
C ****************************
C
V( I,2,K) = - 2.D0 - V( I,3,K)
C
C ----- Right - N0-Slip Wall ---
C
U( I,JJ+1,K) = - U( I,JJ-1,K)
V( I,JJ+1,K) = - V( I,JJ,K)
C
200 CONTINUE
C --------------------------------------
C
U( II,1,K) = - U( II,3,K)
U( II,JJ+1,K) = - U( II,JJ-1,K)
C
C ----- Left - Up ---
C
U( 2,2,K) = U( 3,2,K)
V( 2,2,K) = V( 2,3,K)
C
C ----- Left - Low ---
C
U( II+1,2,K) = U( II,2,K)
V( II,2,K) = V( II,3,K)
C
C ----- Right - UP ---
C
U( 2,JJ,K) = - U( 3,JJ,K)
V( 2,JJ+1,K) = V( 2,JJ,K)
C
C ----- Right - Low ---
C
U( II+1,JJ,K) = - U( II,JJ,K)
V( II,JJ+1,K) = V( II,JJ,K)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VCST ( AAA, V, IS, IE )
C
IMPLICIT REAL*8(A-H,O-Z)
C
C ----- Attention ( Max.Case ) ---
C
C DIMENSION AAA( IX* IY* 2):
C
C Max.: For 40 x 40 ==> 44* 44* 2 = 3872 ---
C
C Max.: For 50 x 50 ==> 54* 54* 2 = 5832 ---
C
DIMENSION AAA( 5832)
C
C ------------------------
DO 100 I = IS, IE
C
AAA( I) = V
C
100 CONTINUE
C ------------------------
C
RETURN
END
C *********************************************************************