戻る

C
PROGRAM MAIN
C *******************************************************************C
C C
C Thermal Fluid Analysis Software by Finite Element Method C
C C
C ( Rectangular-Linear Element ) C
C C
C FIG - F 2 D - F E R ( SV ) C
C C
C ( Version : 3.5 ) C
C C
C Copyright : Yasuhiro MATSUDA C
C C
C *******************************************************************C
C
IMPLICIT REAL*8 (A-H,O-Z)
C
C ----- Attention : ( 10 x 10 ) ---
C
DIMENSION X ( 121), Y( 121), UN( 121), VN( 121), PSI( 121),
1 VOR( 121)
C
INTEGER*2 NDE
DIMENSION NDE( 100,4)
C
COMMON /CNST/ RE, KPLT
COMMON /TIMC/ DT, LPM, LP1, LOP, ALP
C
C ----- Attention --- < = 500 ---
COMMON /VRPT/ VARV( 500), VRUV( 500), VARP( 500)
C
OPEN ( 8, FILE='GH100.DAT')
OPEN ( 10, FILE='F2D_FER_PLT.DAT')
C
READ(10,1000) LOP, DT, RE, NNP, NEL
1000 FORMAT(I4, F6.3, F8.1, 2I4)
C
C -------------------------------------------------------------
C
CALL PRPS ( PSI, VOR, UN, VN, X, Y, NDE, NNP, NEL )
C
C -------------------------------------------------------------
IF ( KPLT.NE.1) GO TO 900
C
C -------------------------------------------
C
CALL PMSH ( X, Y )
C
C -------------------------------------------
C
CALL PLTN ( X, Y, UN, VN, NNP )
C
C -------------------------------------------
C
CALL PLVR ( VRUV, LOP, 1 )
C
C -------------------------------------------
C
CALL PLVR ( VARP, LOP, 4 )
C
C -------------------------------------------
C
CALL PLVR ( VARV, LOP, 3 )
C
C -------------------------------------------
C
C ----- Attention : Change Direction ---
C
DO 100 I = 1, NNP
C
SX = X( I)
X( I) = 1.D0 - Y( I)
Y( I) = SX
C
100 CONTINUE
C --------------------------------------------------
C
CALL PCNT ( X, Y, NDE, PSI, NNP, NEL, 0 )
C
C --------------------------------------------------
GO TO 990
C
900 CONTINUE
C ----------------------------------------
C
CALL CMPT ( UN, VN, X, Y, NNP )
C
C ----------------------------------------
990 CONTINUE
C
CLOSE ( 8)
CLOSE (10)
C
STOP
END
C **********************************************************************
C
SUBROUTINE PRPS ( PSI, VOR, UN, VN, X, Y, NDE, NNP, NEL)
C
IMPLICIT REAL*8 (A-H, O-Z)
C
INTEGER*2 ICR, NDE
C
DIMENSION NDE( NEL,4)
C
DIMENSION PSI( NNP), VOR( NNP), UN( NNP), VN( NNP), X( NNP),
1 Y ( NNP)
C
COMMON /CNST/ RE, KPLT
COMMON /TIMC/ DT, LPM, LP1, LOP, ALP
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
COMMON /VRPT/ VARV( 500), VRUV( 500), VARP( 500)
C
C --------------------
C
CALL TTLE
C
C --------------------
C
WRITE(6,2000) RE, DT, LOP, NNP, NEL
2000 FORMAT(/,' * Re =',F6.1,' DT =',F8.4,' Loop =',I4,' NNP =',
1 I4,' NEL =',I4 )
C
ICR = 15
C
READ(10,1000) ( X( I), Y( I),I = 1, NNP )
1000 FORMAT(10F8.3)
C
READ(10,1100) ( ( NDE(I,J), I = 1, NEL ), J = 1, 4 )
1100 FORMAT(20I4)
C
C ----- At the Last LOP ---
C
READ(10,1000) ( VRUV( I), VARP( I), I = 1, LOP )
READ(10,1000) ( VARV( I), I = 1, LOP )
READ(10,1000) ( UN ( I), VN ( I), I = 1, NNP )
READ(10,1000) ( PSI( I), VOR( I), I = 1, NNP )
C
C ----- To set NDE numbers for a comparison with Giha's Data ---
C
NY = 0
C
DO 100 I = 1, NNP
C
IF ( DABS( X( I)-0.5D0 ).LE.0.01D0 ) THEN
NY = NY + 1 ! NY < = 200
KSY( NY) = I
END IF
C
100 CONTINUE
C
C ------------
C
NX = 0
C
DO 200 I = 1, NNP
C
IF ( DABS( Y( I)-0.5D0 ).LE.0.01D0 ) THEN
NX = NX + 1 ! NX < = 200
KSX( NX) = I
END IF
C
200 CONTINUE
C
WRITE(6,2100) NY, NX
2100 FORMAT(/,' *** NY =', I3, ' NX =',I3)
C
WRITE(6,*) ' '
WRITE(6,*) ' * KPLT = 1 ( Meshes / Str.Fn./Vel.Dis. ) or 2 ( Compa
1rison with Ref. )'
READ(5,*) KPLT
*** WRITE(6,*) ' '
*** WRITE(6,*) ' * KPLT =', KPLT
C
CALL VGA@
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLVR ( RAX, NPL, KKK )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX0, IY0
C
C ----- Attention NPL Data .LE. 500 ---
C
DIMENSION RAX( NPL)
C
INTEGER*2 IVAR( 500), IXX( 500)
C
CHARACTER*80 STR3
CHARACTER*6 STR1( 3)
C
DATA STR3 /'* F2D-FER (SV) ( Max.Var_of Str. funct., u-Velocity &
1 Vorticity(%) )'/
DATA STR1/'0 ','NPL/2',' NPL'/
C
C ----- Attention --- Max. Value ---
C
IF ( KKK.GE.2) GO TO 1234
C
ICR = 15
C
CALL SET_TEXT_ATTRIBUTE@( 3, 1., 0., 0. )
C
CALL DRAW_TEXT@( STR3, 70, 15, ICR )
C
ICR = 15
C
CALL DRAW_TEXT@( STR1( 1), 47, 450, ICR )
CALL DRAW_TEXT@( STR1( 2), 305, 450, ICR )
CALL DRAW_TEXT@( STR1( 3), 585, 450, ICR )
IF ( KKK.NE.1) GO TO 133
C
CALL DRAW_TEXT@( '1000', 15, 35, ICR )
CALL DRAW_TEXT@( ' 100', 15, 115, ICR )
CALL DRAW_TEXT@( ' 10', 15, 195, ICR )
C
CALL DRAW_TEXT@( ' 1', 15, 275, ICR )
CALL DRAW_TEXT@( ' 0.1', 15, 355, ICR )
CALL DRAW_TEXT@( '0.01', 15, 435, ICR )
C
CALL DRAW_TEXT@( '(%)', 15, 240, ICR )
C
133 CONTINUE
IX0 = 50
IY0 = 40
CALL RECTANGLE@( IX0, IY0, IX0+580, IY0+400, 15 )
C
1234 CONTINUE
IX0 = 50
IY0 = 40
C
RIWD = 580.D0/ DFLOAT( NPL)
C ================================
C
IF ( KKK.EQ.1) ICR = 15
IF ( KKK.EQ.2) ICR = 12
IF ( KKK.EQ.3) ICR = 9
IF ( KKK.EQ.4) ICR = 11
C
*** WRITE(6,*) ' * V_RAX=', RAX( NPL)
C
C ------------------------------------------------
DO 100 I = 1, NPL
C
IF ( RAX(I).LE.0.D0 ) GO TO 100
RAX( I) = DLOG10( RAX( I))
C
100 CONTINUE
C ------------------------------------------------
DO 200 I = 1, NPL
C
IXX( I) = IX0 + IDNINT( DFLOAT( I)* RIWD )
C
IVAR( I) = 280 - IDNINT( 80.D0* RAX( I) )
C
200 CONTINUE
C ------------------------------------------------
C
CALL POLYLINE@( IXX, IVAR, NPL, ICR )
C
C ------------------------------------------------
C
*** WRITE(6,*) ' * RAX=', RAX( NPL), ' IVAR=', IVAR( NPL)
C
RETURN
END
C *********************************************************************
C
SUBROUTINE CMPT ( UN, VN, X, Y, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IX1, IY1, IX2, IY2
C
DIMENSION UN( NNP), VN( NNP), X( NNP), Y( NNP)
C
COMMON /TIMC/ DT, LPM, LP1, LOP, ALP
C
NPL = 0
C ---------------
C
IX1 = 0
IY1 = 0
IX2 = 639
IY2 = 479
C
ICR = 7
CALL CLEAR_SCREEN_AREA@( IX1, IY1, IX2, IY2, ICR )
C
KEY = 0
CALL REFCMP ( X, Y, UN, VN, NNP, KEY, NPL )
KEY = KEY + 1
C
RETURN
END
C **********************************************************************
C
SUBROUTINE REFCMP ( X, Y, U, V, NNP, KEY, NPL )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
DIMENSION X( NNP), Y( NNP), U( NNP), V( NNP)
C
DIMENSION GU( 17), GV( 17), GUY( 17), GVX( 17)
C
COMMON /CMPD/ KSY( 200), NY, KSX( 200), NX
C
INTEGER*2 IX1, IY1, IX2, IY2, ICR, IICR
C
CHARACTER*80 STR0
C
DATA STR0 /'** Results of F2D-FER (SV) ( Velocity Distribution co
1mpared with Ghia''s Data )'/
C
IF ( KEY.EQ.0) ICR = 11
IF ( KEY.EQ.1) ICR = 13
IF ( KEY.EQ.2) ICR = 14
C
IICR = 16
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
DO 100 I = 1, 17
C
READ(8,1000) GU( I), GUY( I), GV( I), GVX( I)
1000 FORMAT(4F10.5)
C
100 CONTINUE
REWIND ( 8)
C
C ----- Attention : For Velocity direction change ---
C
DO 200 I = 1, 17
C
GU ( I) = - GU( I)
GVX( I) = 1.D0 - GVX( I)
C
200 CONTINUE
C *****************************
C
X1 = 345
Y1 = 50
X0 = 20
Y0 = 50
C
SAL = 260.D0
X2 = SAL* 0.5D0
Y2 = SAL* 0.5D0
C
C ----- U by Ghia ---
C
IX1 = X0
IY1 = Y0
IX2 = X0 + SAL
IY2 = Y0 + SAL
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X0 + X2
IY1 = Y0 + SAL
IX2 = X0 + X2
IY2 = Y0
C
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
C ----------------------------------------------------
DO 300 I = 1, 16
C
IC = I + 1
C
IX1 = X0 + X2 + GU ( I)* X2
IY1 = Y0 + SAL - GUY( I)* SAL
IX2 = X0 + X2 + GU (IC)* X2
IY2 = Y0 + SAL - GUY(IC)* SAL
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, IICR )
C
300 CONTINUE
C ----------------------------------------------------
C
IX1 = 40
IY1 = 80 + NPL*20
C
IX2 = 80
IY2 = IY1
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, IICR )
CALL DRAW_TEXT@( ' Ref.', IX2+5, IY1-7, IICR )
C
NPL = NPL + 1
C -------------------
C
C ----- Computed U -------------------------------
C
DO 400 J = 1, NY
C
NN = KSY( J)
IX1 = X0 + X2 + U( NN)* X2
IY1 = Y0 + SAL - Y( NN)* SAL
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 + SAL
IY2 = Y1 + SAL
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
IX1 = X1
IY1 = Y1 + Y2
IX2 = X1 + SAL
IY2 = Y1 + Y2
CALL RECTANGLE@( IX1, IY1, IX2, IY2, IICR )
C
C ------------------------------------------------
DO 500 I = 1, 16
C
IC = I + 1
IX1 = X1 + ( 1 - GVX( I))* SAL
IY1 = Y1 + Y2 + GV ( I) * Y2
IX2 = X1 + ( 1 - GVX( IC))* SAL
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 J = 1, NX
C
NN = KSX( J)
IX1 = X1 + X( NN)* SAL
IY1 = Y1 + Y2 - V( NN)* Y2
IX2 = 2
IY2 = 2
C
CALL ELLIPSE@( IX1, IY1, IX2, IY2, ICR )
C
600 CONTINUE
C --------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PMSH ( X, Y )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, ICR2, XIT, YIT, IWD, IWE, JWE
C
DIMENSION X( 121), Y( 121)
C
ICR = 15
ICR2 = 15
C
IWD = 160
XIT = 100
YIT = 70
C
C ------------------------------------------------
DO 100 I = 2, 10
C
JWE = NINT( Y( 1 + 12* ( I-1))* DFLOAT( IWD) )
CALL DRAW_LINE@( XIT, YIT+JWE, XIT+IWD, YIT+JWE, ICR2 )
C
100 CONTINUE
C ------------------------------------------------
DO 200 J = 1, 9
C
IWE = NINT( X( J+1)* DFLOAT( IWD) )
CALL DRAW_LINE@( XIT+IWE, YIT, XIT+IWE, YIT+IWD, ICR2 )
C
200 CONTINUE
C ----------------------------------------------------------
C
CALL RECTANGLE@( XIT, YIT, XIT+IWD, YIT+IWD, ICR )
C
C ----------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PCNT ( X, Y, NDE, TEP, NNP, NEL, KCR )
C
C ----- To Plot Contour Lines for Rectangular Meshes ---
C
C ( Linear Interpolation )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, NDE, X0, Y0, IWD, ICR2
COMMON /ABC / XINT, YINT, IWD
C
DIMENSION X ( NNP), Y ( NNP), NDE( NEL,4), TEP( NNP),
1 PX( 4), PY( 4), HVN( 21)
C
IWD = 160
ICR = 15
X0 = 280
Y0 = 70
C
CALL RECTANGLE@( X0, Y0, X0+IWD, Y0+IWD, ICR )
C
SCX = 160.D0
SCY = 160.D0
NM = 10
C
C ---------------------------------------------------------------
TMIN = 100.D0
TMAX = -100.D0
C ----------------------------------------------------------
DO 100 I = 1, NNP
C
IF ( TEP( I).GE.TMAX ) TMAX = TEP( I)
IF ( TEP( I).LE.TMIN ) TMIN = TEP( I)
C
100 CONTINUE
C ----------------------------------------------------------
WRITE(6,2000) TMAX, TMIN
2000 FORMAT(/,' * Cmax.=',F7.3,' * Cmin.=',F7.3)
C
C ----------------------------------------------------------
RTW = ( TMAX - TMIN ) / DFLOAT( NM)
DO 200 I = 1, NM + 1
C
HVN( I) = TMIN + RTW* DFLOAT( I-1)
IF ( HVN( I).LT.0.D0 ) NMZ = I
C
200 CONTINUE
C ----------------------------------------------------------
C
C ----- To Plot Contour Lines ---
C
IF ( KCR.EQ.0) ICR2 = 11
IF ( KCR.EQ.1) ICR2 = 12
IF ( KCR.EQ.2) ICR2 = 14
C
PXV = 0.D0
PYV = 0.D0
C
C --------------------------------------------------------------------
DO 300 L = 1, NEL
C
C ----------------------------------------------------------
DO 320 IH = 1, NM
C
HV = HVN( IH)
IL = 0
C
C ------------------------------------------------
DO 330 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
IF ( DABS( TEP( I1)-TEP( I2) ).LE.1.D-4 ) THEN
C
IF ( DABS( TEP( I1)-HV ).LE.1.D-4 ) THEN
IL = IL + 1
PX( IL) = X( I1)* SCX
PY( IL) = Y( I1)* SCY
END IF
C
ELSE IF ( ( ( TEP( I1).LE.HV ) .AND. ( HV.LT.TEP( I2) ) )
1 .OR. ( ( TEP( I1).GE.HV ) .AND. ( HV.GT.TEP( I2) ) ) ) THEN
C
IL = IL + 1
RT = ( HV - TEP( I1) ) / ( TEP( I2) - TEP( I1) )
C
PX( IL) = ( X( I1) + ( X( I2) - X( I1) )* RT )* SCX
PY( IL) = ( Y( I1) + ( Y( I2) - Y( I1) )* RT )* SCY
END IF
C
330 CONTINUE
C ------------------------------------------------
C
IF ( IL.GE.2) THEN
C ------------------------------------------------
DO 340 J = 1, IL
C
JM = J
C ----- Attention : D2 = ( PXV-PX( J) )**2 + ( PYV-PY( J) )**2 ---

