PROGRAM POST3DTAUXX
C=======================================================================
C FINITE ELEMENT POST-TREATMENT FOR 3-DIMENSIONAL SOLID SOLUTIONS
C APPLIED TO TORSION PROBLEM
C ELEMENT TYPE: 8-NODED ISOPARAMETRIC HEXAGONAL ELEMENT
C PROGRAM NAME OF 3-DIMENSIONAL SOLID ANALYSIS:3DSTATIC8QFXCOMBINE.FOR
C ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C 23JUNE 2010, 09/NOV./2024
C U(1,I) = U(I), U(2,I) = V(I), U(3,I) = W(I)
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
INCLUDE 'PARAM.DAT'
CCCCC PARAMETER ( MXE=40000,MXN=53000,MXB=5000, ND=8 )
C=======================================================================
C ARRAYS
CHARACTER*14 INPFILE
DIMENSION F1(ND),F2(ND),BX(3,ND),E(3,ND),
* TAUND(7,ND),BPP(3,ND,ND)
DIMENSION NODEX(MXE,ND)
DIMENSION XCOORD(3,MXN),U(3,MXN),
* ICONNECT(MXN), GLBTAU(7,MXN),NEUTRAL(MXN)
DIMENSION NODECTOP(MXN)
C=======================================================================
WRITE (*,*) ND
INPFILE = 'NONE'
IF ( ND .EQ. 8 ) THEN
INPFILE = 'STATIC3D08.DAT'
END IF
IF ( ND .EQ. 27 ) THEN
INPFILE = 'STATIC3D27.DAT'
END IF
IF ( INPFILE .EQ. 'NONE' ) STOP 'NO INPUT FILE AT MAIN'
C=======================================================================
CALL DERIVEPT ( ND, E )
C=======================================================================
CALL DERIVA ( ND, E, F1, F2, BPP )
C=======================================================================
WRITE(*,*)' STATIC STRESSES COMPUTATION PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( INPFILE,MXE,MXN,MXB,ND,NE,NNODE,
* YOUNG,POISSON, NODEX,XCOORD,NNEUT, NEUTRAL, NCTOP,NODECTOP )
WRITE (*,*)' NE=', NE, ' NNODE=',NNODE
C=======================================================================
CALL INPUTUV ( MXN,NNODE,U )
C=======================================================================
WRITE(*,*)' YOUNG MODULUS =', YOUNG
WRITE(*,*)' POISSON RATIO =', POISSON
VISCO = YOUNG / (2.*( 1. + POISSON ))
FLMDA = YOUNG*POISSON / ( (1.+ POISSON)*(1.- 2.* POISSON) )
WRITE(*,*)' SHEAR MODULUS = ', VISCO
WRITE(*,*)' LAMBDA = ', FLMDA
C=======================================================================
WRITE (*,*) 'RELATION'
CALL RELATION ( MXE, MXN, ND, NE, NNODE, NODEX, ICONNECT )
WRITE (*,*) 'STRESS '
CALL STRESS ( MXE,MXN,ND,BX,VISCO,NE,NNODE,XCOORD,
* NODEX,U,TAUND,GLBTAU, ICONNECT,BPP,FLMDA )
WRITE (*,*) 'FILEMK '
CALL FILEMK ( ND,MXN, NNODE, GLBTAU, U,NNEUT,NEUTRAL,XCOORD )
C=======================================================================
CALL ANGLE ( MXN, NNEUT, NEUTRAL, NCTOP,NODECTOP,XCOORD,U )
STOP' NORMAL TERMINATION'
END
C
C
SUBROUTINE FILEMK ( ND,MXN,NNODE,GLBTAU,U,NNEUT,NEUTRAL,XCOORD )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION GLBTAU(7,MXN),U(3,MXN),NEUTRAL(MXN),
* XCOORD(3,MXN), TAU(3)
CHARACTER OUTFILE*12, NEUTRALX*13, SUMMARYX*13
LOGICAL YES
C========> GLBTAU(1,I)=TAUXX = MU*(2DU/DX)+LAMBDA*DIVV
C========> GLBTAU(2,I)=TAUYY = MU*(2DV/DY)+LAMBDA*DIVV
C========> GLBTAU(3,I)=TAUZZ = MU*(2DW/DZ)+LAMBDA*DIVV
C========> GLBTAU(4,I)=TAUXY = MU*(DU/DY+DU/DX)
C========> GLBTAU(5,I)=TAUXZ = MU*(DU/DZ+DW/DX)
C========> GLBTAU(6,I)=TAUYZ = MU*(DV/DZ+DW/DY)
C========> GLBTAU(7,I)= DIVV
C--------0---------0---------0---------0---------0---------0---------0--
C
IF ( ND .EQ. 8 ) THEN
OUTFILE = 'STRESS08.OUT'
NEUTRALX = 'NEUTRAL08.OUT'
SUMMARYX = 'SUMMARY08.OUT'
END IF
IF ( ND .EQ. 27 ) THEN
OUTFILE = 'STRESS27.OUT'
NEUTRALX = 'NEUTRAL27.OUT'
SUMMARYX = 'SUMMARY27.OUT'
END IF
C
OPEN ( 1, FILE=OUTFILE, STATUS='UNKNOWN' )
WRITE (1,*) '# X Y Z U V W TXX TYY TZZ TXY TXZ TYZ DIVV',
* ' CK T1 T2 T3 MSS'
C
ROOT22 = DSQRT(2.D0)/2.D0
DO I = 1 , NNODE
C----- PRINCIPAL STRESSES TAU(3) ----------
CALL PRINCIPAL (GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I),
* GLBTAU(4,I),GLBTAU(5,I),GLBTAU(6,I), IC, TAU )
C----- MISES STRESS -------
YMISES = DSQRT ((TAU(1)-TAU(2))**2 + (TAU(1)-TAU(3))**2 +
* (TAU(2)-TAU(3))**2 ) * ROOT22
WRITE (1,*) I, (XCOORD(J,I),J=1,3), (U(K,I),K=1,3),
* (GLBTAU(J,I),J=1,7),IC, (TAU(J),J=1,3), YMISES
END DO
CLOSE (1)
C-----------------------
OPEN (1,FILE=NEUTRALX, STATUS='UNKNOWN')
WRITE(1,*) 'I X Y Z U V W TXX TYY TZZ'
DO J = 1 , NNEUT
I = NEUTRAL(J)
WRITE(1,*) I, (XCOORD(K,I),K=1,3), (U(K,I),K=1,3),
* GLBTAU(1,I),GLBTAU(2,I),GLBTAU(3,I)
END DO
CLOSE (1)
C--------
OPEN (1,FILE=SUMMARYX, STATUS='UNKNOWN')
CALL COMPUTER (MXN,NNEUT,NEUTRAL,XCOORD,U)
CLOSE (1)
RETURN
END
C
C
SUBROUTINE DERIVEPT ( ND, E )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION E(3,ND)
IF ( ND .EQ. 8 ) THEN
E(1, 1) =-1.D0
E(1, 2) = 1.D0
E(1, 3) = 1.D0
E(1, 4) =-1.D0
E(1, 5) =-1.D0
E(1, 6) = 1.D0
E(1, 7) = 1.D0
E(1, 8) =-1.D0
E(2, 1) =-1.D0
E(2, 2) =-1.D0
E(2, 3) = 1.D0
E(2, 4) = 1.D0
E(2, 5) =-1.D0
E(2, 6) =-1.D0
E(2, 7) = 1.D0
E(2, 8) = 1.D0
E(3, 1) =-1.D0
E(3, 2) =-1.D0
E(3, 3) =-1.D0
E(3, 4) =-1.D0
E(3, 5) = 1.D0
E(3, 6) = 1.D0
E(3, 7) = 1.D0
E(3, 8) = 1.D0
RETURN
END IF
IF ( ND .EQ. 27 ) THEN
E(1, 1) =-1.D0
E(1, 2) = 0.D0
E(1, 3) = 1.D0
E(1, 4) = 1.D0
E(1, 5) = 1.D0
E(1, 6) = 0.D0
E(1, 7) =-1.D0
E(1, 8) =-1.D0
E(1, 9) = 0.D0
E(1,10) =-1.D0
E(1,11) = 0.D0
E(1,12) = 1.D0
E(1,13) = 1.D0
E(1,14) = 1.D0
E(1,15) = 0.D0
E(1,16) =-1.D0
E(1,17) =-1.D0
E(1,18) = 0.D0
E(1,19) =-1.D0
E(1,20) = 0.D0
E(1,21) = 1.D0
E(1,22) = 1.D0
E(1,23) = 1.D0
E(1,24) = 0.D0
E(1,25) =-1.D0
E(1,26) =-1.D0
E(1,27) = 0.D0
E(2, 1) =-1.D0
E(2, 2) =-1.D0
E(2, 3) =-1.D0
E(2, 4) = 0.D0
E(2, 5) = 1.D0
E(2, 6) = 1.D0
E(2, 7) = 1.D0
E(2, 8) = 0.D0
E(2, 9) = 0.D0
E(2,10) =-1.D0
E(2,11) =-1.D0
E(2,12) =-1.D0
E(2,13) = 0.D0
E(2,14) = 1.D0
E(2,15) = 1.D0
E(2,16) = 1.D0
E(2,17) = 0.D0
E(2,18) = 0.D0
E(2,19) =-1.D0
E(2,20) =-1.D0
E(2,21) =-1.D0
E(2,22) = 0.D0
E(2,23) = 1.D0
E(2,24) = 1.D0
E(2,25) = 1.D0
E(2,26) = 0.D0
E(2,27) = 0.D0
E(3, 1) =-1.D0
E(3, 2) =-1.D0
E(3, 3) =-1.D0
E(3, 4) =-1.D0
E(3, 5) =-1.D0
E(3, 6) =-1.D0
E(3, 7) =-1.D0
E(3, 8) =-1.D0
E(3, 9) =-1.D0
E(3,10) = 0.D0
E(3,11) = 0.D0
E(3,12) = 0.D0
E(3,13) = 0.D0
E(3,14) = 0.D0
E(3,15) = 0.D0
E(3,16) = 0.D0
E(3,17) = 0.D0
E(3,18) = 0.D0
E(3,19) = 1.D0
E(3,20) = 1.D0
E(3,21) = 1.D0
E(3,22) = 1.D0
E(3,23) = 1.D0
E(3,24) = 1.D0
E(3,25) = 1.D0
E(3,26) = 1.D0
E(3,27) = 1.D0
RETURN
END IF
STOP 'ND ERROR AT DERIVEPT'
END
C
C
SUBROUTINE DERIVA ( ND, E, F1, F2, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION BPP(3,ND,ND),F1(ND),F2(ND),E(3,ND)
DATA DL / 0.1D0 /
NUMPOINT = ND
C------- COMPUTATION OF BPP(J) = D N(J) / D ETA1
DO K = 1 , NUMPOINT
CALL ISOPARA ( ND , E(1,K)+DL , E(2,K) , E(3,K), F1 )
CALL ISOPARA ( ND , E(1,K)-DL , E(2,K) , E(3,K), F2 )
DO I = 1 , ND
BPP(1,I,K) = (F1(I) - F2(I))/(2.D0*DL)
END DO
C------- COMPUTATION OF BPP(J) = D N(J) / D ETA2
CALL ISOPARA ( ND , E(1,K) , E(2,K)+DL , E(3,K), F1 )
CALL ISOPARA ( ND , E(1,K) , E(2,K)-DL , E(3,K), F2 )
DO I = 1 , ND
BPP(2,I,K) = (F1(I) - F2(I))/(2.D0*DL)
END DO
C------- COMPUTATION OF BPP(J) = D N(J) / D ETA3
CALL ISOPARA ( ND , E(1,K) , E(2,K) , E(3,K)+DL, F1 )
CALL ISOPARA ( ND , E(1,K) , E(2,K) , E(3,K)-DL, F2 )
DO I = 1 , ND
BPP(3,I,K) = (F1(I) - F2(I))/(2.D0*DL)
END DO
END DO
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND, E1 , E2 , E3 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
DATA C / 1.D0 /
IF ( ND .EQ. 8 ) THEN
F(1) = 0.125D0*(1.D0- E1)*(1.D0- E2)*(1.D0- E3)
F(2) = 0.125D0*(1.D0+ E1)*(1.D0- E2)*(1.D0- E3)
F(3) = 0.125D0*(1.D0+ E1)*(1.D0+ E2)*(1.D0- E3)
F(4) = 0.125D0*(1.D0- E1)*(1.D0+ E2)*(1.D0- E3)
F(5) = 0.125D0*(1.D0- E1)*(1.D0- E2)*(1.D0+ E3)
F(6) = 0.125D0*(1.D0+ E1)*(1.D0- E2)*(1.D0+ E3)
F(7) = 0.125D0*(1.D0+ E1)*(1.D0+ E2)*(1.D0+ E3)
F(8) = 0.125D0*(1.D0- E1)*(1.D0+ E2)*(1.D0+ E3)
RETURN
END IF
IF ( ND .EQ. 27 ) THEN
C--------0---------0---------0---------0---------0---------0---------0--
C------ BOTTOM
F( 1)=-0.125D0*(C-E1)*(C-E2)*(C-E3)*E1*E2*E3
F( 2)= 0.25D0*(C-E1)*(C-E2)*(C-E3) *E2*E3*(C+E1)
F( 3)= 0.125D0 *(C-E2)*(C-E3)*E1*E2*E3*(C+E1)
F( 4)= -0.25D0 *(C-E2)*(C-E3)*E1 *E3*(C+E1)*(C+E2)
F( 5)=-0.125D0 *(C-E3)*E1*E2*E3*(C+E1)*(C+E2)
F( 6)= -0.25D0*(C-E1) *(C-E3) *E2*E3*(C+E1)*(C+E2)
F( 7)= 0.125D0*(C-E1) *(C-E3)*E1*E2*E3 *(C+E2)
F( 8)= 0.25D0*(C-E1)*(C-E2)*(C-E3)*E1 *E3 *(C+E2)
F( 9)= -0.5D0*(C-E1)*(C-E2)*(C-E3) *E3*(C+E1)*(C+E2)
C------ MID
F(10)= 0.25D0*(C-E1)*(C-E2)*(C-E3)*E1*E2 *(C+E3)
F(11)= -0.5D0*(C-E1)*(C-E2)*(C-E3) *E2 *(C+E1) *(C+E3)
F(12)= -0.25D0 *(C-E2)*(C-E3)*E1*E2 *(C+E1) *(C+E3)
F(13)= 0.5D0 *(C-E2)*(C-E3)*E1 *(C+E1)*(C+E2)*(C+E3)
F(14)= 0.25D0 *(C-E3)*E1*E2 *(C+E1)*(C+E2)*(C+E3)
F(15)= 0.5D0*(C-E1) *(C-E3) *E2 *(C+E1)*(C+E2)*(C+E3)
F(16)= -0.25D0*(C-E1) *(C-E3)*E1*E2 *(C+E2)*(C+E3)
F(17)= -0.5D0*(C-E1)*(C-E2)*(C-E3)*E1 *(C+E2)*(C+E3)
F(18)= 1.D0*(C-E1)*(C-E2)*(C-E3) *(C+E1)*(C+E2)*(C+E3)
C------ TOP
F(19)= 0.125D0*(C-E1)*(C-E2) *E1*E2*E3 *(C+E3)
F(20)= -0.25D0*(C-E1)*(C-E2) *E2*E3*(C+E1) *(C+E3)
F(21)=-0.125D0 *(C-E2) *E1*E2*E3*(C+E1) *(C+E3)
F(22)= 0.25D0 *(C-E2) *E1 *E3*(C+E1)*(C+E2)*(C+E3)
F(23)= 0.125D0 *E1*E2*E3*(C+E1)*(C+E2)*(C+E3)
F(24)= 0.25D0*(C-E1) *E2*E3*(C+E1)*(C+E2)*(C+E3)
F(25)=-0.125D0*(C-E1) *E1*E2*E3 *(C+E2)*(C+E3)
F(26)= -0.25D0*(C-E1)*(C-E2) *E1 *E3 *(C+E2)*(C+E3)
F(27)= 0.5D0*(C-E1)*(C-E2) *E3*(C+E1)*(C+E2)*(C+E3)
RETURN
END IF
STOP 'NO ELEMENT SELECTED'
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,
* NODEX,U,TAUND,GLBTAU, ICONNECT,BPP,FLMDA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION YI(3,3),YAC(3,3),DU(3,3)
DIMENSION BX(3,ND),TAUND(7,ND), BPP(3,ND,ND)
DIMENSION NODEX(MXE,ND)
DIMENSION XCOORD(3,MXN),U(3,MXN),GLBTAU(7,MXN),ICONNECT(MXN)
C--------0---------0---------0---------0---------0---------0---------0--
C-------- COMPUTATION OF STRESS ======== RETURN VALUES ====> GLBTAU
C-------- DERIVATIVES ARE EVALUATED AT NODAL POINTS.
C========> GLBTAU(1,I)=TAUXX = MU*(2DU/DX)+LAMBDA*DIVV
C========> GLBTAU(2,I)=TAUYY = MU*(2DV/DY)+LAMBDA*DIVV
C========> GLBTAU(3,I)=TAUZZ = MU*(2DW/DZ)+LAMBDA*DIVV
C========> GLBTAU(4,I)=TAUXY = MU*(DU/DY+DU/DX)
C========> GLBTAU(5,I)=TAUXZ = MU*(DU/DZ+DW/DX)
C========> GLBTAU(6,I)=TAUYZ = MU*(DV/DZ+DW/DY)
C-----------------------------------
NUMPOINT = ND
C-----------------------------------
DO K = 1 , 7
DO I = 1 , NNODE
GLBTAU(K,I) = 0.D0
END DO
END DO
C-----------------------------------
DO 500 IEL = 1 , NE
DO 400 K = 1 , NUMPOINT
DO I = 1 , 3
DO J = 1 , 3
YAC(I,J) = 0.D0
END DO
END DO
C
DO J = 1 , ND
NODE = NODEX(IEL,J)
DO I = 1 , 3
YAC(I,1) = YAC(I,1) + BPP(I,J,K) * XCOORD(1,NODE)
YAC(I,2) = YAC(I,2) + BPP(I,J,K) * XCOORD(2,NODE)
YAC(I,3) = YAC(I,3) + BPP(I,J,K) * XCOORD(3,NODE)
END DO
END DO
C
DETJ = YAC(1,1) * ( YAC(2,2)*YAC(3,3)-YAC(2,3)*YAC(3,2) )
* - YAC(1,2) * ( YAC(2,1)*YAC(3,3)-YAC(2,3)*YAC(3,1) )
* + YAC(1,3) * ( YAC(2,1)*YAC(3,2)-YAC(2,2)*YAC(3,1) )
IF ( DETJ .EQ. 0.D0 ) DETJ = 1.0D-15
C------- INVERSE OF THE JACOBIAN MATRIX
YI(1,1) = (YAC(2,2)*YAC(3,3) - YAC(2,3)*YAC(3,2))/DETJ
YI(2,1) = (YAC(2,3)*YAC(3,1) - YAC(2,1)*YAC(3,3))/DETJ
YI(3,1) = (YAC(2,1)*YAC(3,2) - YAC(2,2)*YAC(3,1))/DETJ
YI(1,2) = (YAC(1,3)*YAC(3,2) - YAC(1,2)*YAC(3,3))/DETJ
YI(2,2) = (YAC(1,1)*YAC(3,3) - YAC(1,3)*YAC(3,1))/DETJ
YI(3,2) = (YAC(1,2)*YAC(3,1) - YAC(1,1)*YAC(3,2))/DETJ
YI(1,3) = (YAC(1,2)*YAC(2,3) - YAC(1,3)*YAC(2,2))/DETJ
YI(2,3) = (YAC(1,3)*YAC(2,1) - YAC(1,1)*YAC(2,3))/DETJ
YI(3,3) = (YAC(1,1)*YAC(2,2) - YAC(1,2)*YAC(2,1))/DETJ
DO J = 1 , ND
DO I = 1 , 3
BX(I,J) = YI(I,1)*BPP(1,J,K) + YI(I,2)*BPP(2,J,K)
* + YI(I,3) * BPP(3,J,K)
END DO
END DO
C-------- DU(I,J) REPRESENTS DERIVATIVE OF I COMPONENT WITH RESPECT TO
C-------- J INDEPENDENT VARIABLE.
DO I = 1 , 3
DO J = 1 , 3
DU(I,J) = 0.D0
END DO
END DO
DO L = 1 , ND
NODE = NODEX(IEL,L)
DO I = 1 , 3
DO J = 1 , 3
DU(I,J) = DU(I,J) + U(I,NODE)*BX(J,L)
END DO
END DO
END DO
DIVV = DU(1,1)+DU(2,2)+DU(3,3)
TAUND(1,K)= VISCO*(DU(1,1)+DU(1,1)) + FLMDA*DIVV
TAUND(2,K)= VISCO*(DU(2,2)+DU(2,2)) + FLMDA*DIVV
TAUND(3,K)= VISCO*(DU(3,3)+DU(3,3)) + FLMDA*DIVV
TAUND(4,K)= VISCO*(DU(1,2)+DU(2,1))
TAUND(5,K)= VISCO*(DU(1,3)+DU(3,1))
TAUND(6,K)= VISCO*(DU(2,3)+DU(3,2))
TAUND(7,K)= DIVV
400 CONTINUE
C
DO I = 1 , ND
NODE = NODEX(IEL,I)
DO K = 1 , 7
GLBTAU(K,NODE) = GLBTAU(K,NODE) + TAUND(K,I)
END DO
END DO
500 CONTINUE
DO I = 1 , NNODE
DO K = 1 , 7
GLBTAU(K,I) = GLBTAU(K,I) / ICONNECT(I)
END DO
END DO
RETURN
END
C
C
SUBROUTINE INPUT ( INPFILE,MXE,MXN,MXB,ND,NE,NNODE,
* YOUNG,POISSON, NODEX,XCOORD,NNEUT, NEUTRAL,NCTOP,NODECTOP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(3,MXN),NEUTRAL(MXN)
DIMENSION NODECTOP(MXN)
CHARACTER*14 INPFILE
LOGICAL YES
C========> FILENAME INPFILE SEE MAIN PROGRAM
WRITE (*,*) INPFILE
OPEN ( 1, FILE=INPFILE, STATUS='UNKNOWN' )
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(1,NODE),XCOORD(2,NODE),XCOORD(3,NODE)
END DO
CLOSE (1)
C---------- NEUTRAL NODAL NUMBERS
OPEN ( 1, FILE='NEUTRAL.DAT', STATUS = 'UNKNOWN')
READ (1,*) NNEUT
DO I = 1 , NNEUT
READ (1,*) NEUTRAL(I)
END DO
CLOSE (1)
C---------- NODAL NUMBERS OF ALONG THE CENTER OF THE TOP FACE(+Z)
OPEN ( 1, FILE='NODECTOP.DAT', STATUS = 'UNKNOWN')
READ (1,*) NCTOP
DO I = 1 , NCTOP
READ (1,*) NODECTOP(I)
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INPUTUV ( MXN,NNODE,U )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION U(3,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(1,I) , I = 1 , NNODE )
READ (1) ( U(2,I) , I = 1 , NNODE )
READ (1) ( U(3,I) , I = 1 , NNODE )
CLOSE (1)
WRITE (*,*) ' DISPLACE.MNT READING DONE'
RETURN
END
C
C
SUBROUTINE PRINCIPAL (TXX,TYY,TZZ,TXY,TXZ,TYZ, IC, TAU )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION TAU(3),DRV(3,3),DI(3,3),F(3)
DATA EPS / 1.D-14 /
DATA MAXITA / 100 /
C--------------- STRESS INVARIANTS ---------------
A = TXX+TYY+TZZ
B = TXX*TYY+TYY*TZZ+TZZ*TXX-TXY**2-TYZ**2-TXZ**2
C = TXX*TYY*TZZ+2.D0*TXY*TXZ*TYZ-
* TXY**2*TZZ-TYZ**2*TXX-TXZ**2*TYY
C--------- ENERGY BY COMPUTAED GIVEN STRESSES ------
ENERGY0 = (TXX-TYY)**2+(TYY-TZZ)**2+(TZZ-TXX)**2+
* 6.D0*(TXY**2+TYZ**2+TXZ**2)
C-------------- CASE OF A=B=C=0 -------------
IC = 0
TAU(1) = 0.D0
TAU(2) = 0.D0
TAU(3) = 0.D0
IF ( A.EQ.0.D0 .AND. B.EQ.0.D0 .AND. C.EQ.0.D0 ) RETURN
IF ( ENERGY0 .EQ. 0.D0 ) RETURN
C------------ CASE OF B=C=0 -----------------
IC = -1
IF ( B .EQ. 0.D0 .AND. C .EQ. 0.D0 ) THEN
TAU(2) = 0.D0
IF ( A .GT. 0.D0 ) THEN
TAU(1) = A
TAU(3) = 0.D0
ELSE
TAU(1) = 0.D0
TAU(3) = A
END IF
RETURN
END IF
C---------------- CASE OF C=0 -----------------
IC = -2
IF ( C .EQ. 0.D0 ) THEN
TAU(1) = (A+DSQRT(A*A-4.D0*B))/2.D0
TAU(2) = 0.D0
TAU(3) = (A-DSQRT(A*A-4.D0*B))/2.D0
IF ( TAU(3) .GT. 0.D0 ) THEN
TAU(2) = TAU(3)
TAU(3) = 0.D0
END IF
RETURN
END IF
C----------------- CASE OF A, B, C NOT ZERO -----------
IC = -3
C------------ INITIAL GUESSES FOR TAU'S -------------
TAU(1) = DSQRT(2.D0) + A/3.D0
TAU(2) = DSQRT(3.D0) + B/TAU(1)
TAU(3) = DATAN(1.D0) + C/(TAU(1)*TAU(2))
C
C============== NEWTON-RAPHSON METHOD ITERATION ==============
DO I = 1 , MAXITA
F(1) = TAU(1)+TAU(2)+TAU(3)-A
F(2) = TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1)-B
F(3) = TAU(1)*TAU(2)*TAU(3)-C
C--------- DERIVATIVE MATRIX [DRV] ------
DRV(1,1) = 1.D0
DRV(1,2) = 1.D0
DRV(1,3) = 1.D0
DRV(2,1) = TAU(2)+TAU(3)
DRV(2,2) = TAU(1)+TAU(3)
DRV(2,3) = TAU(1)+TAU(2)
DRV(3,1) = TAU(2)*TAU(3)
DRV(3,2) = TAU(1)*TAU(3)
DRV(3,3) = TAU(1)*TAU(2)
C------- DETERMINANT OF [DRV]
DETJ = DRV(1,1) * ( DRV(2,2)*DRV(3,3)-DRV(2,3)*DRV(3,2) )
* - DRV(1,2) * ( DRV(2,1)*DRV(3,3)-DRV(2,3)*DRV(3,1) )
* + DRV(1,3) * ( DRV(2,1)*DRV(3,2)-DRV(2,2)*DRV(3,1) )
C------- [DI] = INVERSE OF THE MATRIX [DRV]
DI(1,1) = (DRV(2,2)*DRV(3,3) - DRV(2,3)*DRV(3,2))/DETJ
DI(2,1) = (DRV(2,3)*DRV(3,1) - DRV(2,1)*DRV(3,3))/DETJ
DI(3,1) = (DRV(2,1)*DRV(3,2) - DRV(2,2)*DRV(3,1))/DETJ
DI(1,2) = (DRV(1,3)*DRV(3,2) - DRV(1,2)*DRV(3,3))/DETJ
DI(2,2) = (DRV(1,1)*DRV(3,3) - DRV(1,3)*DRV(3,1))/DETJ
DI(3,2) = (DRV(1,2)*DRV(3,1) - DRV(1,1)*DRV(3,2))/DETJ
DI(1,3) = (DRV(1,2)*DRV(2,3) - DRV(1,3)*DRV(2,2))/DETJ
DI(2,3) = (DRV(1,3)*DRV(2,1) - DRV(1,1)*DRV(2,3))/DETJ
DI(3,3) = (DRV(1,1)*DRV(2,2) - DRV(1,2)*DRV(2,1))/DETJ
C-------- NEW VALUES FOR TAU'S
TAU(1) = TAU(1) - ( DI(1,1)*F(1)+DI(1,2)*F(2)+DI(1,3)*F(3) )
TAU(2) = TAU(2) - ( DI(2,1)*F(1)+DI(2,2)*F(2)+DI(2,3)*F(3) )
TAU(3) = TAU(3) - ( DI(3,1)*F(1)+DI(3,2)*F(2)+DI(3,3)*F(3) )
C--------- ENERGY BY COMPUTAED PRINCIPAL STRESSES ------
ENERGY1 = 2.D0*(TAU(1)+TAU(2)+TAU(3))**2-
* 6.D0*(TAU(1)*TAU(2)+TAU(2)*TAU(3)+TAU(3)*TAU(1))
C--------- DECISION --------
IF ( DABS(ENERGY1-ENERGY0)/ENERGY0 .LE. EPS ) THEN
C------- MAKING TAU(1), TAU(2), AND TAU(3) DESCENDING ORDER
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
END IF
IF ( TAU(3) .GT. TAU(2) ) THEN
TEMP = TAU(2)
TAU(2) = TAU(3)
TAU(3) = TEMP
END IF
IF ( TAU(2) .GT. TAU(1) ) THEN
TEMP = TAU(1)
TAU(1) = TAU(2)
TAU(2) = TEMP
END IF
IC = I
RETURN
END IF
END DO
C============= CASE OF NO CONVERGENCE =============
C-- ONE OF OR ALL MAY BE COMPLEX.
TAU(1) = 0.D0
TAU(2) = 0.D0
TAU(3) = 0.D0
IC = -4
RETURN
END
C
C
SUBROUTINE COMPUTER (MXN,NNEUT,NEUTRAL,XCOORD,U)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NEUTRAL(MXN), XCOORD(3,MXN), U(3,MXN)
LAST = NEUTRAL(NNEUT)
WRITE (1,*) 'END DISPLACEMENT = ', U(3,LAST)
C-------- INITIAL VALUE FOR R ---------
R = DABS((XCOORD(1,LAST)**2-U(3,LAST)**2)/(2.D0*U(3,LAST)))
WRITE (1,*) 'R = ', R
DR = 1.D0
N = 0.9*NNEUT
DO K = 1 , 5
F =0.D0
DO J = 1 , N
Z = R-DSQRT(R*R-XCOORD(1,NEUTRAL(J))**2)
F = F + (Z-U(3,NEUTRAL(J)))*Z/(R-Z)
END DO
G = 0.D0
DO J = 1 , N
Z = (R+DR)-DSQRT((R+DR)*(R+DR)-XCOORD(1,NEUTRAL(J))**2)
G = G + (Z-U(3,NEUTRAL(J)))*Z/((R+DR)-Z)
END DO
DFDR= (G-F)/DR
R = R - F/DFDR
WRITE (1,*) 'R = ', R
END DO
RETURN
END
C
C
SUBROUTINE ANGLE ( MXN, NNEUT, NEUTRAL, NCTOP,NODECTOP,XCOORD,U )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NEUTRAL(MXN), XCOORD(3,MXN), U(3,MXN),NODECTOP(MXN)
OPEN ( 1, FILE='ANGLE.DAT', STATUS = 'UNKNOWN')
WRITE (1,*) 'X-COORDINATE ARM DISPLACEMENT-TOP-CENTER',
* ' DISPLACEMENT-NUTRAL TWIST RATE-OF-TWIST'
I = 1
ARM = XCOORD(3,NODECTOP(I)) - XCOORD(3,NEUTRAL(I))
DISPVTOP = U(2,NODECTOP(I))
DISPNTRL = U(2,NEUTRAL(I))
DISTANCE = XCOORD(1,NEUTRAL(I))
TWIST = DATAN ( DABS(DISPVTOP)/ARM )
WRITE (1,*) DISTANCE, ARM, DISPVTOP, DISPNTRL, TWIST
X0 = DISTANCE
V0 = DISPVTOP
DO I = 2 , NNEUT
ARM = XCOORD(3,NODECTOP(I)) - XCOORD(3,NEUTRAL(I))
DISPVTOP = U(2,NODECTOP(I))
DISPNTRL = U(2,NEUTRAL(I))
DISTANCE = XCOORD(1,NEUTRAL(I))
TWIST = DATAN ( DABS(DISPVTOP)/ARM )
X1 = DISTANCE
V1 = DISPVTOP
RTWIST = ( V1 - V0 )/( X1 - X0 )
X0 = X1
V0 = V1
WRITE (1,*) DISTANCE, ARM, DISPVTOP,DISPNTRL,TWIST,ABS(RTWIST)
END DO
CLOSE (1)
RETURN
END