戻る

C ******************************************************************** C
C C
PROGRAM MAIN
C C
C Graphics for F2D-BSG C
C C
C ( Fourth 0rdered Finite Difference Method ) C
C C
C ( Version 3.8 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C ******************************************************************** C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ***** Attention *** No. of Var. of Vorticity Data = < 100,000 ***
C
C ----- DATA ( A ) : Mesh size -----------------------------------------
C
PARAMETER ( IX= 181, IY= 61,
C
*** PARAMETER ( IX= 45, IY= 15,
C
C ---------------------------------------------------------------------
C
1 IX1 = IX-1, IY1 = IY-1, IXY = IX*IY, IXY1 = IX1*IY1 )
C
C ---------------------------------------------------------------------
C
C ----- IXY: NNP, IXY1: NEL ---
C
INTEGER*2 NDE
C
DIMENSION X ( IXY), Y ( IXY), PSV( IXY), VRV( IXY),
1 XX( IXY), YY( IXY), HX ( IXY), HY ( IXY),
2 NDE( IXY1,4)
C
LOGICAL*1 PEP( IXY,4), EEP( IXY,4)
C
DIMENSION RMA( 100000), II( 100000), JJ( 100000), X1R( 100000)
C
C ----- DATA ( B ) ----------------------------------------------------
C
OPEN ( 8, FILE='BSG_OUT_R800.DATA' )
OPEN (17, FILE='BSG_VXR_R800.DATA' )
C
*** OPEN ( 8, FILE='BSG_OUT_R800_F.DATA' )
*** OPEN (17, FILE='BSG_VXR_R800_F.DATA' )
C
*** OPEN ( 8, FILE='BSG_OUT_R1000_F.DATA' )
*** OPEN (17, FILE='BSG_VXR_R1000_F.DATA' )
C
C ---------------------------------------------------------------------
C
CALL TTLE
C
CALL SBINP ( RMA, II, JJ, XX, YY, PEP, EEP, X, Y,
1 NDE, PSV, VRV, X1R, LOP, IX1, IY1, RLL, IX,
2 IS, JS, IE, IXY, IXY1, NEL )
C
CALL VGA@
C
CALL VRPLT ( RMA, II, JJ, LOP, 1 )
C
CALL MESHF ( IX1, IY1, RLL, IXY )
C
CALL TRBDY ( X, Y, XX, YY, HX, HY, PEP, EEP, IX, IE, JS, IXY, 1 )
C
CALL TRBDY ( X, Y, XX, YY, HX, HY, PEP, EEP, IX, IE, JS, IXY, 2 )
C
C ----- Contour for Stream Function ---
C
CALL TRCNT ( X, Y, NDE, PSV, IX, IS, JS, IE, IXY, NEL, 1 )
C
C ----- Contour for Vorticity ---
C
CALL TRCNT ( X, Y, NDE, VRV, IX, IS, JS, IE, IXY, NEL, 2 )
C
C ----- Plot of Xir ---
C
CALL VRPLT ( X1R, II, JJ, LOP, 2 )
C
STOP
END
C **********************************************************************
C
SUBROUTINE MESHF ( IX1, IY1, RLL, IXY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 IP1, IP2, IP3, IP4
C
NEL = IX1*IY1
C
WRITE(6,*) ' '
WRITE(6,2000) IX1, IY1, NEL
2000 FORMAT(' * Meshes :',I4,' X', I4,' NEL =', I6 )
WRITE(6,*) ' '
C
X0 = 75.D0
Y0 = 10.D0
C
XL = RLL*15.D0
C
C ***** Attention : Area size: 2 X RLL ***
C
C Actual height : 2 ==> YL = 60
C
YL = 60.D0
C
DXL = XL/ DFLOAT( IX1)
DYL = YL/ DFLOAT( IY1)
IF ( IX1.LE.50 ) GO TO 333
C ***** Attention ****************
C
WRITE(6,*) ' *** Attention: Sparse Plotting for Fig. MESH ***'
WRITE(6,*) ' '
C
C ----- Attention ---- Sparse Plotting for meshes in X-direction ----
C *****************
C
DO 100 I = 1, IX1+1, 2
C ----- Attention ---
C
IP1 = IDNINT( X0 + DXL* DFLOAT( I-1))
IP2 = IDNINT( Y0 + YL )
IP3 = IP1
IP4 = IDNINT( Y0 + YL + YL )
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, 15 )
C
100 CONTINUE
C
C ----- Attention ---- Sparse Plotting for meshes in Y-direction ----
C *****************
C
C ------------------------------------------------
DO 200 J = 1, IY1+1, 2
C ----- Attention ---
C
IP1 = IDNINT( X0 )
IP2 = IDNINT( Y0 + DYL* DFLOAT( J-1) + YL )
IP3 = IDNINT( X0 + XL )
IP4 = IP2
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, 15 )
C
200 CONTINUE
C ------------------------------------------------
RETURN
C
333 CONTINUE
C
C ------------------------------------------------
DO 300 I = 1, IX1+1
C
IP1 = IDNINT( X0 + DXL* DFLOAT( I-1))
IP2 = IDNINT( Y0 + YL )
IP3 = IP1
IP4 = IDNINT( Y0 + YL + YL )
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, 15 )
C
300 CONTINUE
C -----------------------------------------------
DO 400 J = 1, IY1+1
C
IP1 = IDNINT( X0 )
IP2 = IDNINT( Y0 + DYL* DFLOAT( J-1) + YL )
IP3 = IDNINT( X0 + XL )
IP4 = IP2
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, 15 )
C
400 CONTINUE
C -----------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE SBINP ( RMA, II, JJ, XX, YY, PEP, EEP, X, Y,
1 NDE, PSV, VRV, X1R, LOP, IX1, IY1, RLL, IX,
2 IS, JS, IE, IXY, IXY1, NEL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION XX ( IXY), YY( IXY), X( IXY), Y( IXY), PSV( IXY),
1 VRV( IXY )
C
DIMENSION RMA( 100000), II( 100000), JJ( 100000), X1R( 100000)
C
INTEGER*2 NDE
DIMENSION NDE( IXY1,4)
C
LOGICAL*1 PEP( IXY,4), EEP( IXY,4)
C
COMMON /CNST/ RE
C
ICR = 15
C
READ(8,1000) RE, DT, LOPN, NZU, NPR
1000 FORMAT(F9.2,F6.3,I6,I2,I3)
C
WRITE(6,*) ' '
WRITE(6,2000) RE, DT, LOPN, NZU, NPR
2000 FORMAT(' ** F2D-BSG ** Re =',F9.2,' DT =',F6.3,' LOP =',I6,
1 ' NZU =',I2,' NPR =',I3 )
C
READ(8,1100) LOP, LPFL
1100 FORMAT(2I6)
WRITE(6,2100) LOP
2100 FORMAT(/,' *** NPLT =',I6,/)
C
C ----- Last Loop --- Var. of Vorticity etc. ---
C
DO 100 I = 1, LOP
C
READ(17,1200) IL, RMA( I), II( I), JJ( I), X1R( I)
WRITE(6,2200) IL, RMA( I), II( I), JJ( I), X1R( I)
C
100 CONTINUE
C
1200 FORMAT(I6,F9.3,2I4,F6.3 )
2200 FORMAT(' * LOP =', I6,' Var.of(W)=', F9.3,'(%) Max. at I, J =',
1 2I4,' X1=',F6.3 )
C
READ(8,1300) IX1, IY1, RLL, IX, IS, JS, IE, IXY, NEL
1300 FORMAT(2I4,F7.2,4I4,2I6)
C
WRITE(6,2300) LPFL
2300 FORMAT(/,' *** Final LOP =', I6)
C
READ(8,1400) ( XX( I), I = 1, IXY)
READ(8,1400) ( YY( I), I = 1, IXY)
READ(8,1500) ( ( PEP( I,J), J = 1, 4), I = 1, IXY)
READ(8,1500) ( ( EEP( I,J), J = 1, 4), I = 1, IXY)
1500 FORMAT(80L1)
READ(8,1400) ( X( I), I = 1, IXY)
READ(8,1400) ( Y( I), I = 1, IXY)
1400 FORMAT(10F8.4)
READ(8,1600) ( ( NDE( I,J), J =1, 4), I = 1, NEL )
1600 FORMAT(20I6)
READ(8,1400) ( PSV( I), I = 1, IXY)
READ(8,1400) ( VRV( I), I = 1, IXY)
C
C ----- Attention ------------
C
DO 200 I = 1 , IXY
C
Y( I) = 1.D0 - Y( I)
200 CONTINUE
C ----------------------------
NEL = IX1*IY1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TRBDY ( X, Y, XX, YY, HX, HY, PEP, EEP, IX, IE,
1 JS, IXY, KY )
C
C To plot Boundaries
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X ( IXY), Y ( IXY), XX( IXY), YY( IXY), HX( IXY),
1 HY( IXY)
C
LOGICAL*1 PEP( IXY,4), EEP( IXY,4)
C
C ----- Attention ---
C
INTEGER*2 X0, Y0, IP1, IP2, IP3, IP4, IP5, IWK
C
C ----- Attention ---
C
C ----------------------------
DO 100 I = 1, IXY
C
HX( I) = 0.15D0
HY( I) = 0.15D0
C
100 CONTINUE
C ----------------------------
C
IWK = 15
X0 = 75
C
IF ( KY.EQ.1) Y0 = 180
IF ( KY.EQ.2) Y0 = 320
C
SCLX = 15.D0
SCLY = 50.D0
C
C ----- Boundary conditions ----------------------
DO 200 I = 1, IXY
C
IF ( PEP( I,1)) GO TO 70
IF ( PEP( I,2)) GO TO 80
IF ( PEP( I,3)) GO TO 90
GO TO 60
C
C ----- Continuative Boundary ---
C
70 CONTINUE
C
C --------------------------------------
DO 300 J = 1, 4
C
IF ( EEP( I,J)) GO TO 71
GO TO 300
C
71 IF ( PEP( I,4)) THEN
PX = X0 + X( I)* SCLX
PY = Y0 + Y( I)* SCLY
ENDIF
C
300 CONTINUE
C --------------------------------------
GO TO 200
C
C ----- Slip boundary ---
C
80 CONTINUE
C
PX = X0 + X( I)* SCLX
PY = Y0 + Y( I)* SCLY
C
IP1 = IDNINT( PX - XX( I)* SCLX / 2.D0 )
IP2 = IDNINT( PY - YY( I)* SCLY / 2.D0 )
IP3 = IDNINT( PX + XX( I)* SCLX / 2.D0 )
IP4 = IDNINT( PY + YY( I)* SCLY / 2.D0 )
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, IWK )
C
WRITE(6,*) '1', PX - XX( I)*SCLX/2, PY + YY( I)*SCLY/2
GO TO 200
C
C ----- Move Boundary ---
C
90 CONTINUE
C
PX = X0 + X( I)* SCLX
PY = Y0 + Y( I)* SCLY
C
IP1 = IDNINT( PX - XX( I)* SCLX / 2.D0 )
IP2 = IDNINT( PY - YY( I)* SCLY / 2.D0 )
IP3 = IDNINT( PX + XX( I)* SCLX / 2.D0 )
IP4 = IDNINT( PY + YY( I)* SCLY / 2.D0 )
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, IWK )
C
GO TO 200
C
C ----- No-slip Boundary ---
C
60 CONTINUE
PX = X0 + XX( I)* SCLX
PY = Y0 + YY( I)* SCLY
C
IJ = ( I-1)/IX
IF ( IJ .EQ. (JS-1)) GO TO 65
C
IP1 = IDNINT( PX)
IP2 = IDNINT( PY - HY( I)* SCLY / 2.D0 )