D2 = ( PXV - PX( J) )**2 + ( PYV - PY( J) )**2
C
C --------------------------------------
DO 360 JJ = J+1, IL
C
C --- Attention ---
C
C ----- D2V = ( PXV - PX( JJ) )**2 + ( PYV - PY( JJ) )**2 ---
C
D2V = ( PXV - PX( JJ) )** 2 +( PYV - PY( JJ) )** 2
IF ( D2V.LT.D2) THEN
JM = JJ
D2 = D2V
END IF
C
360 CONTINUE
C --------------------------------------
PXV = PX( JM)
PYV = PY( JM)
C
PX( JM) = PX( J)
PY( JM) = PY( J)
C
PX( J) = PXV
PY( J) = PYV
C
340 CONTINUE
C ------------------------------------------------
PXV = PX( IL)
PYV = PY( IL)
C
END IF
C ------------------------------------------------
C
IF ( IL.GT.0) THEN
C ------------------------------------------------
DO 380 J = 2, IL
C
CALL DRAW_LINE@
C
1 ( X0+IDNINT( PY( J-1)), Y0+IDNINT( PX( J-1)),
2 X0+IDNINT( PY( J)), Y0+IDNINT( PX( J)), ICR2 )
C
380 CONTINUE
C -----------------------------------------------
END IF
C
320 CONTINUE
C ----------------------------------------------------------
C
300 CONTINUE
C --------------------------------------------------------------------
C
RETURN
END
C **********************************************************************
C
SUBROUTINE TTLE
C
IMPLICIT REAL*8 (A-H,O-Z)
C
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
WRITE(6,*) ' '
WRITE(6,*) ' Fluid Analysis Software by FEM '
WRITE(6,*) ' '
WRITE(6,*) ' F I G - F 2 D - F E R ( SV ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Version : 3.5 ) '
WRITE(6,*) ' '
WRITE(6,*) ' ( Modified Galerkin Method - Rectangular Linear ) '
WRITE(6,*) ' '
WRITE(6,*) ' Copyright 2011 : Yasuhiro MATSUDA '
WRITE(6,*) ' '
WRITE(6,*) '****************************************************'
C
RETURN
END
C **********************************************************************
C
SUBROUTINE PLTN ( X, Y, UN, VN, NNP )
C
IMPLICIT REAL*8 (A-H,O-Z)
C
INTEGER*2 ICR, IWD, IX0, IY0, IX1, IY1, IX2, IY2, GX, GY
C
DIMENSION X( NNP), Y( NNP), UN( NNP), VN( NNP)
C
COMMON /TIMC/ DT, LPM, LP1, LOP, ALP
C
SCX = 240.D0
SCY = 240.D0
SCL = 20.D0
C
ICR = 15
IX0 = 460
IY0 = 70
IWD = 160
c
TL = DT* DFLOAT( LOP)
C
IX1 = IX0 + IWD
IY1 = IY0 + IWD
C
CALL RECTANGLE@( IX0, IY0, IX1, IY1, ICR )
C
ICR = 11
C
C ------------------------------------------------
DO 100 I = 1, NNP
C
GX = IX0 + X( I)* IWD
GY = IY0 + IWD - Y( I)* IWD
C
IF ( GX.EQ.IX0 ) GO TO 100
IF ( GY.EQ.IY0 ) GO TO 100
IF ( GX.EQ.IX0+IWD ) GO TO 100
IF ( GY.EQ.IY0+IWD ) GO TO 100
IX1 = GX
IY1 = GY
IX2 = GX + UN( I)* SCL
IY2 = GY - VN( I)* SCL
C
CALL DRAW_LINE@( IX1, IY1, IX2, IY2, ICR )
C
100 CONTINUE
C ------------------------------------------------
C
RETURN
END
C **********************************************************************