戻る

C
PROGRAM MAIN
C
C ******************************************************************** C
C C
C Graphics for F2D-D42E C
C C
C ( Cavity Flow Analysis : Equal mesh sizes ) C
C C
C ( V.4.0 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
C ----- ( 1 ) -----------------------------------------------
C
PARAMETER ( IX= 11, IY= 11, IX1= IX-1, IY1= IY-1 )
C
*** PARAMETER ( IX= 21, IY= 21, IX1= IX-1, IY1= IY-1 )
*** PARAMETER ( IX= 41, IY= 41, IX1= IX-1, IY1= IY-1 )
C
C -----------------------------------------------------------
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
C
DIMENSION X ( IX1*IY1), Y ( IX1*IY1), UE( IX1*IY1),
1 VE( IX1*IY1), XY( 2,IX*IY), PE( IX*IY)
C
DIMENSION U( IX,IY), V( IX,IY)
C
COMMON /DEF/ LOP, DT, RE, FSM, WEI
C
INTEGER*2 ICR
C
C ----- ( 2 ) ------------------------------------------
C
OPEN ( 7, FILE='F2DD42_10_FIG.DAT' ) ! 10 X 10
C
*** OPEN ( 7, FILE='F2DD42_20_FIG.DAT' ) ! 20 X 20
*** OPEN ( 7, FILE='F2DD42_40_FIG.DAT' ) ! 40 X 40
C
OPEN ( 8, FILE='GH100.DAT' ) ! Re = 100
C
*** OPEN ( 8, FILE='GH400.DAT' ) ! Re = 400
*** OPEN ( 8, FILE='GH1000.DAT' ) ! Re = 1000
C
C ----------------------------
C
CALL TTLE
C
C ----------------------------
READ(7,1000) FSM, WEI, RE, DT, LOP, NEL, NNP, NBC
1000 FORMAT(A8,3F9.3,3I4,I2)
C ---------------------------------------------------------------------
C
CALL PRSB ( UE, VE, U, V, X, Y, XY, PE, IX, IY,
1 NNP, NEL, ICR, KKK, NBC )
C
C ---------------------------------------------------------------------
C
NPL = 0
C
CALL VGA@
C
IF ( KKK.EQ.1) GO TO 100
IF ( KKK.EQ.2) GO TO 200
C
C ----- At the Last LOP ---
C
100 CONTINUE
C ------------------------------------------------
C
CALL MSHFG
C
C ------------------------------------------------
C
CALL PLTN ( X, Y, UE, VE, NEL )
C
C ------------------------------------------------
C
ICR = 15
C ------------------------------------------------
C
CALL VRPT ( WMX, ICR )
C
C ------------------------------------------------
C
ICR = 12
C ------------------------------------------------
C
CALL VRPT ( DMX, ICR )
C
C ------------------------------------------------
C
CALL CCNT ( XY, PE, NDE, NNP, NEL )
C
C ------------------------------------------------
GO TO 900
200 CONTINUE
C ------------------------------------------------
C
CALL CMPT ( U, V, X, Y, IX, IY, NNP, NPL )
C
C ------------------------------------------------
C
900 CONTINUE
CLOSE ( 7)
CLOSE ( 8)
C
STOP
END
C **********************************************************************
C
SUBROUTINE CMPT ( UT, VT, X, Y, IM, JM, NNP, NPL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
DIMENSION UT( IM,JM), VT( IM,JM), X( NNP), Y( NNP)
C
COMMON /DEF/ LOP, DT, RE, FSM, WEI
C
IX1 = 0
IY1 = 0
IX2 = 639
IY2 = 479
ICR = 7
C
CALL CLEAR_SCREEN_AREA@( IX1, IY1, IX2, IY2, ICR )
C
KEY = 0
C
WRITE(6,2000) RE, DT, LOP
2000 FORMAT(/,' ** Re =',F7.1,' * DELT =',F7.4,' * Last LOP =',I4 )
C
C ----------------------------------------------------------
C
CALL RFCMP ( X, Y, UT, VT, IM, JM, NNP, KEY, NPL )
C
C ----------------------------------------------------------
C
KEY = KEY + 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE RFCMP ( X, Y, U, V, IM, JM, NNP, KEY, NPL )
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION X ( NNP), Y ( NNP), U ( IM,JM), V ( IM,JM)
DIMENSION GU( 17), GV( 17), GUY( 17), GVX( 17)
C
INTEGER*2 IX1, IY1, IX2, IY2, ICR, IICR
C
CHARACTER*80 STR0
C
DATA STR0 /'** Results of F2D-D42 ** ( Velocity Distribution comp
1ared with Ghia''s Data )'/
C
II = IM
IM1 = II - 1
C
IF ( KEY.EQ.0) ICR = 11
IF ( KEY.EQ.1) ICR = 13
IF ( KEY.EQ.2) ICR = 14
C
IICR = 16
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 10, 5, IICR )
C
C ----- Ghia's Data --- ( Re = 100 or 400 or 1000)---
C
C ----------------------------------------------------------
DO 100 I = 1, 17
C
READ(8,1000) GU( I), GUY( I), GV( I), GVX( I)
1000 FORMAT(4F10.5)
C
100 CONTINUE
C ----------------------------------------------------------
REWIND ( 8)
C
C ----- Attention : For Velocity direction change ---
C
DO 200 I = 1, 17
C
GV ( I) = - GV( I)
GVX( I) = 1.D0 - GVX( I)
C
200 CONTINUE
C ----------------------------------------------------------
C
X1 = 345
Y1 = 50
X0 = 20
Y0 = 50
C
SCL = 260.D0
X2 = SCL* 0.5D0
Y2 = SCL* 0.5D0
C
C ----- U by Ghia ---
C
IX1 = X0
IY1 = Y0
IX2 = X0 + SCL
IY2 = Y0 + SCL
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X0 + X2
IY1 = Y0 + SCL
IX2 = X0 + X2
IY2 = Y0
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
C ----------------------------------------------------
DO 300 I = 1, 16
C
IC = I + 1
IX1 = X0 + X2 + GU ( I)* X2
IY1 = Y0 + SCL - GUY( I)* SCL
IX2 = X0 + X2 + GU (IC)* X2
IY2 = Y0 + SCL - GUY(IC)* SCL
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, IICR )
C
300 CONTINUE
C ----------------------------------------------------
C
IX1 = 50
IY1 = 80 + NPL* 20
C
IX2 = 100
IY2 = IY1
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
NPL = NPL + 1
C ----------------------------------------------------
C
C ----- Computed U ---
C
C ------------------------------------------------
DO 400 J = 2, II-1
C
N = IM/ 2
UP = ( U( N,J) + U( N+1,J))/2.D0
I = J - 2
N1 = IM1* I + II/2
N2 = IM1* ( I+1) + II / 2
YP = ( Y( N1) + Y( N2)) / 2.D0
IX1 = X0 + X2 + UP* X2
IY1 = Y0 + SCL - YP* SCL
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
400 CONTINUE
C ------------------------------------------------
C
C ----- V by Ghia ---
C
IX1 = 345
IY1 = 50
IX2 = X1 + SCL
IY2 = Y1 + SCL
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X1
IY1 = Y1 + Y2
IX2 = X1 + SCL
IY2 = Y1 + Y2
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
C ------------------------------------------------
DO 500 I = 1, 16
C
IC = I + 1
IX1 = X1 + ( 1 - GVX( I))* SCL
IY1 = Y1 + Y2 + GV ( I)* Y2
IX2 = X1 + ( 1 - GVX( IC))* SCL
IY2 = Y1 + Y2 + GV ( IC)* Y2
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, IICR )
C
500 CONTINUE
C ------------------------------------------------
C
C ----- Computed V ---
C
C ------------------------------------------------
DO 600 I = 2, II-1
C
N = JM / 2
VP = - ( V( I,N) + V( I,N+1))/ 2.D0
J = ( II/2 - 1)* IM1 + I - 1
XP = ( X( J) + X( J+1))/ 2.D0
IX1= X1 + XP* SCL
IY1= Y1 + Y2 + VP* Y2
C
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
600 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PRSB ( UE, VE, U, V, X, Y, XY, PE, IX,
1 IY, NNP, NEL, ICR, KKK, NBC )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
COMMON /DEF/ LOP, DT, RE, FSM, WEI
C
CHARACTER*8 FSM
INTEGER*2 ICR
C
DIMENSION UE( NEL), VE( NEL), X( NEL), Y( NEL), PE( NNP),
1 XY( 2,NNP), U ( IX,IY), V( IX,IY)
C
COMMON /HIST/ WMX( 5000), DMX( 5000)
C
ICR = 15
C
WRITE(6,2000) FSM, WEI, RE, DT, LOP, NNP, NEL, IX, IY, NBC
2000 FORMAT(/,1X,A8,' * WEI =',F7.2,' * Re =',F8.3,' * DT =',F7.4,
1 ' * Loop =',I4,/,' * NNP =',I4,
2 ' * NEL =',I4,' IX* IY=',I3,' X',I3,' * NBC =',I2)
C
READ(7,1000) ( X ( I), Y ( I), I = 1, NEL)
READ(7,1000) ( UE( I), VE( I), I = 1, NEL )
1000 FORMAT(10F8.4)
READ(7,1200) ( WMX( I), I= 1, LOP-1 )
READ(7,1200) ( DMX( I), I= 1, LOP-1 )
1200 FORMAT(10F8.2)
READ(7,1000) ( ( XY( I,J), I = 1, 2), J = 1, NNP )
READ(7,1000) ( PE( I), I = 1, NNP )
READ(7,1000) ( ( U( I,J), I = 1, IX), J = 1, IY )
READ(7,1000) ( ( V( I,J), I = 1, IX), J = 1, IY )
C
WRITE(6,*) ' '
WRITE(6,*) ' * KKK = 1 ( Out-1 ), = 2 ( Out-2 )'
C
READ(5,*) KKK
WRITE(6,2100) KKK
2100 FORMAT(' * KKK =',I2)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,20)
WRITE(6,21)'*****************************************************'
WRITE(6,21)' '
WRITE(6,21)' ##### ### #### ### # ### '
WRITE(6,21)' # # # # # # # ## # # '
WRITE(6,21)' #### # # # #### # # # # # '
WRITE(6,21)' # ## # # # # ###### ## '
WRITE(6,21)' ## ##### #### ### # ##### '
WRITE(6,21)' '
WRITE(6,21)' Graphics Output '
WRITE(6,21)' '
WRITE(6,21)' ( Version 4.0 ) '
WRITE(6,21)' '
WRITE(6,21)' for Two Dimensional Viscous Fluid Flow '
WRITE(6,21)' '
WRITE(6,21)' by the 4th-ordered and the 2nd-ordered FDM '
WRITE(6,21)' '
WRITE(6,21)' ( ** Equal mesh size : DX = DY = DS ** ) '
WRITE(6,21)' '
WRITE(6,21)' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,21)' '
WRITE(6,21)'*****************************************************'
C
20 FORMAT(/)
21 FORMAT(A70)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE MSHFG
C
INTEGER*2 ICR, XIT, YIT, IWD, IWE, IW1, IW2
C
ICR = 15
C
IWD = 150
XIT = 90
YIT = 100
C
IW1 = XIT + IWD
IW2 = YIT + IWD
C
CALL RECTANGLE@( XIT, YIT, IW1, IW2, ICR )
C
IWE = IWD / 10
C
C -----------------------------------------------------
DO 100 I = 1, 9
C
IW2 = YIT + IWE* I
CALL DRAW_LINE@( XIT, IW2, IW1, IW2, ICR )
C
100 CONTINUE
C -----------------------------------------------------
C
IW2 = YIT + IWD
C -----------------------------------------------------
DO 200 J = 1, 9
C
IW1 = XIT + IWE* J
CALL DRAW_LINE@( IW1, YIT, IW1, IW2, ICR )
C
200 CONTINUE
C -----------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VRPT ( RMAX, ICR )
C
IMPLICIT REAL*8 (A-H, O-Z)
C
INTEGER*2 ICR, IX0, IY0, IX1, IY1
C
C ----- Attention N.LE.500 ---
C
DIMENSION RMAX( 5000)
INTEGER*2 IVR( 5000), IXX( 5000)
C
COMMON /DEF/ LOP, DT, RE, FSM, WEI
C
CHARACTER*8 FSM
CHARACTER*80 STR0
CHARACTER*80 STR01
CHARACTER* 4 STR1( 3)
CHARACTER* 3 STR2( 4)
C
DATA STR0
1/'*** Results for F2D-D42 *** ( Max_rel_var_Vorticity & Velocity'/
DATA STR01
1 /' Meshes & Stream fUEction contours & Velocity distribution'/
DATA STR1 /'0 ', 'LP/2', 'LOP'/
DATA STR2 /'400','300','200','100'/
C
IF ( ICR.NE.15) GO TO 90
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 120, 10, ICR )
CALL DRAW_TEXT@( STR01,120, 30, ICR )
CALL DRAW_TEXT@( STR1( 1), 45, 450, ICR )
CALL DRAW_TEXT@( STR1( 2), 315, 450, ICR )
CALL DRAW_TEXT@( STR1( 3), 605, 450, ICR )
C
CALL DRAW_TEXT@( STR2( 1), 20, 50, ICR )
CALL DRAW_TEXT@( STR2( 2), 20, 150, ICR )
CALL DRAW_TEXT@( STR2( 3), 20, 250, ICR )
CALL DRAW_TEXT@( STR2( 4), 20, 350, ICR )
C
IX0 = 50
IY0 = 50
IX1 = IX0 + 580
IY1 = IY0 + 400
C
CALL RECTANGLE@( IX0, IY0, IX1, IY1, ICR )
C
ICR = 11
90 CONTINUE
C
RIWD = 580.D0/ DFLOAT( LOP)
C
C -------------------------------------------
DO 100 I = 1, LOP
C
IXX( I) = IX0 + DFLOAT( I)* RIWD
IVR( I) = IY0 + 400 - RMAX( I)
C
100 CONTINUE
C -------------------------------------------
C
CALL POLYLINE@( IXX, IVR, LOP, ICR )
C -------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTN ( X, Y, UE, VE, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
CHARACTER*8 FSM
C
INTEGER*2 X0, Y0, X1, Y1, X2, Y2, SCX, SCY, GX, GY
C
COMMON /DEF/ LOP, DT, RE, FSM, WEI
C
DIMENSION X( NEL), Y( NEL), UE( NEL), VE( NEL)
C
X0 = 450
Y0 = 100
SCX = 150
SCY = 150
C
SCL = 50
X1 = X0 + SCX
Y1 = Y0 + SCY
C
CALL RECTANGLE@( X0, Y0, X1, Y1, 15 )
C
TL = DT* DFLOAT( LOP)
C
C ------------------------------------------------
DO 100 I = 1, NEL
C
GX = X0 + X( I)* SCX
GY = Y0 + SCY - Y( I)* SCY
C
IF ( GX.EQ.X0 ) GO TO 100
IF ( GY.EQ.Y0 ) GO TO 100
IF ( GX.EQ.X0 + SCX) GO TO 100
IF ( GY.EQ.Y0 + SCY) GO TO 100
C
X1 = GX
Y1 = GY
GGX = UE( I)* SCL
GGY = VE( I)* SCL
X2 = GX + GGX
Y2 = GY - GGY
C
CALL DRAW_LINE@( X1, Y1, X2, Y2, 11 )
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE CCNT ( X, CTV, NODE, NNP, NEL )
C
C ----- Plot Contour Lines for Rectangular Meshes ---
C ( Linear Interpolation )
C
IMPLICIT REAL*8(A-H,O-Z)
C
INTEGER*2 X00, Y00, X1, Y1, X2, Y2, SCLX, SCLY, X0, Y0
C
DIMENSION X( 2,NNP), CTV( NNP), NODE( NEL,4)
DIMENSION PX( 4), PY( 4), HVE( 100)
C
X00 = 270
Y00 = 100
C
C ----- SCLX, SCLY : Width ---
C
SCLX = 150
SCLY = 150
X1 = X00 + SCLX
Y1 = Y00 + SCLY
C
CALL RECTANGLE@ ( X00, Y00, X1, Y1, 15 )
C
NQ = DSQRT( DFLOAT( NNP+1)) - 1.D0
L = 0
K = 0
C
C ----------------------------------
DO 100 J = 1, NQ
C
DO 150 I = 1, NQ
C ---------------------------------
K = K + 1
L = L + 1
NODE( K,1) = L
NODE( K,2) = L + 1
NODE( K,3) = L + 1 + NQ + 1
NODE( K,4) = L + NQ + 1
C
150 CONTINUE
L = L + 1
C
100 CONTINUE
C ----------------------------------
X0 = 0
Y0 = 0
C
C ----------------------------------------------------------
DO 200 I = 1, NQ* NQ
C
X1 = ( X( 1, NODE( I,4)) + X0 )* SCLX + X00
Y1 = ( X( 2, NODE( I,4)) + Y0 )* SCLY + Y00
C
C ------------------------------------------------
DO 250 J = 1, 4
C
X1 = ( X( 1, NODE( I,J)) + X0 )* SCLX + X00
Y1 = ( X( 2, NODE( I,J)) + Y0 )* SCLY + Y00
C
250 CONTINUE
C ------------------------------------------------
C
200 CONTINUE
C ---------------------------------------------------------
C
RMI = 100.D0
RMX = 0.D0
C
DO 300 I = 1, NNP
C
IF ( CTV( I).LE.RMI ) RMI = CTV( I)
IF ( CTV( I).GE.RMX ) RMX = CTV( I)
300 CONTINUE
C
WRITE(6,2000) RMI, RMX
2000 FORMAT(' ** Rmin.=',F8.3,5X,' Rmax.=',F8.3 )
C
NM = 12
C ============
C
DRM = ( RMX - RMI)/ DFLOAT( NM - 1)
C
DO 400 I = 1, NM
C
HVE( I) = RMI + DRM* DFLOAT( I-1)
400 CONTINUE
C
C ----- Plot Contour Lines ---
C
PXV = 0.D0
PYV = 0.D0
C
C --------------------------------------------
DO 500 L = 1, NEL
C
DO 550 IH = 2, NM-1
C
HV = HVE( IH)
C ******************
C
ICL = 0
C
C -----------------------------------------
DO 570 J = 1, 4
C
I1 = NODE( L, J)
C
IF ( J .NE. 4 ) THEN
I2 = NODE ( L, J+1)
ELSE
I2 = NODE ( L, 1)
END IF
C
IF ( ABS ( CTV( I1) - CTV( I2)).LE. 0.0001) THEN
IF ( ABS ( CTV( I1) - HV ) .LE.0.0001) THEN
C
ICL = ICL + 1
PX( ICL ) = X( 1,I1)
PY( ICL ) = X( 2,I1)
END IF
C
ELSE IF ( ( ( CTV( I1) .LE. HV ) .AND. ( HV .LT. CTV( I2)) )
1 .OR. ( ( CTV( I1) .GE. HV ) .AND. ( HV .GT. CTV( I2)) ))
2 THEN
C
ICL = ICL + 1
RT = ( HV - CTV( I1)) / ( CTV( I2) - CTV( I1))
C
PX( ICL ) = ( X( 1,I1) + ( X( 1,I2) - X( 1,I1))* RT )
PY( ICL ) = ( X( 2,I1) + ( X( 2,I2) - X( 2,I1))* RT )
END IF
C
570 CONTINUE
C ----------------------------------------------
C
IF ( ICL .GE. 2) THEN
C -----------------------------------------------
DO 580 J = 1, ICL
C
JM = J
D2 = ( PXV-PX( J))**2 + ( PYV-PY( J))**2
C
C -------------------------------------
DO 590 JJ = J + 1, ICL
C
D2V = ( PXV - PX( JJ))**2 + ( PYV - PY( JJ))**2
C
IF ( D2V .LT. D2) THEN
JM = JJ
D2 = D2V
END IF
C
590 CONTINUE
C --------------------------------------
C
PXV = PX( JM)
PYV = PY( JM)
PX( JM) = PX( J)
PY( JM) = PY( J)
PX( J) = PXV
PY( J) = PYV
C
580 CONTINUE
C ------------------------------------------------
PXV = PX( ICL)
PYV = PY( ICL)
C
END IF
C
IF ( ICL .GT.0) THEN
C
X1 = ( X0 + PX( 1))* SCLX + X00
Y1 = ( Y0 + PY( 1))* SCLY + Y00
C
DO 595 J = 2, ICL
C
X2 = ( X0 + PX( J))* SCLX + X00
Y2 = ( Y0 + PY( J))* SCLY + Y00
CALL DRAW_LINE@ ( X1, Y1, X2, Y2, 15 )
C
595 CONTINUE
C
IF ( ICL.EQ.4) X2 = ( X0+PX( 1))* SCLX + X00
IF ( ICL.EQ.4) Y2 = ( Y0+PY( 1))* SCLX + Y00
IF ( ICL.EQ.4 ) CALL DRAW_LINE@ ( X1, Y1, X2, Y2, 15 )
C
END IF
C
550 CONTINUE
C
500 CONTINUE
C -------------------------------------------
C
RETURN
END
C **********************************************************************