IP3 = IDNINT( PX + HX( I)* SCLX / 2.D0 )
IP4 = IDNINT( PY + HY( I)* SCLY / 2.D0 )

IP5 = PX - HX( I)* SCLX / 2.D0
C
C ----- Inlet boundary ( Omit ) ---
C
IJ = ( I-IE) / IX * IX
C
65 CONTINUE
C
C ----- Plot of Top & Bottom ***
C
IF ( DABS( YY( I) - YY( 1) ).GT.0.001D0 .AND.
1 DABS( YY( I) - YY( IXY)).GT.0.001D0 ) GO TO 66
C
IP1 = IDNINT( PX - HX( I)* SCLX / 2.D0 )
IP2 = IDNINT( PY - HY( I)* SCLY / 2.D0 )
C
IP3 = IDNINT( PX + HX( I)* SCLX / 2.D0 )
IP4 = IDNINT( PY + HY( I)* SCLY / 2.D0 )
C
CALL DRAW_LINE@ ( IP1, IP2, IP3, IP4, IWK )
C
66 CONTINUE

200 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TRCNT ( X, Y, NDE, CTV, IX, IS, JS, IE, IXY, NEL, KY )
C
C ----- To plot contour lines for rectangular meshes ---
C
C ( Linear Interpolation )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 NDE
INTEGER*2 ICR, IWK1, IWK2, IWK3, IWK4
C
DIMENSION X ( IXY), Y ( IXY), NDE( NEL,4), CTV( IXY), PX( 4),
1 PY( 4), HVA( 50)
C
RMI= 100.D0
RMX = 0.D0
C
C ------------------------------------------------
DO 100 I = 1, IXY
C
IF ( CTV( I).LE.RMI ) RMI = CTV( I)
IF ( CTV( I).GE.RMX ) RMX = CTV( I)
C
100 CONTINUE
C ------------------------------------------------
WRITE(6,1231) KY, RMI, RMX
1231 FORMAT(' ** KY =',I2,' Rmin.=',F8.3,5X,' RMA.=',F8.3 )
C
NM = 31
C ===========
C
DRM = ( RMX - RMI)/ DFLOAT( NM - 1)
DO 200 I = 1, NM
C
HVA( I) = RMI + DRM* DFLOAT( I-1)
C
200 CONTINUE
C
IF ( KY.EQ.0) GO TO 555
C
C ------- For Reference with JSME-Paper (2001-9) ---
C
IF ( KY.NE.1) GO TO 333
C
C ----- For Stream Function ---
C
ICR = 15
X0 = 75.D0
Y0 = 230.D0
C
NM = 20
C
HVA( 1) = -1.12D0
HVA( 2) = -1.10D0
HVA( 3) = -1.08D0
HVA( 4) = -1.06D0
HVA( 5) = -1.04D0
C
HVA( 6) = -1.02D0
HVA( 7) = -1.D0
HVA( 8) = -0.8D0
HVA( 9) = -0.6D0
HVA(10) = -0.4D0
C
HVA(11) = -0.2D0
HVA(12) = 0.D0
HVA(13) = 0.2D0
HVA(14) = 0.4D0
HVA(15) = 0.6D0
C
HVA(16) = 0.8D0
HVA(17) = 0.96D0
HVA(18) = 1.D0
HVA(19) = 1.008D0
HVA(20) = 1.016
GO TO 555
C
333 CONTINUE
C
C ----- For Vorticity ---
C
ICR = 15
X0 = 75.D0
Y0 = 370.D0
C
NM = 10
C
HVA( 1) = 10.D0
HVA( 2) = 8.D0
HVA( 3) = 6.D0
HVA( 4) = 4.D0
HVA( 5) = 2.D0
C
HVA( 6) = -2.D0
HVA( 7) = -4.D0
HVA( 8) = -6.D0
HVA( 9) = -8.D0
HVA(10) = 0.00000001D0
C
555 CONTINUE
C
C ----- Attention ---
C
SCLX = 15.D0
SCLY = 50.D0
C
PXV = 0.D0
PYV = 0.D0
C
C ----------------------------------------------------------
DO 700 L = 1, NEL
C
C ------------------------------------------------
DO 800 IH = 1, NM
C
HV = HVA( IH)
ICL = 0
C
C --------------------------------------
DO 900 J = 1, 4
C
I1 = NDE( L,J)
IF ( J .NE.4) THEN
I2 = NDE( L,J+1)
ELSE
I2 = NDE( L,1)
END IF
C
C ----- Attention ---
C
IF ( DABS ( CTV( I1) - CTV( I2)).LE.1.D-30 ) THEN
C ===========
IF ( DABS ( CTV( I1) - HV ) .LE.1.D-30 ) THEN
C ===========
ICL = ICL + 1
PX( ICL ) = X( I1)* SCLX
PY( ICL ) = Y( I1)* SCLY
END IF
C
ELSE IF ( ( ( CTV( I1) .LE. HV ).AND.( HV .LT. CTV( I2)) )
1 .OR. ( ( CTV( I1) .GE. HV )
2 .AND.( HV .GT. CTV( I2 )) )) THEN
C
ICL = ICL + 1
C
RT = ( HV - CTV( I1)) / ( CTV( I2) - CTV( I1))
PX( ICL ) = ( X( I1) + ( X( I2) - X( I1))* RT )* SCLX
PY( ICL ) = ( Y( I1) + ( Y( I2) - Y( I1))* RT )* SCLY
C
END IF
C
900 CONTINUE
C --------------------------------------
C
IF ( ICL.GE.2 ) THEN
C
C --------------------------------------
DO 400 J = 1, ICL
C
JM = J
D2 = ( PXV - PX( J))** 2 + ( PYV - PY( J))** 2
C
C ----------------------------
DO 500 JJ = J + 1, ICL
C
D2V = ( PXV - PX( JJ))** 2 + ( PYV - PY( JJ))** 2
IF ( D2V .LT. D2 ) THEN
JM = JJ
D2 = D2V
END IF
C
500 CONTINUE
C ----------------------------
C
PXV = PX( JM)
PYV = PY( JM)
PX( JM) = PX( J)
PY( JM) = PY( J)
PX( J) = PXV
PY( J) = PYV
C
400 CONTINUE
C --------------------------------------
C
PXV = PX( ICL)
PYV = PY( ICL)
END IF
C
IF ( ICL.GT.0 ) THEN
IWK1 = IDNINT( X0 + PX( 1))
IWK2 = IDNINT( Y0 + PY( 1))
C
C --------------------------------------
DO 600 J = 2, ICL
C
IWK3 = IDNINT( X0 + PX( J))
IWK4 = IDNINT( Y0 + PY( J))
C
CALL DRAW_LINE@ ( IWK1, IWK2, IWK3, IWK4, 11 )
C
600 CONTINUE
C --------------------------------------
C
END IF
C
800 CONTINUE
C ------------------------------------------------
C
700 CONTINUE
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
WRITE(6,2100) '*************************************************'
WRITE(6,2100) ' '
WRITE(6,2100) ' Graphics for F2D-BSG '
WRITE(6,2100) ' '
WRITE(6,2100) ' Backstep Fluid Flow Analysis '
WRITE(6,2100) ' '
WRITE(6,2100) ' by 4th ordered FDM '
WRITE(6,2100) ' '
WRITE(6,2100) ' ( Version : 3.8 ) '
WRITE(6,2100) ' '
WRITE(6,2100) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,2100) ' '
WRITE(6,2100) '*************************************************'
2100 FORMAT(A55)
C
RETURN
END
C **********************************************************************
C
SUBROUTINE VRPLT ( RMA, II, JJ, LOP, KY )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
C ----- Attention N.LE.30000 ---
C
DIMENSION RMA( 100000), II( 100000), JJ( 100000)
C
INTEGER*2 IVAR( 100000), IXX( 100000)
C
COMMON / CNST/ RE
C
CHARACTER*70 STR0
CHARACTER*6 STR1( 3)
CHARACTER*3 STR3( 4)
C
DATA STR0/'*** Results of F2D-BSG *** '/
DATA STR1/'0 ','LOP/2',' LOP'/
DATA STR3/' 40',' 30',' 20',' 10'/
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@ ( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR0, 70, 20, ICR )
C
IF ( DABS( RE - 1000.D0).LE.0.01D0 )
C
1 CALL DRAW_TEXT@ ( '( Re = 1000 )', 370, 20, ICR )
C
IF ( DABS( RE - 800.D0).LE.0.01D0 )
C
1 CALL DRAW_TEXT@ ( '( Re = 800 )', 370, 20, ICR )
C
ICR = 15
C
CALL DRAW_TEXT@ ( STR1( 1), 45, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 2), 300, 450, ICR )
CALL DRAW_TEXT@ ( STR1( 3), 580, 450, ICR )
C
IF ( KY.NE.1 ) GO TO 133
C
CALL DRAW_TEXT@ ( '1000', 15, 38, ICR )
CALL DRAW_TEXT@ ( ' 100', 15, 118, ICR )
CALL DRAW_TEXT@ ( ' 10', 15, 198, ICR )
CALL DRAW_TEXT@ ( ' 1', 15, 278, ICR )
CALL DRAW_TEXT@ ( ' 0.1', 15, 358, ICR )
CALL DRAW_TEXT@ ( '0.01', 15, 438, ICR )
C
CALL DRAW_TEXT@ ( '(%)', 15, 245, ICR )
C
GO TO 125
C
133 CONTINUE
CALL DRAW_TEXT@ ( STR3( 1), 604, 38, ICR )
CALL DRAW_TEXT@ ( STR3( 2), 604, 138, ICR )
CALL DRAW_TEXT@ ( STR3( 3), 604, 238, ICR )
CALL DRAW_TEXT@ ( STR3( 4), 604, 338, ICR )
C
125 CONTINUE
IX0 = 50
IY0 = 40
C
CALL RECTANGLE@( IX0, IY0, IX0+580, IY0+400, 15 )
C
RIWD = 580.D0/ DFLOAT( LOP)
C
IF ( KY.EQ.1) THEN
C
DO 150 I = 1, LOP
C
RMA( I) = DLOG10( RMA( I))
C ===============================
C
150 CONTINUE
C
SCL = 80.D0
C
END IF
C
IF ( KY.EQ.2) SCL = 10.D0
C
DO 100 I = 1, LOP
C
IXX ( I) = IX0 + IDNINT( DFLOAT( I)* RIWD )
IVAR( I) = 280 - IDNINT( RMA( I)* SCL )
IF ( KY.EQ.2 ) IVAR( I) = 440 - IDNINT( RMA( I)* SCL )
C
100 CONTINUE
C
IF ( KY.EQ.1 ) ICR = 14
IF ( KY.EQ.2 ) ICR = 12
C
CALL POLYLINE@ ( IXX, IVAR, LOP, ICR )
C
RETURN
END
C **********************************************************************