PROGRAM SETAXISYMM4_DISC
C=======================================================================
C ********** 4-NODED ELEMENT **********
C------------------- ELEMENT NUMBERING: HORIZONTAL SCAN ----------------
C DATA GENERATING PROGRAM FOR STATIC8QFXAXISYMM-NEW.FOR
C DOMAIN SIZE: TLR BY TLZ
C BOUNDARY CONDITIONS
C NOV. 28, 2012
C EIJI FUKUMORI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4, INTEPT=2, MXE=31300, MXN=34410, MXB=21000 )
PARAMETER ( TLR0=0.D0, TLR1=5.D0, TLZ=0.05D0 )
PARAMETER ( YOUNG=210.D9, POISSON=0.3D0 )
PARAMETER ( NER=400,NEZ=40 )
PARAMETER ( P0=1.D0 )
C=======================================================================
DIMENSION NODEX(MXE,ND), RCOORD(MXN), ZCOORD(MXN),
* IBND1R(MXB),IBND1Z(MXB),BV1R(MXB),BV1Z(MXB),
* IBND2R(MXB),IBND2Z(MXB),BV2R(MXB),BV2Z(MXB)
DIMENSION NEUTRAL(MXN)
CHARACTER PROJECT*12,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'STATIC4.DAT' /
C=======================================================================
C NER : NEMBER OF ELEMENTS IN RADIAL DIRECTION
C NEZ: NEMBER OF ELEMENTS IN circumferential DIRECTION
C P0, P1, P2: APPLIED HYDRAULIC PRESSURE
C=======================================================================
THETA0 = DATAN ( 1.D0 )
THETA1 = 1.5D0*DATAN ( 1.D0 )
C=======================================================================
DZ = TLZ/NEZ
DR = (TLR1-TLR0) / NER
NDR = NER + 1
NDZ = NEZ + 1
C=======================================================================
WRITE (*,*)' YOUNG MODULUS = ',YOUNG
WRITE (*,*)' POISSON RATIO = ',POISSON
C=======================================================================
C ELEMENT CREATION
NE = 0
DO I = 1 , NER
DO J = 1 , NEZ
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = (I-1)*NDZ + J
NODEX(NE,2) = NODEX(NE,1) + NDZ
NODEX(NE,3) = NODEX(NE,2) + 1
NODEX(NE,4) = NODEX(NE,1) + 1
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
NNODE = 0
DO I = 1 , NDR
DO J = 1 , NDZ
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
RCOORD(NNODE) = (I-1)*DR+TLR0
ZCOORD(NNODE) = (J-1)*DZ
END DO
END DO
C=======================================================================
PI = 4.D0* ATAN (1.D0)
C BOUNDARY CONDITIONS
C==== FIRST KIND
C--------- NAVIER EQUATIONS
NB1R = 0
NB1Z = 0
NB2R = 0
NB2Z = 0
C
C------- BOTTOM RIGHT END
NODE = NDZ*NER+1
NB1Z = NB1Z + 1
IBND1Z(NB1Z) = NODE
BV1Z(NB1Z) = 0.D0
C
C
C------- LEFT END
DO I = 1 , NDZ
NB1R = NB1R + 1
IBND1R(NB1R) = I
BV1R(NB1R) = 0.D0
END DO
C==== SECOND KIND
C-------RESET
DO I = 1 , NDR
BV2Z(I) = 0.D0
END DO
C--------- LOAD DISTRIBUTION FOR PARABOLIC ELEMENT UNDER UNIFORM LOAD
DO I = 1 , NER
NODE1 = I*NDZ
NODE2 = (I+1)*NDZ
R1 = RCOORD(NODE1)
R3 = RCOORD(NODE2)
R2 = (R1+R3)/2.D0
A1 = PI*(R2**2 - R1**2)
A2 = PI*(R3**2 - R2**2)
C--------Z
NB2Z = I
IBND2Z(NB2Z) = NODE1
BV2Z (NB2Z) = BV2Z(NB2Z) + A1*P0
NB2Z = I + 1
IBND2Z(NB2Z) = NODE2
BV2Z (NB2Z) = A2*P0
END DO
C=======================================================================
WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
WRITE (*,*) ' NUMBER OF DIRICHLET R (NB1R) = ',NB1R
WRITE (*,*) ' NUMBER OF DIRICHLET Z (NB1Z) = ',NB1Z
WRITE (*,*) ' NUMBER OF NEUMANN R (NB2R) = ',NB2R
WRITE (*,*) ' NUMBER OF NEUMANN Z (NB2Z) = ',NB2Z
C=======================================================================
C DATA FILE INQUIRY
EXFILE = 'NEW'
INQUIRE ( FILE = PROJECT, EXIST = YES )
IF ( YES ) EXFILE='OLD'
C=======================================================================
C MAKING DATA FILES
C---------- 'PROJECT'.JNK
IR = 1
C---------- PARAMETERS
OPEN ( IR, FILE=PROJECT, STATUS = EXFILE )
WRITE(1,*) YOUNG , POISSON
C---------- ELEMENTS
C---------- ELEMENTS
WRITE(1,*) NE
DO I = 1 , NE
WRITE (1,*) I, (NODEX(I,J), J = 1 , ND )
END DO
C---------- COORDINATES OF NONAL POINTS
WRITE(1,*) NNODE
DO I = 1 , NNODE
WRITE(1,*) I,RCOORD(I), ZCOORD(I)
END DO
C---------- DIRICHLET TYPE BOUNDARY CONDITIONS
WRITE(1,*) NB1R
IF ( NB1R .GT. 0 ) THEN
DO I = 1 , NB1R
WRITE (1,*) IBND1R(I), BV1R(I)
END DO
END IF
WRITE(1,*) NB1Z
IF ( NB1Z .GT. 0 ) THEN
DO I = 1 , NB1Z
WRITE (1,*) IBND1Z(I), BV1Z(I)
END DO
END IF
C---------- NUEMANN TYPE BOUNDARY CONDITIONS
WRITE(1,*) NB2R
IF ( NB2R .GT. 0 ) THEN
DO I = 1 , NB2R
WRITE(1,*) IBND2R(I), BV2R(I)
END DO
ENDIF
WRITE(1,*) NB2Z
IF ( NB2Z .GT. 0 ) THEN
DO I = 1 , NB2Z
WRITE(1,*) IBND2Z(I), BV2Z(I)
END DO
ENDIF
C---------- FINAL
CLOSE (1)
C---------- ELEMENT DRAWING
OPEN ( 1, FILE='ELEMENT4.DAT', STATUS = 'UNKNOWN')
DO I = 1, NE
DO J = 1, ND
WRITE (1,*) RCOORD(NODEX(I,J)), ZCOORD(NODEX(I,J))
END DO
WRITE (1,*) RCOORD(NODEX(I,1)), ZCOORD(NODEX(I,1))
WRITE (1,*)
END DO
CLOSE (1)
C=======================================================================
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
WRITE (1,*) ' PARAMETER ( ND=',ND,', INTEPT=',INTEPT,' )'
WRITE (1,*) ' PARAMETER ( MXE=',NE,', MXN=',NNODE,
* ', MXB=',MAX0(NBFX,NBFY,NFORCEX,NFORCEY), ', MXW=',NBW, ' )'
CLOSE (1)
C=======================================================================
OPEN ( 1, FILE='PRINTDISP.DAT', STATUS='UNKNOWN' )
WRITE(1,*) NDR
DO I = 1 , NDR
WRITE(1,*) I*NDZ
END DO
CLOSE (1)
C=======================================================================
CALL EXACT ( MXE,MXN,ND,P0,NNODE,NE,RCOORD,ZCOORD,
* YOUNG,POISSON )
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
DIMENSION NODEX(MXE,ND)
C------- RETURN VALUE: NBW
NBW = 0
DO I = 1 , NE
DO J = 1 , ND-1
DO K = J+1 , ND
NBW = MAX0 ( NBW , IABS(NODEX(I,J)-NODEX(I,K)) )
END DO
END DO
END DO
NBW = NBW + 1
WRITE (*,*) ' HALH BANDWIDTH =', NBW
RETURN
END
C
C
SUBROUTINE EXACT ( MXE,MXN,ND,P0,NNODE,NE,RCOORD,ZCOORD,
* YOUNG,POISSON )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), RCOORD(MXN), ZCOORD(MXN)
C-------- Q = P0 ---------
C = (ZCOORD(NNODE) - ZCOORD(1))/2.D0
A = RCOORD(NNODE) - RCOORD(1)
C1 = (2.D0+POISSON)/(8.D0*C**3)
C2 = -3.D0*(3.D0+POISSON)/(32.D0*C**3)
C3 = -(3.D0/8.D0)*(2.D0+POISSON)/(5.D0*C)
C4 = -C2*A**2
C------- EXACT AT SELECTED PNODAL POINTS BY PRINTDISP.DAT ---
OPEN(3,FILE='EXACT.OUT',STATUS='UNKNOWN')
WRITE(3,*) 'NODE R Z TAURR TAUZZ TAURZ'
OPEN ( 1, FILE='PRINTDISP.DAT', STATUS='UNKNOWN' )
PREVIOUS = ZCOORD(NNODE)
READ(1,*) N
DO K = 1 , N
READ(1,*) I
Z = -(ZCOORD(I)-C)
R = RCOORD(I)
IF ( PREVIOUS .NE. ZCOORD(I) ) WRITE(3,*)
WRITE(3,*) I,R,ZCOORD(I),-TAURR(P0,C1,C2,C3,C4,R,Z),
* -TAUZZ(P0,C,Z),TAURZ(P0,R,C,Z)
PREVIOUS = ZCOORD(I)
END DO
CLOSE (3)
CLOSE (1)
C-------------------- EXACT AT ALL NODAL POINT --------------
OPEN(3,FILE='EXACTALLNODALPOINTS.OUT',STATUS='UNKNOWN')
WRITE(3,*) 'NODE R Z TAURR TAUZZ TAURZ'
DO I = 1 , NNODE
Z = -(ZCOORD(I)-C)
R = RCOORD(I)
WRITE(3,*) I,R,ZCOORD(I),-TAURR(P0,C1,C2,C3,C4,R,Z),
* -TAUZZ(P0,C,Z),TAURZ(P0,R,C,Z)
PREVIOUS = ZCOORD(I)
END DO
CLOSE (3)
RETURN
END
C
C
FUNCTION TAURR (Q,C1,C2,C3,C4,R,Z)
IMPLICIT REAL*8 ( A-H , O-Z )
TAURR = Q*(C1*Z**3 + C2*R**2*Z + C3*Z + C4*Z)
RETURN
END
C
C
FUNCTION TAUZZ (Q,C,Z)
IMPLICIT REAL*8 ( A-H , O-Z )
TAUZZ = Q*(-Z**3/(4.D0*C**3)+(3.D0*Z/(4.D0*C))-0.5D0)
RETURN
END
C
C
FUNCTION TAURZ (Q,R,C,Z)
IMPLICIT REAL*8 ( A-H , O-Z )
TAURZ = -3.D0*Q*R/(8.D0*C**3)*(C**2-Z**2)
RETURN
END