PROGRAM POSTTAU
C=======================================================================
C FINITE ELEMENT POST-TREATMENT FOR TWO-DIMENSIONAL SOLID SOLUTIONS
C ELEMENT TYPE: FOUR-NODED ISOPARAMETRIC ELEMENT
C ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C=======================================================================
INCLUDE 'PARAM.DAT'
CC PARAMETER ( MXE=20000,MXN=23000,MXB=5000, ND=4,NP=ND, NP0=1 )
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
C ARRAYS
CHARACTER*12 INPFILE
DIMENSION X(ND),Y(ND),BX(2,ND), E1(ND),E2(ND), TAUND(3,ND)
DIMENSION BP(2,ND,ND)
DIMENSION NODEX(MXE,ND)
DIMENSION IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB)
DIMENSION XCOORD(MXN),YCOORD(MXN),U(MXN),V(MXN),
* DIVV(MXN), ICONNECT(MXN), GLBTAU(3,MXN)
C=======================================================================
DATA E1 / -1.D0, +1.D0, +1.D0, -1.D0 /
DATA E2 / -1.D0, -1.D0, +1.D0, +1.D0 /
DATA INPFILE /'STATIC04.DAT'/
C=======================================================================
CALL DERIVA ( ND, E1, E2, X, Y, BP )
C=======================================================================
WRITE(*,*)' STATIC STRESSES COMPUTATION PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( INPFILE,MXE,MXN,MXB,ND,NE,NNODE,NBFX,NBFY,YOUNG,
* POISSON,NODEX,XCOORD,YCOORD, IBNDFX,IBNDFY,BVX,BVY )
WRITE (*,*)' NE=', NE, ' NNODE=',NNODE
C=======================================================================
CALL INPUTUV ( MXN,NNODE,U,V )
C=======================================================================
WRITE(*,*)' YOUNG MODULUS =', YOUNG
WRITE(*,*)' POISSON RATIO =', POISSON
VISCO = YOUNG / (2.*( 1. + POISSON ))
FLMDA = YOUNG*POISSON / ( (1.+ POISSON)*(1.- 2.* POISSON) )
WRITE(*,*)' VISCO = ', VISCO
WRITE(*,*)' FLMDA = ', FLMDA
C=======================================================================
CALL RELATION ( MXE, MXN, ND, NE, NNODE, NODEX, ICONNECT )
WRITE (*,*) 'AFTER RELATION'
CALL DIVERG ( MXE,MXN,ND,BP,FLMDA,NE,NNODE,XCOORD,YCOORD,
* NODEX,BX,U,V, DIVV, ICONNECT )
WRITE (*,*) 'AFTER DIVERG'
CALL STRESS ( MXE,MXN,ND,BX,VISCO,NE,NNODE,XCOORD,YCOORD,NODEX,
* U, V, DIVV,TAUND,GLBTAU, ICONNECT,BP )
WRITE (*,*) 'AFTER STRESS'
CALL FILEMK ( MXN, NNODE, DIVV, GLBTAU )
STOP' NORMAL TERMINATION'
END
C
C
SUBROUTINE FILEMK ( MXN, NNODE, DIVV, GLBTAU )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION DIVV(MXN), GLBTAU(3,MXN)
CHARACTER OUTFILE*12, EXFILE*3
LOGICAL YES
C========> DIVV(I)=DIVERGENCE,
C========> GLBTAU(1,I)=TAUXX = MU*(2DU/DX)+LAMBDA*DIVV
C========> GLBTAU(1,I)=TAUYY = MU*(2DV/DY)+LAMBDA*DIVV
C========> GLBTAU(1,I)=TAUXY = MU*(DU/DY+DU/DX)
OUTFILE = 'STRESS04.DAT'
EXFILE = 'NEW'
INQUIRE ( FILE=OUTFILE, EXIST=YES )
IF ( YES ) EXFILE='OLD'
OPEN ( 1, FILE=OUTFILE, STATUS=EXFILE )
WRITE (1,*) ' ----- STRESSES -----'
WRITE (1,110)
WRITE (1,120)
TAUXXMAX = GLBTAU(1,1)
TAUXXMIN = GLBTAU(1,1)
TAUYYMAX = GLBTAU(2,1)
TAUYYMIN = GLBTAU(2,1)
TAUXYMAX = GLBTAU(3,1)
TAUXYMIN = GLBTAU(3,1)
DO I = 1 , NNODE
TAUXXMAX = DMAX1 ( GLBTAU(1,I), TAUXXMAX )
TAUXXMIN = DMIN1 ( GLBTAU(1,I), TAUXXMIN )
TAUYYMAX = DMAX1 ( GLBTAU(2,I), TAUYYMAX )
TAUYYMIN = DMIN1 ( GLBTAU(2,I), TAUYYMIN )
TAUXYMAX = DMAX1 ( GLBTAU(3,I), TAUXYMAX )
TAUXYMIN = DMIN1 ( GLBTAU(3,I), TAUXYMIN )
END DO
C
ROOT2 = 1.D0 / DSQRT(2.D0)
YMAX = 0.D0
YMIN = 0.D0
DO I = 1 , NNODE
TAUAVE = ( GLBTAU(1,I) + GLBTAU(2,I) ) / 2.D0
TAUDIFF= ( GLBTAU(1,I) - GLBTAU(2,I) ) / 2.D0
ADDI = DSQRT ( TAUDIFF*TAUDIFF + GLBTAU(3,I)*GLBTAU(3,I) )
TAU1 = TAUAVE + ADDI
TAU2 = TAUAVE - ADDI
YMISES = DSQRT ( (TAU1-TAU2)**2 + TAU1**2 + TAU2**2 ) * ROOT2
WRITE (1,100) I, GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I), YMISES,
* TAU1, TAU2
YMAX = DMAX1 ( YMAX, YMISES )
YMIN = DMIN1 ( YMIN, YMISES )
IF ( I .EQ. 1 ) THEN
TAU1MAX = TAU1
TAU1MIN = TAU1
TAU2MAX = TAU2
TAU2MIN = TAU2
ELSE
TAU1MAX = DMAX1 ( TAU1MAX, TAU1 )
TAU1MIN = DMIN1 ( TAU1MIN, TAU1 )
TAU2MAX = DMAX1 ( TAU2MAX, TAU2 )
TAU2MIN = DMIN1 ( TAU2MIN, TAU2 )
END IF
END DO
WRITE (1,120)
WRITE (1,121) TAUXXMIN,TAUYYMIN,TAUXYMIN,
* YMIN,TAU1MIN,TAU2MIN
WRITE (1,122) TAUXXMAX, TAUYYMAX, TAUXYMAX,
* YMAX,TAU1MAX,TAU2MAX
CLOSE (1)
100 FORMAT ( 1X,I5,6G18.7 )
110 FORMAT ( 1X,'NODE#',' TAUXX',' TAUYY',
* ' TAUXY', ' MISES',
* ' TAU1 TAU2')
120 FORMAT ( 1X,5('-'), 6( 5X,9('-'),4X ) )
121 FORMAT ( 1X,' MIN ',6G18.7 )
122 FORMAT ( 1X,' MAX ',6G18.7 )
RETURN
END
C
C
SUBROUTINE DERIVA ( ND, E1, E2, F1, F2, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION BPP(2,ND,ND), F1(ND), F2(ND), E1(ND), E2(ND)
C------- COMPUTATION OF BP(J) = D N(J) / D ETA1
DO K = 1 , ND
X = E1(K) + 0.5
CALL ISOPARA ( ND , X , E2(K) , F1 )
X = E1(K) - 0.5
CALL ISOPARA ( ND , X , E2(K) , F2 )
DO I = 1 , ND
BPP(1,I,K) = F1(I) - F2(I)
END DO
C------- COMPUTATION OF BP(J) = D N(J) / D ETA2
Y = E2(K) + 0.5
CALL ISOPARA ( ND , E1(K) , Y , F1 )
Y = E2(K) - 0.5
CALL ISOPARA ( ND , E1(K) , Y , F2 )
DO I = 1 , ND
BPP(2,I,K) = F1(I) - F2(I)
END DO
END DO
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND , E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.25 * (1.-E1)*(1.-E2)
F(2) = 0.25 * (1.+E1)*(1.-E2)
F(3) = 0.25 * (1.+E1)*(1.+E2)
F(4) = 0.25 * (1.-E1)*(1.+E2)
RETURN
END
C
C
SUBROUTINE RELATION ( MXE, MXN, ND, NE, NNODE, NODEX, ICONNECT )
DIMENSION NODEX(MXE,ND), ICONNECT(MXN)
DO I = 1 , NNODE
ICONNECT(I) = 0
END DO
DO IEL = 1 , NE
DO I = 1 , ND
ICONNECT(NODEX(IEL,I)) = ICONNECT(NODEX(IEL,I)) + 1
END DO
END DO
RETURN
END
C
C
SUBROUTINE STRESS ( MXE,MXN,ND,BX,VISCO,NE,NNODE,XCOORD,YCOORD,
* NODEX,U,V,DIVV,TAUND,GLBTAU, ICONNECT,BP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),BX(2,ND),
* U(MXN), V(MXN),DIVV(MXN),TAUND(3,ND), GLBTAU(3,MXN),
* ICONNECT(MXN), BP(2,ND,ND)
C-------- COMPUTATION OF STRESS ======== RETURN VALUES ====> GLBTAU
C-------- DERIVATIVES ARE EVALUATED AT NODAL POINTS(E1,E2).
DO K = 1 , 3
DO I = 1 , NNODE
GLBTAU(K,I) = 0.
END DO
END DO
DO 500 IEL = 1 , NE
DO 400 K = 1 , ND
YAC11 = 0.
YAC12 = 0.
YAC21 = 0.
YAC22 = 0.
DO I = 1 , ND
NODE = NODEX(IEL,I)
YAC11 = YAC11 + BP(1,I,K) * XCOORD(NODE)
YAC12 = YAC12 + BP(1,I,K) * YCOORD(NODE)
YAC21 = YAC21 + BP(2,I,K) * XCOORD(NODE)
YAC22 = YAC22 + BP(2,I,K) * YCOORD(NODE)
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
Y11 = YAC22 / DETJ
Y12 = -YAC12 / DETJ
Y21 = -YAC21 / DETJ
Y22 = YAC11 / DETJ
DO J = 1 , ND
BX(1,J) = Y11 * BP(1,J,K) + Y12 * BP(2,J,K)
BX(2,J) = Y21 * BP(1,J,K) + Y22 * BP(2,J,K)
END DO
DUDX = 0.
DUDY = 0.
DVDX = 0.
DVDY = 0.
DO I = 1 , ND
NODE = NODEX(IEL,I)
DUDX = DUDX + BX(1,I)*U(NODE)
DVDX = DVDX + BX(1,I)*V(NODE)
DVDY = DVDY + BX(2,I)*V(NODE)
DUDY = DUDY + BX(2,I)*U(NODE)
END DO
TAUND(1,K)= VISCO*(DUDX+DUDX)
TAUND(2,K)= VISCO*(DVDY+DVDY)
TAUND(3,K)= VISCO*(DUDY+DVDX)
400 CONTINUE
DO I = 1 , ND
NODE = NODEX(IEL,I)
GLBTAU(1,NODE) = GLBTAU(1,NODE) + TAUND(1,I)
GLBTAU(2,NODE) = GLBTAU(2,NODE) + TAUND(2,I)
GLBTAU(3,NODE) = GLBTAU(3,NODE) + TAUND(3,I)
END DO
500 CONTINUE
DO I = 1 , NNODE
GLBTAU(1,I) = GLBTAU(1,I) / ICONNECT(I) + DIVV(I)
GLBTAU(2,I) = GLBTAU(2,I) / ICONNECT(I) + DIVV(I)
GLBTAU(3,I) = GLBTAU(3,I) / ICONNECT(I)
END DO
RETURN
END
C
C
SUBROUTINE DIVERG ( MXE,MXN,ND,BP,FLMDA,NE,NNODE,XCOORD,
* YCOORD,NODEX,BX,U,V, DIVV, ICONNECT )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),ICONNECT(MXN),
* BX(2,ND),BP(2,ND,ND),DIVV(MXN), U(MXN), V(MXN)
C-------- COMPUTATION OF LAMBDA*DIVERGENCE OF STRAIN
C-------- RETURN VALUES DIVV
DO I = 1 , NNODE
DIVV(I) = 0.
END DO
DO 500 IEL = 1 , NE
DO 400 K = 1 , ND
YAC11 = 0.
YAC12 = 0.
YAC21 = 0.
YAC22 = 0.
DO I = 1 , ND
NODE = NODEX(IEL,I)
YAC11 = YAC11 + BP(1,I,K) * XCOORD(NODE)
YAC12 = YAC12 + BP(1,I,K) * YCOORD(NODE)
YAC21 = YAC21 + BP(2,I,K) * XCOORD(NODE)
YAC22 = YAC22 + BP(2,I,K) * YCOORD(NODE)
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
Y11 = YAC22 / DETJ
Y12 = -YAC12 / DETJ
Y21 = -YAC21 / DETJ
Y22 = YAC11 / DETJ
DO J = 1 , ND
BX(1,J) = Y11 * BP(1,J,K) + Y12 * BP(2,J,K)
BX(2,J) = Y21 * BP(1,J,K) + Y22 * BP(2,J,K)
END DO
C-------- COMPUTATION OF DU/DX AND DV/DY
DUDX = 0.
DVDY = 0.
DO I = 1 , ND
DUDX = DUDX + BX(1,I)*U(NODEX(IEL,I))
DVDY = DVDY + BX(2,I)*V(NODEX(IEL,I))
END DO
400 CONTINUE
DO I = 1 , ND
NODE = NODEX(IEL,I)
DIVV(NODE) = DIVV(NODE) + FLMDA * (DUDX+DVDY) / NP0
END DO
500 CONTINUE
DO I = 1 , NNODE
DIVV(I) = DIVV(I) / ICONNECT(I)
END DO
RETURN
END
C
C
SUBROUTINE INPUT ( INPFILE,MXE,MXN,MXB,ND,NE,NNODE,NBFX,NBFY,
* YOUNG,POISSON, NODEX,XCOORD,YCOORD,IBNDFX,IBNDFY,BVX,BVY )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB)
CHARACTER*12 INPFILE
LOGICAL YES
C========> FILENAME INPFILE SEE MAIN PROGRAM
INQUIRE ( FILE=INPFILE, EXIST=YES )
IF ( YES ) THEN
OPEN ( 1, FILE=INPFILE, STATUS='OLD' )
ELSE
WRITE (*,*)' INPUT FILE DOES NOT EXIST'
STOP
ENDIF
C========> PARAMETERS
READ (1,*) YOUNG, POISSON
C========> ELEMENTS
READ (1,*) NE
IF ( NE .GT. MXE) STOP'ERROR #1'
IF ( NE .LE. 0 ) STOP'NE =0'
DO I = 1 , NE
READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
END DO
C========> FILENAME COORDINATES OF NODAL POINTS
READ (1,*) NNODE
IF ( NNODE .GT. MXN) STOP'ERROR #2'
IF ( NNODE .LE. 0 ) STOP'NNODE =0'
DO I = 1 , NNODE
READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
END DO
C========> DIRICHLET TYPE BOUNDARY CONDITIONS
READ (1,*) NBFX
IF ( NBFX .GT. MXB) STOP'NBFX>MXB'
IF ( NBFX .LT. 1 ) STOP'NBFX<1'
DO I = 1 , NBFX
READ (1,*) IBNDFX(I) , BVX(I)
END DO
READ (1,*) NBFY
IF ( NBFY .GT. MXB) STOP'NBFY .GT. MXB'
IF ( NBFY .LT. 1 ) STOP'NBFY<1'
DO I = 1 , NBFY
READ (1,*) IBNDFY(I) , BVY(I)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INPUTUV ( MXN,NNODE,U,V )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION U(MXN) , V(MXN)
CHARACTER INPFILE*12
LOGICAL YES
C========> FILENAME DISPLACE.MNT
INPFILE = 'DISPLACE.MNT'
INQUIRE ( FILE=INPFILE, EXIST=YES )
IF ( .NOT. YES ) STOP 'NO DISPLACE.MNT FILE EXIST'
OPEN (1,FILE=INPFILE,STATUS='OLD',FORM='UNFORMATTED')
READ (1) ( U(I) , I = 1 , NNODE )
READ (1) ( V(I) , I = 1 , NNODE )
CLOSE (1)
WRITE (*,*) ' DISPLACE.MNT READING DONE'
RETURN
END