! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE B_MATRIX_AXISYM_ELASTIC (NOD_PER_EL, N_SPACE, & N_G_DOF, GDH, H, R, NS, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! AXISYMMETRIC ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF, NS REAL(DP), INTENT(IN) :: R, GDH (N_SPACE, NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(IN) :: H (NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF) INTEGER :: J, K, L ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! GDH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NS = NUMBER OF STRAINS (ROWS IN B): XX, YY, XY, HOOP ! N_SPACE = DIMENSION OF SPACE ! R,Z,T DENOTE RADIAL, AXIAL, CIRCUMFERENCE B = 0.0d0 DO J = 1, NOD_PER_EL ! ROW NUMBER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U L = K + 1 ! SECOND COLUMN, V B (1, K) = GDH (1, J) ! DU/DX FOR XX NORMAL B (3, K) = GDH (2, J) ! DU/DY FOR XY SHEAR IF ( R <= 0.0d0 ) STOP 'R=0, B_MATRIX_AXISYM_ELASTIC' B (4, K) = H (J) / R ! U/R HOOP, ZZ NORMAL B (2, L) = GDH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = GDH (1, J) ! DV/DX FOR XY SHEAR END DO END SUBROUTINE B_MATRIX_AXISYM_ELASTIC SUBROUTINE B_MATRIX_BAR_ELASTIC (NOD_PER_EL, N_SPACE, N_G_DOF, & GDH, NS, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTIC BAR STRAIN-DISPLACEMENT RELATION (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF, NS REAL(DP), INTENT(IN) :: GDH (N_SPACE, NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF) INTEGER :: J, K ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! GDH = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NS = NUMBER OF STRAINS (ROWS IN B) = 1 ! N_SPACE = DIMENSION OF SPACE DO J = 1, NOD_PER_EL ! ROW NUMBER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U B (1, K) = GDH (1, J) ! DU/DX FOR XX NORMAL END DO END SUBROUTINE B_MATRIX_BAR_ELASTIC SUBROUTINE B_MATRIX_CYL_ELASTIC (NOD_PER_EL, GDH, H, R, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTIC CYLINDER STRAIN-DISPLACEMENT RELATIONS (B) ! STRESS & STRAIN COMPONENT ORDER: RR, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL REAL(DP), INTENT(IN) :: GDH (1, NOD_PER_EL) REAL(DP), INTENT(IN) :: H (NOD_PER_EL), R REAL(DP), INTENT(OUT) :: B (2, NOD_PER_EL) INTEGER :: J ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! GDH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! R,T DENOTE RADIAL, CIRCUMFERENCE B = 0.d0 DO J = 1, NOD_PER_EL ! ROW NUMBER B (1, J) = GDH (1, J) ! DU/DR FOR RR NORMAL IF ( R <= 0.d0 ) STOP 'R=0, B_MATRIX_CYL_ELASTIC' B (2, J) = H (J) / R ! U/R HOOP END DO END SUBROUTINE B_MATRIX_CYL_ELASTIC SUBROUTINE B_MATRIX_ELASTIC (I_OPT, NOD_PER_EL, N_SPACE, & N_G_DOF, GDH, H, R, NS, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: I_OPT, NOD_PER_EL, N_SPACE, N_G_DOF, NS REAL(DP), INTENT(IN) :: R, GDH (N_SPACE, NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(IN) :: H (NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF) INTEGER :: J, K, L, M ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! GDH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! I_OPT = ELASTICITY CLASS ! = 1, AXIAL BAR, N_G_DOF = 1, NS = 1 ! = 2, PLANE STRESS, N_G_DOF = 2, NS = 3 ! = 3, PLANE STRAIN, N_G_DOF = 2, NS = 3 ! = 4, AXISYMMETRIC, N_G_DOF = 2, R = RADIUS, NS = 4 ! = 5, 3-D SOLID, N_G_DOF = 3, NS = 6 ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NS = NUMBER OF STRAINS (ROWS IN B) ! XX, YY, XY, ZZ, XZ, YZ, 1 OR 3 OR 4 OR 6 ! N_SPACE = DIMENSION OF SPACE B = 0.d0 IF ( I_OPT < 1 .OR. I_OPT > 5 ) STOP 'I_OPT ERROR, B_MATRIX_ELASTIC' DO J = 1, NOD_PER_EL ! ROW NUBMER K = N_G_DOF * (J - 1) + 1 ! "FIRST" COLUMN, U !--> ONE-DIMENSIONAL AXIAL BAR, I_OPT = 1 B (1, K) = GDH (1, J) ! DU/DX FOR XX NORMAL IF ( I_OPT == 1 ) CYCLE ! TO NEXT J VALUE !-> PLANE STRESS, PLANE STRAIN, AXISYMMETRIC, 3D ! STRESS & STRAIN COMPONENT ORDER: XX, YY, AND XY L = K + 1 ! "SECOND" COLUMN, V B (3, K) = GDH (2, J) ! DU/DY FOR XY SHEAR B (2, L) = GDH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = GDH (1, J) ! DV/DX FOR XY SHEAR IF ( I_OPT == 2 .OR. I_OPT == 3 ) CYCLE ! TO NEXT J VALUE !-> AXISYMMETRIC ONLY ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, AND TT IF ( I_OPT == 4 ) THEN IF ( R <= 0.d0 ) STOP 'R=0, BELAST' B (4, K) = H (J) / R ! U/R HOOP, ZZ NORMAL CYCLE ! TO NEXT J VALUE END IF ! AXISYMMETRIC !-> 3D SOLID ONLY ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, ZZ, XZ, YZ M = L + 1 ! "THIRD" COLUMN, W B (5, K) = GDH (3, J) ! DU/DZ FOR XZ SHEAR B (6, L) = GDH (3, J) ! DV/DZ FOR YZ SHEAR B (4, M) = GDH (3, J) ! DW/DZ FOR ZZ NORMAL B (5, M) = GDH (1, J) ! DW/DX FOR XZ SHEAR B (6, M) = GDH (2, J) ! DW/DY FOR YZ SHEAR END DO END SUBROUTINE B_MATRIX_ELASTIC SUBROUTINE B_MATRIX_PLANE_ELASTIC (NOD_PER_EL, N_SPACE, & N_G_DOF, GDH, NS, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! 2-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF, NS !b xxx REAL(DP), INTENT(IN) :: GDH (N_SPACE, NOD_PER_EL * N_G_DOF) REAL(DP), INTENT(IN) :: GDH (N_SPACE, NOD_PER_EL) !b REAL(DP), INTENT(OUT) :: B (NS, NOD_PER_EL * N_G_DOF) INTEGER :: J, K, L ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! GDH = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NS = NUMBER OF STRAINS (ROWS IN B) XX, YY, XY ! N_SPACE = DIMENSION OF SPACE !b print *,'size (GDH) ', size (GDH, 1), size (GDH, 2) !b !b print *,'size (B ) ', size (B , 1), size (B , 2) !b DO J = 1, NOD_PER_EL ! ROW NUBMER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U L = K + 1 ! SECOND COLUMN, V B (1, K) = GDH (1, J) ! DU/DX FOR XX NORMAL B (2, K) = 0.d0 B (3, K) = GDH (2, J) ! DU/DY FOR XY SHEAR B (1, L) = 0.d0 B (2, L) = GDH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = GDH (1, J) ! DV/DX FOR XY SHEAR END DO END SUBROUTINE B_MATRIX_PLANE_ELASTIC SUBROUTINE HOOKE (E, STRAIN, STRAIN_0, STRESS, N_R_B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! STRESSES DUE TO INITIAL STRAINS ! STRESS(L) = E(L,M) * ( STRAIN(M) - STRAIN_0(M) ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: ZERO = 0.d0 INTEGER, INTENT(IN) :: N_R_B REAL(DP), INTENT(IN) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: STRAIN (N_R_B), STRAIN_0 (N_R_B) REAL(DP), INTENT(OUT) :: STRESS (N_R_B) ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B MATRIX ! STRAIN = MECHANICAL STRAIN VECTOR ! STRAIN_0 = INITIAL STRAIN VECTOR ! STRESS = STRESS VECTOR STRESS = MATMUL ( E, (STRAIN - STRAIN_0)) END SUBROUTINE HOOKE SUBROUTINE VON_MISES_STRESS (STRESS, PR, I_OPT, EFFECTIVE, NS) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! VON MISES EFFECTIVE STRESS FROM STRESS VECTOR ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module INTEGER, INTENT(IN) :: NS, I_OPT REAL(DP), INTENT(IN) :: STRESS (NS), PR REAL(DP), INTENT(OUT) :: EFFECTIVE REAL(DP), PARAMETER :: A = 0.7071067812d0, SIX = 6.d0 REAL(DP) :: STRESS_4 ! I_OPT = Options: 1=1D, 2=plane stress, 3=plane strain, ! 4=axisymmetric stress, 5=3D solid ! NS = Number of stress components ! PR = Poisson ratio ! STRESS = Stress vector components ! I_OPT = 1, Sxx, Axial ! I_OPT = 2, Sxx Syy Sxy, Plane stress ! I_OPT = 3, Sxx Syy Sxy Szz, Plane strain ! I_OPT = 4, Sxx Syy Sxy Szz, Axisymmetric ! I_OPT = 5, Sxx Syy Sxy Szz Sxz Syz, Full 3-D ! EFFECTIVE = Von Mises Effective Stress SELECT CASE (I_OPT) CASE (1) ! Stress = Sxx, Axial EFFECTIVE = ABS (STRESS (1)) CASE (2) ! Stress = Sxx Syy Sxy, Plane stress EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2 & + STRESS (2) **2 & + STRESS (1) **2 + SIX * STRESS (3) **2) CASE (3) ! Stress = Sxx Syy Sxy Szz, Plane strain STRESS_4 = PR * (STRESS (1) + STRESS (2)) EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2 & + (STRESS (2) - STRESS_4) **2 & + (STRESS_4 - STRESS (1)) **2 & + SIX * STRESS (3) **2) CASE (4) ! Stress = Sxx Syy Sxy Szz, Axisymmetric EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2 & + (STRESS (2) - STRESS (4)) **2 & + (STRESS (4) - STRESS (1)) **2 & + SIX * STRESS (3) **2) CASE (5) ! Stress = Sxx Syy Sxy Szz Sxz Syz, Full 3-D EFFECTIVE = A * SQRT ( (STRESS (1) - STRESS (2)) **2 & + (STRESS (2) - STRESS (4)) **2 & + (STRESS (4) - STRESS (1)) **2 & + SIX * (STRESS (3) **2 + STRESS (5) **2 & + STRESS (6) **2)) CASE DEFAULT STOP 'INVALID I_OPT, VON_MISES_STRESS' END SELECT END SUBROUTINE VON_MISES_STRESS SUBROUTINE AXISYMMETRIC_E_MATRIX (Y_M, P_R, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! AXISYMMETRIC CONSTITUTIVE MATRIX DEFINITION, N_R_B = 4 ! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_R_B IMPLICIT NONE REAL(DP), INTENT(IN) :: Y_M ! YOUNG'S MODULUS REAL(DP), INTENT(IN) :: P_R ! POISSON'S RATIO REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! CONSTITUTIVE REAL(DP) :: C1, C2, C3, S_M ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! S_M = SHEAR MODULUS IF ( N_R_B /= 4 ) STOP 'INVALID b_rows IN AXISYMMETRIC_E_MATRIX' S_M = 0.5d0 * Y_M / ( 1.0d0 + P_R) C1 = Y_M /(1.d0 + P_R)/(1.d0 - P_R - P_R) C2 = C1 * P_R ; C3 = C1 * (1.d0 - P_R) E (:, 3) = 0.d0 ; E (3, :) = 0.d0 E (1, 1) = C3 ; E (2, 1) = C2 ; E (4, 1) = C2 E (1, 2) = C2 ; E (2, 2) = C3 ; E (4, 2) = C2 E (1, 4) = C2 ; E (2, 4) = C2 ; E (4, 4) = C3 E (3, 3) = S_M !b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b END SUBROUTINE AXISYMMETRIC_E_MATRIX SUBROUTINE AXISYMMETRIC_E_MATRIX_GEN (Y_M, P_R, S_M, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! AXISYMMETRIC CONSTITUTIVE MATRIX, INDEPENDENT SHEAR MODULUS ! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_R_B IMPLICIT NONE REAL(DP), INTENT(IN) :: Y_M ! YOUNG'S MODULUS REAL(DP), INTENT(IN) :: P_R ! POISSON'S RATIO REAL(DP), INTENT(IN) :: S_M ! SHEAR MODULUS REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! CONSTITUTIVE REAL(DP) :: C1, C2, C3 ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES IF ( N_R_B /= 4 ) STOP 'INVALID b_rows IN AXISYMMETRIC_E_MATRIX_GEN' C1 = Y_M /(1.d0 + P_R)/(1.d0 - P_R - P_R) C2 = C1 * P_R ; C3 = C1 * (1.d0 - P_R) E (:, 3) = 0.d0 ; E (3, :) = 0.d0 E (1, 1) = C3 ; E (2, 1) = C2 ; E (4, 1) = C2 E (1, 2) = C2 ; E (2, 2) = C3 ; E (4, 2) = C2 E (1, 4) = C2 ; E (2, 4) = C2 ; E (4, 4) = C3 E (3, 3) = S_M !b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b END SUBROUTINE AXISYMMETRIC_E_MATRIX_GEN SUBROUTINE PLANE_STRESS_E_MATRIX (Y, PR, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_R_B IMPLICIT NONE REAL(DP), INTENT(IN) :: Y, PR REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: CONST, SM ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! PR = POISSON'S RATIO ! SM = SHEAR MODULUS ! Y = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN PLANE_STRESS_E_MATRIX' SM = 0.5d0 * Y / ( 1.0d0 + PR) CONST = Y / (1.0d0 - PR * PR) E (1, 1) = CONST ; E (2, 1) = CONST * PR ; E (3, 1) = 0.d0 E (1, 2) = CONST * PR ; E (2, 2) = CONST ; E (3, 2) = 0.d0 E (1, 3) = 0.d0 ; E (2, 3) = 0.d0 ; E (3, 3) = SM !b if ( I_BUG > 0 ) print *, 'plane stress E ' !b !b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b END SUBROUTINE PLANE_STRESS_E_MATRIX SUBROUTINE PLANE_STRESS_E_MATRIX_GEN (Y, PR, SM, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION ! WITH INDEPENDENT SHEAR MODULUS, SM ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_R_B IMPLICIT NONE REAL(DP), INTENT(IN) :: Y, PR, SM REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: CONST ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! PR = POISSON'S RATIO ! SM = SHEAR MODULUS ! Y = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN PLANE_STRESS_E_MATRIX_GEN' CONST = Y / (1.0d0 - PR * PR) E (1, 1) = CONST ; E (2, 1) = CONST * PR ; E (3, 1) = 0.d0 E (1, 2) = CONST * PR ; E (2, 2) = CONST ; E (3, 2) = 0.d0 E (1, 3) = 0.d0 ; E (2, 3) = 0.d0 ; E (3, 3) = SM !b if ( I_BUG > 0 ) print *, 'plane stress E generalized ' !b !b if ( I_BUG > 0 ) call rprint (E, N_R_B, N_R_B, 1) !b END SUBROUTINE PLANE_STRESS_E_MATRIX_GEN SUBROUTINE PLANE_STRAIN_E_MATRIX (Y, PR, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PLANE_STRAIN CONSTITUTIVE MATRIX DEFINITION ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_R_B IMPLICIT NONE REAL(DP), INTENT(IN) :: Y, PR REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: C_1, C_2, SM ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! PR = POISSON'S RATIO ! SM = SHEAR MODULUS ! Y = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN PLANE_STRAIN_E_MATRIX' IF ( PR == 0.5d0 ) STOP 'INVALID PROPERTY IN PLANE_STRAIN_E_MATRIX' C_1 = Y * (1.d0 - PR) / ((1.d0 + PR)*(1.d0 - PR - PR)) C_2 = Y * PR / ((1.d0 + PR)*(1.d0 - PR - PR)) SM = 0.5d0 * Y / ( 1.0d0 + PR) E (1, 1) = C_1 ; E (2, 1) = C_2 ; E (3, 1) = 0.d0 E (1, 2) = C_2 ; E (2, 2) = C_1 ; E (3, 2) = 0.d0 E (1, 3) = 0.d0 ; E (2, 3) = 0.d0 ; E (3, 3) = SM END SUBROUTINE PLANE_STRAIN_E_MATRIX SUBROUTINE ELASTIC_B_AXISYMMETRIC (DGH, R, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! AXISYMMETRIC ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE Use Elem_Type_Data ! for LT_FREE, LT_N, H (LT_N) IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (N_SPACE, LT_N) ! Gradients REAL(DP), INTENT(IN) :: R ! Radius REAL(DP), INTENT(OUT) :: B (N_R_B, LT_N * N_G_DOF) INTEGER :: J, K, L ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE = 2 here ! N_R_B = NUMBER OF STRAINS (ROWS IN B) = 4: XX, YY, XY, HOOP ! N_SPACE = DIMENSION OF SPACE = 2 here ! R,Z,T DENOTE RADIAL, AXIAL, CIRCUMFERENCE B = 0.d0 DO J = 1, LT_N ! ROW NUMBER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U L = K + 1 ! SECOND COLUMN, V B (1, K) = DGH (1, J) ! DU/DX FOR XX NORMAL B (3, K) = DGH (2, J) ! DU/DY FOR XY SHEAR IF ( R <= 0.d0 ) STOP 'R=0, IN ELASTIC_B_AXISYMMETRIC' B (4, K) = H (J) / R ! U/R HOOP, ZZ NORMAL B (2, L) = DGH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = DGH (1, J) ! DV/DX FOR XY SHEAR END DO END SUBROUTINE ELASTIC_B_AXISYMMETRIC SUBROUTINE ELASTIC_B_BAR (DGH, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTIC BAR STRAIN-DISPLACEMENT RELATION (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE Use Elem_Type_Data ! for LT_FREE, LT_N IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL(DP), INTENT(OUT) :: B (N_R_B, LT_N * N_G_DOF) INTEGER :: J, K ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! N_R_B = NUMBER OF STRAINS (ROWS IN B) = 1 ! N_SPACE = DIMENSION OF SPACE DO J = 1, LT_N ! ROW NUMBER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U B (1, K) = DGH (1, J) ! DU/DX FOR XX NORMAL END DO END SUBROUTINE ELASTIC_B_BAR SUBROUTINE ELASTIC_B_CYLINDER (DGH, R, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTIC CYLINDER STRAIN-DISPLACEMENT RELATIONS (B) ! STRESS & STRAIN COMPONENT ORDER: RR, AND TT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE Use Elem_Type_Data ! for LT_FREE, LT_N, H (LT_N) IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (1, LT_N), R REAL(DP), INTENT(OUT) :: B (2, LT_N) INTEGER :: J ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! R,T DENOTE RADIAL, CIRCUMFERENCE B = 0.d0 DO J = 1, LT_N ! ROW NUMBER B (1, J) = DGH (1, J) ! DU/DR FOR RR NORMAL IF ( R <= 0.d0 ) STOP 'R=0, IN ELASTIC_B_CYLINDER' B (2, J) = H (J) / R ! U/R HOOP END DO END SUBROUTINE ELASTIC_B_CYLINDER SUBROUTINE ELASTIC_B_PLANAR (DGH, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! 2-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE Use Elem_Type_Data ! for LT_FREE, LT_N IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL(DP), INTENT(OUT) :: B (N_R_B, LT_N * N_G_DOF) INTEGER :: J, K, L ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE = 2 ! N_R_B = NUMBER OF STRAINS (ROWS IN B) = 3: XX, YY, XY ! N_SPACE = DIMENSION OF SPACE = 2 here DO J = 1, LT_N ! ROW NUBMER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U L = K + 1 ! SECOND COLUMN, V B (1, K) = DGH (1, J) ! DU/DX FOR XX NORMAL B (2, K) = 0.d0 B (3, K) = DGH (2, J) ! DU/DY FOR XY SHEAR B (1, L) = 0.d0 B (2, L) = DGH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = DGH (1, J) ! DV/DX FOR XY SHEAR END DO END SUBROUTINE ELASTIC_B_PLANAR SUBROUTINE ELASTIC_B_SOLID (DGH, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! 3-D ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE Use Elem_Type_Data ! for LT_FREE, LT_N IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL(DP), INTENT(OUT) :: B (N_R_B, LT_N * N_G_DOF) INTEGER :: J, K, L, M ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF ELEM INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! N_R_B = NUMBER OF STRAINS (ROWS IN B) XX, YY, XY, ZZ, XZ, YZ ! N_SPACE = DIMENSION OF SPACE B = 0.d0 DO J = 1, LT_N ! ROW NUMBER K = N_G_DOF * (J - 1) + 1 ! FIRST COLUMN, U TERMS L = K + 1 ! SECOND COLUMN, V TERMS M = L + 1 ! THIRD COLUMN, W TERMS B (1, K) = DGH (1, J) ! DU/DX FOR XX NORMAL B (3, K) = DGH (2, J) ! DU/DY FOR XY SHEAR B (5, K) = DGH (3, J) ! DU/DZ FOR XZ SHEAR B (2, L) = DGH (2, J) ! DV/DY FOR YY NORMAL B (3, L) = DGH (1, J) ! DV/DX FOR XY SHEAR B (6, L) = DGH (3, J) ! DV/DZ FOR YZ SHEAR B (4, M) = DGH (3, J) ! DW/DZ FOR ZZ NORMAL B (5, M) = DGH (1, J) ! DW/DX FOR XZ SHEAR B (6, M) = DGH (2, J) ! DW/DY FOR YZ SHEAR END DO END SUBROUTINE ELASTIC_B_SOLID SUBROUTINE E_PLANE_STRESS (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PLANE_STRESS CONSTITUTIVE MATRIX DEFINITION ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3 ! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B Use Sys_Properties_Data ! for function GET_REAL_LP IMPLICIT NONE REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: CONST, P_R, S_M, Y_M ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! P_R = POISSON'S RATIO ! S_M = SHEAR MODULUS ! Y_M = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN E_PLANE_STRESS' IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_PLANE_STRAIN' Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2) S_M = 0.5d0 * Y_M / (1 + P_R) CONST = Y_M / (1 - P_R * P_R) E (1, 1) = CONST ; E (2, 1) = CONST * P_R ; E (3, 1) = 0.d0 E (1, 2) = CONST * P_R ; E (2, 2) = CONST ; E (3, 2) = 0.d0 E (1, 3) = 0.d0 ; E (2, 3) = 0.d0 ; E (3, 3) = S_M END SUBROUTINE E_PLANE_STRESS SUBROUTINE E_PLANE_STRAIN (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PLANE_STRAIN CONSTITUTIVE MATRIX DEFINITION ! STRESS & STRAIN COMPONENT ORDER: XX, YY, XY, SO N_R_B = 3 ! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B Use Sys_Properties_Data ! for function GET_REAL_LP IMPLICIT NONE REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: C_1, C_2, P_R, S_M, Y_M ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! P_R = POISSON'S RATIO ! S_M = SHEAR MODULUS ! Y_M = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 3 ) STOP 'INVALID b_rows IN E_PLANE_STRAIN' IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_PLANE_STRAIN' Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2) IF ( P_R == 0.5d0 ) STOP "INVALID POISSON'S RATIO IN E_PLANE_STRAIN" C_1 = Y_M * (1 - P_R) / ((1 + P_R)*(1 - P_R - P_R)) C_2 = Y_M * P_R / ((1 + P_R)*(1 - P_R - P_R)) S_M = 0.5d0 * Y_M / ( 1 + P_R) E (1, 1) = C_1 ; E (2, 1) = C_2 ; E (3, 1) = 0.d0 E (1, 2) = C_2 ; E (2, 2) = C_1 ; E (3, 2) = 0.d0 E (1, 3) = 0.d0 ; E (2, 3) = 0.d0 ; E (3, 3) = S_M END SUBROUTINE E_PLANE_STRAIN SUBROUTINE E_AXISYMMETRIC_STRESS (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! AXISYMMETRIC CONSTITUTIVE MATRIX DEFINITION, N_R_B = 4 ! STRESS & STRAIN COMPONENT ORDER: RR, ZZ, RZ, AND TT ! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B Use Sys_Properties_Data ! for function GET_REAL_LP IMPLICIT NONE REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! CONSTITUTIVE REAL(DP) :: C_1, C_2, C_3, P_R, S_M, Y_M ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! P_R = POISSON'S RATIO ! S_M = SHEAR MODULUS ! Y_M = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 4 ) STOP 'INVALID b_rows IN E_AXISYMMETRIC_STRESS' IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_AXISYMMETRIC_STRESS' Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2) S_M = 0.5d0 * Y_M / ( 1 + P_R) C_1 = Y_M /(1 + P_R)/(1 - P_R - P_R) C_2 = C_1 * P_R ; C_3 = C_1 * (1 - P_R) E (:, 3) = 0.d0 ; E (3, :) = 0.d0 E (1, 1) = C_3 ; E (2, 1) = C_2 ; E (4, 1) = C_2 E (1, 2) = C_2 ; E (2, 2) = C_3 ; E (4, 2) = C_2 E (1, 4) = C_2 ; E (2, 4) = C_2 ; E (4, 4) = C_3 E (3, 3) = S_M END SUBROUTINE E_AXISYMMETRIC_STRESS SUBROUTINE E_SOLID_STRESS (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! SOLID_STRESS CONSTITUTIVE MATRIX DEFINITION ! STRESS COMPONENT ORDER: XX, YY, XY, ZZ, XZ, YZ, SO N_R_B = 6 ! PROPERTY ORDER: 1-YOUNG'S MODULUS, 2-POISSON'S RATIO ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B Use Sys_Properties_Data ! for function GET_REAL_LP IMPLICIT NONE REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) REAL(DP) :: C_1, C_2, P_R, S_M, Y_M ! E = CONSTITUTIVE MATRIX ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! P_R = POISSON'S RATIO ! S_M = SHEAR MODULUS ! Y_M = YOUNG'S MODULUS OF ELASTICITY IF ( N_R_B /= 6 ) STOP 'INVALID b_rows IN E_SOLID_STRESS' IF ( EL_REAL < 2 ) STOP 'KEYWORD el_real < 2 IN E_SOLID_STRESS' Y_M = GET_REAL_LP (1) ; P_R = GET_REAL_LP (2) C_1 = Y_M * (1 - P_R) / (1 + P_R) / (1 - P_R - P_R) C_2 = C_1 * P_R / (1 - P_R) ; S_M = 0.5d0 * Y_M / (1 + P_R) E = 0.d0 ! Mostly E (1, 1) = C_1 ; E (2, 1) = C_2 ; E (4, 1) = C_2 ! normals E (1, 2) = C_2 ; E (2, 2) = C_1 ; E (4, 2) = C_2 ! normals E (1, 4) = C_2 ; E (2, 4) = C_2 ; E (4, 4) = C_1 ! normals E (3, 3) = S_M ; E (5, 5) = S_M ; E (6, 6) = S_M ! shears END SUBROUTINE E_SOLID_STRESS SUBROUTINE ELASTICITY_E_MATRIX (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELASTICITY CONSTITUTIVE MATRIX DEFINITION ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_SPACE, ! AXISYMMETRIC, PLANE_STRAIN IMPLICIT NONE REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! CONSTITUTIVE ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! E = SYMMETRIC CONSTITUTIVE MATRIX IF ( ANISOTROPIC ) THEN STOP 'ELASTICITY_E_MATRIX: ANISOTROPIC ELASTICITY E UNAVAILABLE' ! SELECT CASE (N_SPACE) ! CASE (1) ! 1-D or cylinder ! IF ( AXISYMMETRIC ) THEN ! Components RR & TT ! CALL ANISOTROPIC_CYL_STRESS (E) ! ELSE ! Component XX ! CALL E_BAR_STRESS (E) ; END IF ! CASE (2) ! 2-D or axisymmetric ! IF ( AXISYMMETRIC ) THEN ! Components RR ZZ RZ TT ! CALL ANISOTROPIC_AXI_STRESS (E) ! ELSE IF ( PLANE_STRAIN ) THEN ! Components XX YY XY ! CALL ANISOTROPIC_2D_STRAIN (E) ! ELSE ! plane stress Components XX YY XY ! CALL ANISOTROPIC_2D_STRESS (E) ; END IF ! CASE (3) ! 3-D ! CALL ANISOTROPIC_3D_STRESS (E) ! XX YY XY ZZ ZX ZY ! END SELECT ELSE ! Isotropic SELECT CASE (N_SPACE) CASE (1) ! 1-D or cylinder IF ( AXISYMMETRIC ) THEN ! Components RR & TT ! CALL E_CYLINDER_STRESS (E) ELSE ! Component XX ! CALL E_BAR_STRESS (E) END IF CASE (2) ! 2-D or axisymmetric IF ( AXISYMMETRIC ) THEN ! Components RR ZZ RZ TT CALL E_AXISYMMETRIC_STRESS (E) ELSE IF ( PLANE_STRAIN ) THEN ! Components XX YY XY CALL E_PLANE_STRAIN (E) ELSE ! plane stress Components XX YY XY CALL E_PLANE_STRESS (E) ; END IF CASE (3) ! 3-D CALL E_SOLID_STRESS (E) ! XX YY XY ZZ ZX ZY END SELECT END IF ! anisotropic material END SUBROUTINE ELASTICITY_E_MATRIX SUBROUTINE ELASTICITY_B_MATRIX (DGH, XYZ, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERAL ELASTICITY STRAIN-DISPLACEMENT RELATIONS (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, N_R_B, N_G_DOF, N_SPACE, ! AXISYMMETRIC Use Elem_Type_Data ! for LT_FREE, LT_N, H (LT_N) IMPLICIT NONE REAL(DP), INTENT(IN) :: DGH (N_SPACE, LT_N) ! Gradients REAL(DP), INTENT(IN) :: XYZ (N_SPACE) ! Coordinates REAL(DP), INTENT(OUT) :: B (N_R_B, LT_N * N_G_DOF) ! B = STRAIN-DISPLACEMENT MATRIX (RETURNED) ! DGH = GLOBAL DERIVATIVES OF H ! H = ELEMENT INTERPOLATION FUNCTIONS ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! N_R_B = NUMBER OF STRAINS (ROWS IN B) ! N_SPACE = DIMENSION OF SPACE ! XYZ = CURRENT POINT COORDINATES, RADIUS = XYZ (1) SELECT CASE (N_SPACE) CASE (1) ! 1-D or cylinder IF ( AXISYMMETRIC ) THEN ! Components RR & TT CALL ELASTIC_B_CYLINDER (DGH, XYZ (1), B) ELSE ! Component XX CALL ELASTIC_B_BAR (DGH, B) ; END IF CASE (2) ! 2-D or axisymmetric IF ( AXISYMMETRIC ) THEN ! Components RR ZZ RZ TT CALL ELASTIC_B_AXISYMMETRIC (DGH, XYZ (1), B) ELSE ! plane stress or strain. Components XX YY XY CALL ELASTIC_B_PLANAR (DGH, B); END IF CASE (3) ! 3-D CALL ELASTIC_B_SOLID (DGH, B) ! XX YY XY ZZ ZX ZY END SELECT END SUBROUTINE ELASTICITY_B_MATRIX