From ba1d149ebd92fa91fcc7a39d3c87487478ae4796 Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 6 Oct 2017 17:03:09 -0400 Subject: [PATCH 01/10] Solution to Homework 3 --- .../Output_Files/Abaqus_umat_holeplate_3d.out | 146 ------------------ .../Abaqus_Umat_Linear_elastic_2d.in | 46 +++--- EN234_FEA/src/main.f90 | 43 +++--- EN234_FEA/user_codes/src/abaqus_uel_2D.for | 97 +++++++++++- .../src/abaqus_uel_linelastic_3d.for | 4 +- 5 files changed, 140 insertions(+), 196 deletions(-) delete mode 100644 EN234_FEA/Output_Files/Abaqus_umat_holeplate_3d.out diff --git a/EN234_FEA/Output_Files/Abaqus_umat_holeplate_3d.out b/EN234_FEA/Output_Files/Abaqus_umat_holeplate_3d.out deleted file mode 100644 index 574309e..0000000 --- a/EN234_FEA/Output_Files/Abaqus_umat_holeplate_3d.out +++ /dev/null @@ -1,146 +0,0 @@ - - - The input file has been read successfully - Number of elements: 3075 - Number of nodes: 4104 - Number of constraints: 0 - Total no. degrees of freedom: 12312 - - - Direct Solver Profiler has completed profile computation - Total number of equations 12312 - Maximum half bandwidth: 596 - Mean half bandwidth: 375 - RMS half bandwidth: 390 - Size of stiffness matrix: 4613967 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3562E+04 - Min diagonal in reduced matrix: 0.1000E+01 - Ratio: 0.3562E+04 - Maximum no. diagonal digits lost: 1 - - - Newton Raphson iteration 1 - Correction norm 0.390615D+00 - Generalized Force norm 0.568430D+01 - Out of balance force norm 0.240174D-01 - Ratio 0.395353D-02 - Tolerance 0.100000D-04 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3568E+04 - Min diagonal in reduced matrix: 0.1000E+01 - Ratio: 0.3568E+04 - Maximum no. diagonal digits lost: 1 - - - Newton Raphson iteration 2 - Correction norm 0.124137D-03 - Generalized Force norm 0.568438D+01 - Out of balance force norm 0.601132D-04 - Ratio 0.105749D-04 - Tolerance 0.100000D-04 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3568E+04 - Min diagonal in reduced matrix: 0.1000E+01 - Ratio: 0.3568E+04 - Maximum no. diagonal digits lost: 1 - - - Newton Raphson iteration 3 - Correction norm 0.273510D-06 - Generalized Force norm 0.568438D+01 - Out of balance force norm 0.914709D-07 - Ratio 0.160916D-07 - Tolerance 0.100000D-04 - - - Step number 1 has completed successfully - Steps remaining: 2 - Elapsed time: 1.0000 - Time remaining: 9.0000 - Time step has been adjusted to: 1.0000 - Max time step: 1.0000 - Min time step: 0.10000E-02 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3471E-02 - Min diagonal in reduced matrix: 0.2005E-03 - Ratio: 0.1731E+02 - Maximum no. diagonal digits lost: 0 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3574E+04 - Min diagonal in reduced matrix: 0.1000E+01 - Ratio: 0.3574E+04 - Maximum no. diagonal digits lost: 1 - - - Newton Raphson iteration 1 - Correction norm 0.322128D-03 - Generalized Force norm 0.113544D+02 - Out of balance force norm 0.109670D-03 - Ratio 0.965856D-05 - Tolerance 0.100000D-04 - - - Step number 2 has completed successfully - Steps remaining: 1 - Elapsed time: 2.0000 - Time remaining: 8.0000 - Time step has been adjusted to: 1.0000 - Max time step: 1.0000 - Min time step: 0.10000E-02 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3471E-02 - Min diagonal in reduced matrix: 0.2005E-03 - Ratio: 0.1731E+02 - Maximum no. diagonal digits lost: 0 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3580E+04 - Min diagonal in reduced matrix: 0.1000E+01 - Ratio: 0.3580E+04 - Maximum no. diagonal digits lost: 1 - - - Newton Raphson iteration 1 - Correction norm 0.320594D-03 - Generalized Force norm 0.170101D+02 - Out of balance force norm 0.122905D-03 - Ratio 0.722529D-05 - Tolerance 0.100000D-04 - - - Step number 3 has completed successfully - Steps remaining: 0 - Elapsed time: 3.0000 - Time remaining: 7.0000 - Time step has been adjusted to: 1.0000 - Max time step: 1.0000 - Min time step: 0.10000E-02 - - - Direct Solver has completed LU decomposition - Conditioning information: - Max diagonal in reduced matrix: 0.3471E-02 - Min diagonal in reduced matrix: 0.2005E-03 - Ratio: 0.1731E+02 - Maximum no. diagonal digits lost: 0 diff --git a/EN234_FEA/input_files/Abaqus_Umat_Linear_elastic_2d.in b/EN234_FEA/input_files/Abaqus_Umat_Linear_elastic_2d.in index b6c9176..386c0bb 100644 --- a/EN234_FEA/input_files/Abaqus_Umat_Linear_elastic_2d.in +++ b/EN234_FEA/input_files/Abaqus_Umat_Linear_elastic_2d.in @@ -19,24 +19,24 @@ % Enter x,y,z coords of nodes. The node number is optional, and is ignored in the code. COORDINATES % Coords for 8 noded element - 1, 0.d0, 0.d0 - 2, 1.d0, 0.d0 - 3, 2.d0, 0.d0 - 4, 2.d0, 1.d0 - 5, 2.d0, 2.d0 - 6, 1.d0, 2.d0 - 7, 0.d0, 2.d0 - 8, 0.d0, 1.d0 +% 1, 0.d0, 0.d0 +% 2, 1.d0, 0.d0 +% 3, 2.d0, 0.d0 +% 4, 2.d0, 1.d0 +% 5, 2.d0, 2.d0 +% 6, 1.d0, 2.d0 +% 7, 0.d0, 2.d0 +% 8, 0.d0, 1.d0 % Coords for two 6 noded triangles -% 1, 0.d0, 0.d0 -% 2, 1.d0, 0.d0 -% 3, 2.d0, 0.d0 -% 4, 0.d0, 1.d0 -% 5, 1.d0, 1.d0 -% 6, 2.d0, 1.d0 -% 7, 0.d0, 2.d0 -% 8, 1.d0, 2.d0 -% 9, 2.d0, 2.d0 + 1, 0.d0, 0.d0 + 2, 1.d0, 0.d0 + 3, 2.d0, 0.d0 + 4, 0.d0, 1.d0 + 5, 1.d0, 1.d0 + 6, 2.d0, 1.d0 + 7, 0.d0, 2.d0 + 8, 1.d0, 2.d0 + 9, 2.d0, 2.d0 END COORDINATES END NODES @@ -64,9 +64,9 @@ % The element number (first number in the list) is optional, and is ignored in the code CONNECTIVITY, zone1 % Uncomment lines below for 6 noded triangles, and comment out the 8 noded connectivity -% 1, 1, 3, 9, 2, 6, 5 -% 2, 1, 9, 7, 5, 8, 4 - 1, 1, 3, 5, 7, 2, 4, 6, 8 + 1, 1, 3, 9, 2, 6, 5 + 2, 1, 9, 7, 5, 8, 4 +% 1, 1, 3, 5, 7, 2, 4, 6, 8 END CONNECTIVITY % The PROPERTIES, PARAMETERS, CONNECTIVITY keywords can be repeated here to define more set of elements with different properties @@ -88,10 +88,10 @@ % The NODESET key defines a list of nodes NODESET, left - 1, 7, 8 + 1, 4, 7 END NODESET NODESET, right - 3, 4, 5 + 3, 6, 9 END NODESET NODESET, base 1, 2, 3 @@ -144,7 +144,7 @@ INITIAL TIME STEP, 1.d0 MAX TIME STEP, 1.d0 MIN TIME STEP, 0.001d0 - MAX NUMBER OF STEPS, 3 + MAX NUMBER OF STEPS, 10 STOP TIME, 10.d0 STATE PRINT STEP INTERVAL, 1 USER PRINT STEP INTERVAL, 1 diff --git a/EN234_FEA/src/main.f90 b/EN234_FEA/src/main.f90 index eb5860a..b761ebd 100644 --- a/EN234_FEA/src/main.f90 +++ b/EN234_FEA/src/main.f90 @@ -8,10 +8,10 @@ program en234fea character (len=80) :: VS_root_folder character (len=80) :: Eclipse_root_folder - VS_root_folder = 'H:/Repos/EN234_FEA/EN234_FEA/' ! This should work with Intel Studio on the remote desktop if you follow the instructions for cloning your EN234FEA fork + VS_root_folder = 'C:/Users/DongLi/Source/Repos/EN234_FEA/EN234_FEA/' ! This should work with Intel Studio on the remote desktop if you follow the instructions for cloning your EN234FEA fork Eclipse_root_folder = './' ! This should work with Eclipse - VS_root_folder = 'C:/Users/Bower/Source/Repos/EN234_FEA/EN234_FEA/' + VS_root_folder = 'C:/Users/DongLi/Source/Repos/EN234_FEA/EN234_FEA/' root_directory = VS_root_folder ! @@ -31,24 +31,27 @@ program en234fea ! Demo codes - these provide examples of coding and testing ABAQUS user elements in EN234FEA ! ! Small strain linear elasticity - the UEL is in Abaqus_uel_3d.for - infil = 'input_files/Abaqus_uel_linear_elastic_3d.in' - outfil = 'Output_files/Abaqus_uel_linear_elastic_3d.out' + !infil = 'input_files/Abaqus_uel_linear_elastic_3d.in' + !outfil = 'Output_files/Abaqus_uel_linear_elastic_3d.out' ! Linear elastic plate with a central hole using an ABAQUS UEL -! infil = 'input_files/Abaqus_uel_holeplate_3d.in' -! outfil = 'Output_files/Abaqus_uel_holeplate_3d.out' + !infil = 'input_files/Abaqus_uel_holeplate_3d.in' + !outfil = 'Output_files/Abaqus_uel_holeplate_3d.out' ! Simple 1 element demonstration of an ABAQUS VUEL ! The source code for the user element is in abaqus_vuel.for -! infil = 'input_files/Abaqus_vuel_linear_elastic_3d.in' -! outfil = 'Output_files/Abaqus_vuel_linear_elastic_3d.out' + !infil = 'input_files/Abaqus_vuel_linear_elastic_3d.in' + !outfil = 'Output_files/Abaqus_vuel_linear_elastic_3d.out' ! ! Runs an explicit dynamic simulation of a 3D plate with a central hole with and ABAQUS VUEL ! This simulation will take a few minutes to run (running in release mode will speed it up) ! -! infil = 'input_files/Abaqus_uel_holeplate_3d.in' -! outfil = 'output_files/Abaqus_uel_holeplate_3d.out' + !infil = 'input_files/Abaqus_uel_holeplate_3d.in' + !outfil = 'output_files/Abaqus_uel_holeplate_3d.out' + + !infil = 'input_files/Abaqus_uel_test.in' + !outfil = 'output_files/Abaqus_uel_test.out' ! Tests an ABAQUS format UMAT subroutine (in abaqus_umat_elastic.for) with two 8 noded quadrilateral elements ! infil = 'input_files/Abaqus_umat_linear_elastic_3d.in' @@ -66,23 +69,23 @@ program en234fea ! Homework 3: develop and test an ABAQUS user element implementing 2D linear elasticity with full integration ! Simple test of a 2D plane element -! infil = 'input_files/Abaqus_uel_linear_elastic_2d.in' -! outfil = 'Output_files/Abaqus_uel_linear_elastic_2d.out' + !infil = 'input_files/Abaqus_uel_linear_elastic_2d.in' + !outfil = 'Output_files/Abaqus_uel_linear_elastic_2d.out' ! Solve hole-in-a-plate problem with 4 noded quadrilateral elements -! infil = 'input_files/Abaqus_uel_holeplate_2d_quad4.in' -! outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad4.out' + !infil = 'input_files/Abaqus_uel_holeplate_2d_quad4.in' + !outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad4.out' ! Solve hole-in-a-plate problem with 8 noded quads -! infil = 'input_files/Abaqus_uel_holeplate_2d_quad8.in' -! outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad8.out' + !infil = 'input_files/Abaqus_uel_holeplate_2d_quad8.in' + !outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad8.out' ! Solve hole-in-a-plate problem with 3 noded triangles -! infil = 'input_files/Abaqus_uel_holeplate_2d_tri3.in' -! outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri3.out' + !infil = 'input_files/Abaqus_uel_holeplate_2d_tri3.in' + !outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri3.out' -! infil = 'input_files/Abaqus_uel_holeplate_2d_tri6.in' -! outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri6.out' + infil = 'input_files/Abaqus_uel_holeplate_2d_tri6.in' + outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri6.out' ! HW5 Cantilever beam to test incompatible mode elements diff --git a/EN234_FEA/user_codes/src/abaqus_uel_2D.for b/EN234_FEA/user_codes/src/abaqus_uel_2D.for index 44a0400..908c4f3 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_2D.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_2D.for @@ -11,7 +11,7 @@ ! abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) = defines integration points for 1D line integral !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL_2D(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -112,8 +112,8 @@ double precision :: norm(2) ! Normal to an element face double precision :: dxdxi1(2) ! Derivative of 1D spatial coord wrt normalized areal coord ! - double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12, 2e13, 2e23] - double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12, s13, s23] + double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12] + double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12] double precision :: D(4,4) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) double precision :: B(4,22) ! strain = B*(dof_total) double precision :: ktemp(22,22) ! Temporary stiffness (for incompatible mode elements) @@ -135,15 +135,77 @@ ! PROPS(2) Poisson's ratio - if (NNODE == 3) n_points = 1 ! Linear triangle + if (NNODE == 3) n_points = 3 ! Linear triangle if (NNODE == 4) n_points = 4 ! Linear rectangle if (NNODE == 6) n_points = 4 ! Quadratic triangle if (NNODE == 8) n_points = 9 ! Serendipity rectangle if (NNODE == 9) n_points = 9 ! Quadratic rect ! Write your code for a 2D element below + call abq_UEL_2D_integrationpoints(n_points, NNODE, xi, w) - END SUBROUTINE UEL_2D + if (MLVARX<2*NNODE) then + write(6,*) ' Error in abaqus UEL ' + write(6,*) ' Variable MLVARX must exceed 2*NNODE' + write(6,*) ' MLVARX = ',MLVARX,' NNODE = ',NNODE + stop + endif + + RHS(1:MLVARX,1) = 0.d0 + AMATRX(1:NDOFEL,1:NDOFEL) = 0.d0 + + D = 0.d0 + E = PROPS(1) + xnu = PROPS(2) + d44 = 0.5D0*E/(1+xnu) + d11 = (1.D0-xnu)*E/( (1+xnu)*(1-2.D0*xnu) ) + d12 = xnu*E/( (1+xnu)*(1-2.D0*xnu) ) + D(1:3,1:3) = d12 + D(1,1) = d11 + D(2,2) = d11 + D(3,3) = d11 + D(4,4) = d44 + + ENERGY(1:8) = 0.d0 + + ! -- Loop over integration points + do kint = 1, n_points + call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) + dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) + + B = 0.d0 + B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) + B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) + B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) + B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) + + strain = matmul(B(1:4,1:2*NNODE),U(1:2*NNODE)) + + stress = matmul(D,strain) + RHS(1:2*NNODE,1) = RHS(1:2*NNODE,1) + 1 - matmul(transpose(B(1:4,1:2*NNODE)),stress(1:4))* + 2 w(kint)*determinant + + AMATRX(1:2*NNODE,1:2*NNODE) = AMATRX(1:2*NNODE,1:2*NNODE) + 1 + matmul(transpose(B(1:4,1:2*NNODE)),matmul(D,B(1:4,1:2*NNODE))) + 2 *w(kint)*determinant + + ENERGY(2) = ENERGY(2) + 1 + 0.5D0*dot_product(stress,strain)*w(kint)*determinant ! Store the elastic strain energy + + if (NSVARS>=n_points*4) then ! Store stress at each integration point (if space was allocated to do so) + SVARS(4*kint-3:4*kint) = stress(1:4) + endif + end do + + + PNEWDT = 1.d0 ! This leaves the timestep unchanged (ABAQUS will use its own algorithm to determine DTIME) + + return + + END SUBROUTINE UEL @@ -441,7 +503,32 @@ end if end subroutine abq_UEL_2D_shapefunctions + + subroutine abq_UEL_invert2d(A,A_inverse,determinant) + + double precision, intent(in) :: A(2,2) + double precision, intent(out) :: A_inverse(2,2) + double precision, intent(out) :: determinant + +! Compute inverse and determinant of 2x2 matrix + + determinant = A(1,1)*A(2,2) - A(1,2)*A(2,1) + + IF (determinant==0.d0) THEN + write(6,*) ' Error in subroutine abq_UEL_inver2d' + write(6,*) ' A 2x2 matrix has a zero determinant' + stop + endif + + A_inverse(1,1) = A(2,2) + A_inverse(1,2) = -A(1,2) + A_inverse(2,1) = -A(2,1) + A_inverse(2,2) = A(1,1) + + A_inverse = A_inverse / determinant + + end subroutine abq_UEL_invert2d subroutine abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) diff --git a/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for b/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for index e8a43f5..3a4033d 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for @@ -13,7 +13,7 @@ ! !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL_3D(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -239,7 +239,7 @@ return - END SUBROUTINE UEL + END SUBROUTINE UEL_3D subroutine abq_UEL_3D_integrationpoints(n_points, n_nodes, xi, w) From 9132abafc162d87c15c87a7f6cda42fbc77a561a Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 20 Oct 2017 20:00:18 -0400 Subject: [PATCH 02/10] UEL for incompatible element --- .../input_files/Abaqus_uel_cantilever.in | 2 +- EN234_FEA/src/main.f90 | 8 +- EN234_FEA/user_codes/src/abaqus_uel_2D.for | 118 +++++++++++++----- 3 files changed, 91 insertions(+), 37 deletions(-) diff --git a/EN234_FEA/input_files/Abaqus_uel_cantilever.in b/EN234_FEA/input_files/Abaqus_uel_cantilever.in index c52c328..e041aa4 100644 --- a/EN234_FEA/input_files/Abaqus_uel_cantilever.in +++ b/EN234_FEA/input_files/Abaqus_uel_cantilever.in @@ -247,7 +247,7 @@ INITIAL TIME STEP, 1.d0 MAX TIME STEP, 1.d0 MIN TIME STEP, 0.001d0 - MAX NUMBER OF STEPS, 3 + MAX NUMBER OF STEPS, 10 STOP TIME, 10.d0 STATE PRINT STEP INTERVAL, 1 USER PRINT STEP INTERVAL, 1 diff --git a/EN234_FEA/src/main.f90 b/EN234_FEA/src/main.f90 index b761ebd..11aee89 100644 --- a/EN234_FEA/src/main.f90 +++ b/EN234_FEA/src/main.f90 @@ -84,13 +84,13 @@ program en234fea !infil = 'input_files/Abaqus_uel_holeplate_2d_tri3.in' !outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri3.out' - infil = 'input_files/Abaqus_uel_holeplate_2d_tri6.in' - outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri6.out' + !infil = 'input_files/Abaqus_uel_holeplate_2d_tri6.in' + !outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri6.out' ! HW5 Cantilever beam to test incompatible mode elements -! infil = 'input_files/Abaqus_uel_cantilever.in' -! outfil = 'Output_files/Abaqus_uel_cantilever.out' + infil = 'input_files/Abaqus_uel_cantilever.in' + outfil = 'Output_files/Abaqus_uel_cantilever.out' ! HW6 Porous elasticity UMAT diff --git a/EN234_FEA/user_codes/src/abaqus_uel_2D.for b/EN234_FEA/user_codes/src/abaqus_uel_2D.for index 908c4f3..8bd2989 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_2D.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_2D.for @@ -97,6 +97,7 @@ integer :: face_node_list(3) ! List of nodes on an element face ! double precision :: xi(2,9) ! Area integration points + double precision :: xi0(2,1) double precision :: w(9) ! Area integration weights double precision :: N(9) ! 2D shape functions double precision :: dNdxi(9,2) ! 2D shape function derivatives @@ -112,8 +113,8 @@ double precision :: norm(2) ! Normal to an element face double precision :: dxdxi1(2) ! Derivative of 1D spatial coord wrt normalized areal coord ! - double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12] - double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12] + double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12, 2e13, 2e23] + double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12, s13, s23] double precision :: D(4,4) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) double precision :: B(4,22) ! strain = B*(dof_total) double precision :: ktemp(22,22) ! Temporary stiffness (for incompatible mode elements) @@ -125,6 +126,7 @@ double precision :: alpha(4) ! Internal DOF for incompatible mode element double precision :: dxidx(2,2), determinant, det0 ! Jacobian inverse and determinant double precision :: E, xnu, D44, D11, D12 ! Material properties + double precision :: ualpha(22) ! ! Example ABAQUS UEL implementing 2D linear elastic elements @@ -135,7 +137,7 @@ ! PROPS(2) Poisson's ratio - if (NNODE == 3) n_points = 3 ! Linear triangle + if (NNODE == 3) n_points = 1 ! Linear triangle if (NNODE == 4) n_points = 4 ! Linear rectangle if (NNODE == 6) n_points = 4 ! Quadratic triangle if (NNODE == 8) n_points = 9 ! Serendipity rectangle @@ -167,40 +169,92 @@ D(4,4) = d44 ENERGY(1:8) = 0.d0 + + xi0 = 0.d0 + call abq_UEL_2D_shapefunctions(xi0,NNODE,N,dNdxi) + dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,det0) ! -- Loop over integration points + ktemp = 0.d0 do kint = 1, n_points - call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) - dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) - call abq_UEL_invert2d(dxdxi,dxidx,determinant) - dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) + call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) + dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) + + B = 0.d0 + B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) + B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) + B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) + B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) + B(1,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,1) + B(1,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,1) + B(2,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,2) + B(2,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,2) + B(4,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,2) + B(4,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,1) + B(4,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,2) + B(4,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,1) - B = 0.d0 - B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) - B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) - B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) - B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) - - strain = matmul(B(1:4,1:2*NNODE),U(1:2*NNODE)) - - stress = matmul(D,strain) - RHS(1:2*NNODE,1) = RHS(1:2*NNODE,1) - 1 - matmul(transpose(B(1:4,1:2*NNODE)),stress(1:4))* - 2 w(kint)*determinant - - AMATRX(1:2*NNODE,1:2*NNODE) = AMATRX(1:2*NNODE,1:2*NNODE) - 1 + matmul(transpose(B(1:4,1:2*NNODE)),matmul(D,B(1:4,1:2*NNODE))) - 2 *w(kint)*determinant - - ENERGY(2) = ENERGY(2) - 1 + 0.5D0*dot_product(stress,strain)*w(kint)*determinant ! Store the elastic strain energy - - if (NSVARS>=n_points*4) then ! Store stress at each integration point (if space was allocated to do so) - SVARS(4*kint-3:4*kint) = stress(1:4) - endif + ktemp = ktemp + matmul(transpose(B(1:4,1:2*NNODE+4)), + 1 matmul(D,B(1:4,1:2*NNODE+4)))*w(kint)*determinant end do - - + + kuu = ktemp(1:2*NNODE,1:2*NNODE) + kua = ktemp(1:2*NNODE,2*NNODE+1:2*NNODE+4) + kau = ktemp(2*NNODE+1:2*NNODE+4,1:2*NNODE) + kaa = ktemp(2*NNODE+1:2*NNODE+4,2*NNODE+1:2*NNODE+4) + call abq_inverse_LU(kaa,kaainv,4) + + alpha = - matmul(kaainv,matmul(kau(1:4,1:2*NNODE),U(1:2*NNODE))) + ualpha(1:2*NNODE) = U(1:2*NNODE) + ualpha(2*NNODE+1:2*NNODE+4) = alpha + strain = matmul(B(1:4,1:2*NNODE+4),ualpha(1:2*NNODE+4)) + stress = matmul(D,strain) + + ENERGY(2) = ENERGY(2) + 1 + 0.5D0*dot_product(stress,strain) ! Store the elastic strain energy + + rhs_temp = 0.d0 + do kint = 1, n_points + rhs_temp(1:2*NNODE+4) = rhs_temp(1:2*NNODE+4) + 1 - matmul(transpose(B(1:4,1:2*NNODE+4)),stress(1:4))* + 2 w(kint)*determinant + + call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) + dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) + + B = 0.d0 + B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) + B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) + B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) + B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) + B(1,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,1) + B(1,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,1) + B(2,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,2) + B(2,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,2) + B(4,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,2) + B(4,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,1) + B(4,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,2) + B(4,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,1) + + strain = matmul(B(1:4,1:2*NNODE+4),ualpha(1:2*NNODE+4)) + stress = matmul(D,strain) + + if (NSVARS>=n_points*4) then ! Store stress at each integration point (if space was allocated to do so) + SVARS(4*kint-3:4*kint) = stress(1:4) + endif + end do + + RHS(1:2*NNODE,1) = rhs_temp(1:2*NNODE)-matmul(kua(1:2*NNODE,1:4), + 1 matmul(kaainv,rhs_temp(2*NNODE+1:2*NNODE+4))) + + AMATRX(1:2*NNODE,1:2*NNODE) = kuu(1:2*NNODE,1:2*NNODE) - + 1 matmul(kua(1:2*NNODE,1:4),matmul(kaainv,kau(1:4,1:2*NNODE))) + PNEWDT = 1.d0 ! This leaves the timestep unchanged (ABAQUS will use its own algorithm to determine DTIME) return From efc84e7668869f5a073b6d908f9ba7f4b487922f Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 27 Oct 2017 21:10:36 -0400 Subject: [PATCH 03/10] HW6 - nonlinear user material --- .../input_files/Abaqus_umat_porous_elastic.in | 24 ++-- EN234_FEA/src/main.f90 | 8 +- .../user_codes/src/abaqus_umat_elastic.for | 116 ++++++++---------- 3 files changed, 65 insertions(+), 83 deletions(-) diff --git a/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in b/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in index 1960817..b003a2f 100644 --- a/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in +++ b/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in @@ -74,8 +74,8 @@ END HISTORY HISTORY, dload_history - 0.d0, 20.d0 - 10.d0, 20.d0 + 0.d0, 0.d0 + 10.d0, 4.9d0 END HISTORY % The NODESET key defines a list of nodes @@ -101,9 +101,9 @@ % The syntax is node set name, DOF number, and either a value or a history name. % DEGREES OF FREEDOM - 1, 3, VALUE, 0.d0 - side, 2, VALUE, 0.d0 - left, 1, VALUE, 0.d0 + 1, 3, VALUE, 0.d0 + side, 2, VALUE, 0.d0 + left, 1, VALUE, 0.d0 % right, 1, HISTORY, dof_history END DEGREES OF FREEDOM @@ -115,9 +115,9 @@ % element set, face #, HISTORY,history name, nx,(ny),(nz) (time dependent traction or flux to element face in direction (nx,ny,nz)) % element set, face #, NORMAL, history name (applies time dependent pressure normal to element face) % element set, face #, SUBROUTINE, subroutine parameter name - DISTRIBUTED LOADS - end_element, 4, NORMAL,dload_history - END DISTRIBUTED LOADS + DISTRIBUTED LOADS + end_element, 4, NORMAL,dload_history + END DISTRIBUTED LOADS @@ -139,10 +139,10 @@ STATIC STEP - INITIAL TIME STEP, 1.d0 - MAX TIME STEP, 1.d0 + INITIAL TIME STEP, 0.5d0 + MAX TIME STEP, 0.5d0 MIN TIME STEP, 0.001d0 - MAX NUMBER OF STEPS, 3 + MAX NUMBER OF STEPS, 20 STOP TIME, 10.d0 STATE PRINT STEP INTERVAL, 1 USER PRINT STEP INTERVAL, 1 @@ -168,7 +168,7 @@ % In this example the states are the stresses sxx,syy,szz,sxy,sxz,syz PRINT STATE, Output_files\contourplots.dat DEGREES OF FREEDOM - FIELD VARIABLES, S11,S22,S33,S12,S13,S23 + FIELD VARIABLES, S11,S22,S33,S12,S13,S23,E11,E22,E33 DISPLACED MESH DISPLACEMENT SCALE FACTOR, 1.d0 END PRINT STATE diff --git a/EN234_FEA/src/main.f90 b/EN234_FEA/src/main.f90 index 11aee89..7a22df3 100644 --- a/EN234_FEA/src/main.f90 +++ b/EN234_FEA/src/main.f90 @@ -89,13 +89,13 @@ program en234fea ! HW5 Cantilever beam to test incompatible mode elements - infil = 'input_files/Abaqus_uel_cantilever.in' - outfil = 'Output_files/Abaqus_uel_cantilever.out' + !infil = 'input_files/Abaqus_uel_cantilever.in' + !outfil = 'Output_files/Abaqus_uel_cantilever.out' ! HW6 Porous elasticity UMAT -! infil = 'input_files/Abaqus_umat_porous_elastic.in' -! outfil = 'Output_files/Abaqus_umat_porous_elastic.out' + infil = 'input_files/Abaqus_umat_porous_elastic.in' + outfil = 'Output_files/Abaqus_umat_porous_elastic.out' ! HW7 Hyperelastic user element ! infil = 'input_files/Abaqus_uel_hyperelastic.in' diff --git a/EN234_FEA/user_codes/src/abaqus_umat_elastic.for b/EN234_FEA/user_codes/src/abaqus_umat_elastic.for index 458ff10..92afe90 100644 --- a/EN234_FEA/user_codes/src/abaqus_umat_elastic.for +++ b/EN234_FEA/user_codes/src/abaqus_umat_elastic.for @@ -4,14 +4,14 @@ ! SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, - 1 RPL,DDSDDT,DRPLDE,DRPLDT, - 2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, - 3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, - 4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) + 1RPL,DDSDDT,DRPLDE,DRPLDT, + 2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, + 3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, + 4CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC) ! ! INCLUDE 'ABA_PARAM.INC' ! WARNING - the aba_param.inc file declares - Implicit double precision (a-h,o-z) + Implicit double precision (a-h,o-z) ! This means that, by default, any variables with ! first letter between a-h or o-z are double precision. ! The rest are integers. @@ -20,9 +20,9 @@ ! CHARACTER*80 CMNAME DIMENSION STRESS(NTENS),STATEV(NSTATV), - 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), - 2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), - 3 PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) + 1DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), + 2STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1), + 3PROPS(NPROPS),COORDS(3),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3) ! ! Relevant manual sections: @@ -47,14 +47,14 @@ ! This array is passed in as the stress tensor at the beginning ! of the increment and must be updated in this routine to be the ! stress tensor at the end of the increment. If you specified -! initial stresses (�Initial conditions,� Section 19.2.1), this +! initial stresses (�Initial conditions,?Section 19.2.1), this ! array will contain the initial stresses at the start of the ! analysis. The size of this array depends on the value of NTENS ! as defined below. In finite-strain problems the stress tensor ! has already been rotated to account for rigid body motion in ! the increment before UMAT is called, so that only the corotational ! part of the stress integration should be done in UMAT. The -! measure of stress used is �true� (Cauchy) stress. +! measure of stress used is �true?(Cauchy) stress. ! ! NB When used in ABAQUS with hybrid elements the stress array has different dimensions ! and additional variables must be defined. See the ABAQUS manual for details. @@ -64,11 +64,11 @@ ! An array containing the solution-dependent state variables. ! These are passed in as the values at the beginning of the ! increment unless they are updated in user subroutines USDFLD -! (�USDFLD,� Section 25.2.39) or UEXPAN (�UEXPAN,� Section 25.2.20), +! (�USDFLD,?Section 25.2.39) or UEXPAN (�UEXPAN,?Section 25.2.20), ! in which case the updated values are passed in. In all cases ! STATEV must be returned as the values at the end of the increment. ! The size of the array is defined as described in -! �Allocating space� in �User subroutines: overview,� Section 25.1.1. +! �Allocating space?in �User subroutines: overview,?Section 25.1.1. ! ! In finite-strain problems any vector-valued or tensor-valued ! state variables must be rotated to account for rigid body @@ -78,7 +78,7 @@ ! ! SSE, SPD, SCD ! Specific elastic strain energy, plastic dissipation, and -! �creep� dissipation, respectively. These are passed in as +! �creep?dissipation, respectively. These are passed in as ! the values at the start of the increment and should be ! updated to the corresponding specific energy values at ! the end of the increment. They have no effect on the solution, @@ -107,19 +107,19 @@ ! algorithms in ABAQUS/Standard (if automatic time incrementation is chosen). ! For a quasi-static procedure the automatic time stepping that ABAQUS/Standard ! uses, which is based on techniques for integrating standard creep laws -! (see �Quasi-static analysis,� Section 6.2.5), cannot be controlled from within +! (see �Quasi-static analysis,?Section 6.2.5), cannot be controlled from within ! the UMAT subroutine. ! PNEWDT is set to a large value before each call to UMAT. ! If PNEWDT is redefined to be less than 1.0, ABAQUS/Standard must abandon the ! time increment and attempt it again with a smaller time increment. The ! suggested new time increment provided to the automatic time integration -! algorithms is PNEWDT � DTIME, where the PNEWDT used is the minimum value +! algorithms is PNEWDT ?DTIME, where the PNEWDT used is the minimum value ! for all calls to user subroutines that allow redefinition of PNEWDT for this ! iteration. ! If PNEWDT is given a value that is greater than 1.0 for all calls to user ! subroutines for this iteration and the increment converges in this iteration, ! ABAQUS/Standard may increase the time increment. The suggested new time increment -! provided to the automatic time integration algorithms is PNEWDT � DTIME, where +! provided to the automatic time integration algorithms is PNEWDT ?DTIME, where ! the PNEWDT used is the minimum value for all calls to user subroutines for ! this iteration. ! If automatic time incrementation is not selected in the analysis procedure, @@ -134,7 +134,7 @@ ! strains passed into UMAT are the mechanical strains only (that is, the ! thermal strains computed based upon the thermal expansion coefficient have ! been subtracted from the total strains). These strains are available for output -! as the �elastic� strains. +! as the �elastic?strains. ! ! In finite-strain problems the strain components have been rotated to account for ! rigid body motion in the increment before UMAT is called and are approximations @@ -168,7 +168,7 @@ ! Array of increments of predefined field variables. ! ! CMNAME -! User-defined material name, left justified. Some internal material models are given names starting with the �ABQ_� character string. To avoid conflict, you should not use �ABQ_� as the leading string for CMNAME. +! User-defined material name, left justified. Some internal material models are given names starting with the �ABQ_?character string. To avoid conflict, you should not use �ABQ_?as the leading string for CMNAME. ! ! NDI ! Number of direct stress components at this point. @@ -181,8 +181,8 @@ ! ! NSTATV ! Number of solution-dependent state variables that are associated with -! this material type (defined as described in �Allocating space� in �User -! subroutines: overview,� Section 25.1.1). +! this material type (defined as described in �Allocating space?in �User +! subroutines: overview,?Section 25.1.1). ! ! PROPS(NPROPS) ! User-specified array of material constants associated with this user material. @@ -193,7 +193,7 @@ ! COORDS ! An array containing the coordinates of this point. These are the current ! coordinates if geometric nonlinearity is accounted for during the step -! (see �Procedures: overview,� Section 6.1.1); otherwise, the array contains +! (see �Procedures: overview,?Section 6.1.1); otherwise, the array contains ! the original coordinates of the point. ! ! DROT(3,3) @@ -253,57 +253,39 @@ double precision E,xnu integer i,j + double precision A(ntens,ntens),B(ntens,ntens) + double precision eij(ntens),del(ntens),edev(ntens) ddsdde = 0.d0 + A = 0.d0 + B = 0.d0 + del = 0.d0 - E = props(1) + shmod = props(1) xnu = props(2) - -! for debugging, you can use -! write(6,*) ' Hello ' -! Output is then written to the .dat file - - If (ndi==3 .and. nshr==1) then ! Plane strain or axisymmetry - ddsdde(1,1) = 1.d0-xnu - ddsdde(1,2) = xnu - ddsdde(1,3) = xnu - ddsdde(2,1) = xnu - ddsdde(2,2) = 1.d0-xnu - ddsdde(2,3) = xnu - ddsdde(3,1) = xnu - ddsdde(3,2) = xnu - ddsdde(3,3) = 1.d0-xnu - ddsdde(4,4) = 0.5d0*(1.d0-2.d0*xnu) - ddsdde = ddsdde*E/( (1.d0+xnu)*(1.d0-2.d0*xnu) ) - else if (ndi==2 .and. nshr==1) then ! Plane stress - ddsdde(1,1) = 1.d0 - ddsdde(1,2) = xnu - ddsdde(2,1) = xnu - ddsdde(2,2) = 1.d0 - ddsdde(3,3) = 0.5d0*(1.d0-xnu) - ddsdde = ddsdde*E/( (1.d0+xnu*xnu) ) - else ! 3D - ddsdde(1,1) = 1.d0-xnu - ddsdde(1,2) = xnu - ddsdde(1,3) = xnu - ddsdde(2,1) = xnu - ddsdde(2,2) = 1.d0-xnu - ddsdde(2,3) = xnu - ddsdde(3,1) = xnu - ddsdde(3,2) = xnu - ddsdde(3,3) = 1.d0-xnu - ddsdde(4,4) = 0.5d0*(1.d0-2.d0*xnu) - ddsdde(5,5) = ddsdde(4,4) - ddsdde(6,6) = ddsdde(4,4) - ddsdde = ddsdde*E/( (1.d0+xnu)*(1.d0-2.d0*xnu) ) - endif -! -! NOTE: ABAQUS uses engineering shear strains, -! i.e. stran(ndi+1) = 2*e_12, etc... - do i = 1,ntens - do j = 1,ntens - stress(i) = stress(i) + ddsdde(i,j)*dstran(j) + e0 = props(3) + pt = props(4) + + bumod = pt*(1-2*xnu)/(1+xnu)*(1+e0) + tv = exp(-(1+e0)/bumod*(sum(stran(1:3)+dstran(1:3)))) + CapKb = tv*pt/3*(1+e0)/bumod-2*shmod/3 + + do i = 1, NDI + A(i,i) = 2 + do j = 1, NDI + B(i,j) = 1 + end do end do + do i = NDI+1, NDI+NSHR + A(i,i) = 1 + end do + + ddsdde = shmod*A + CapKb*B + + do i = 1,ntens + do j = 1,ntens + stress(i) = stress(i) + ddsdde(i,j)*dstran(j) + end do end do RETURN From 720aff2303970582252eab607a023594e2681ffd Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 27 Oct 2017 21:37:50 -0400 Subject: [PATCH 04/10] HW6 - nonlinear user material --- .../input_files/Abaqus_umat_porous_elastic.in | 10 ++++----- .../user_codes/src/abaqus_umat_elastic.for | 22 ++++++++++++------- 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in b/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in index b003a2f..e3e593a 100644 --- a/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in +++ b/EN234_FEA/input_files/Abaqus_umat_porous_elastic.in @@ -70,7 +70,7 @@ % The HISTORY key defines a time history that can be applied to DOFs or distributed loads HISTORY, dof_history 0.d0, 0.d0 % Each line gives a time value and then a function value - 10.d0, 0.1d0 + 10.d0, 1.0d0 END HISTORY HISTORY, dload_history @@ -104,7 +104,7 @@ 1, 3, VALUE, 0.d0 side, 2, VALUE, 0.d0 left, 1, VALUE, 0.d0 -% right, 1, HISTORY, dof_history + right, 1, HISTORY, dof_history END DEGREES OF FREEDOM @@ -115,9 +115,9 @@ % element set, face #, HISTORY,history name, nx,(ny),(nz) (time dependent traction or flux to element face in direction (nx,ny,nz)) % element set, face #, NORMAL, history name (applies time dependent pressure normal to element face) % element set, face #, SUBROUTINE, subroutine parameter name - DISTRIBUTED LOADS - end_element, 4, NORMAL,dload_history - END DISTRIBUTED LOADS +% DISTRIBUTED LOADS +% end_element, 4, NORMAL,dload_history +% END DISTRIBUTED LOADS diff --git a/EN234_FEA/user_codes/src/abaqus_umat_elastic.for b/EN234_FEA/user_codes/src/abaqus_umat_elastic.for index 92afe90..280c9e8 100644 --- a/EN234_FEA/user_codes/src/abaqus_umat_elastic.for +++ b/EN234_FEA/user_codes/src/abaqus_umat_elastic.for @@ -3,7 +3,7 @@ ! ! - SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, + SUBROUTINE UMAT_1(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1RPL,DDSDDT,DRPLDE,DRPLDT, 2STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME, 3NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT, @@ -277,16 +277,22 @@ end do end do do i = NDI+1, NDI+NSHR - A(i,i) = 1 + A(i,i) = 0.5 end do ddsdde = shmod*A + CapKb*B - do i = 1,ntens - do j = 1,ntens - stress(i) = stress(i) + ddsdde(i,j)*dstran(j) - end do - end do + edev(1:3) = stran(1:3)+dstran(1:3)-sum(stran(1:3)+dstran(1:3))/3 + edev(4:6) = (stran(4:6)+dstran(4:6))/2 + + stress = 2*shmod*edev + stress(4:6) = pt/3*(1-tv) + + !do i = 1,ntens + ! do j = 1,ntens + ! stress(i) = stress(i) + ddsdde(i,j)*dstran(j) + ! end do + !end do RETURN - END subroutine UMAT + END subroutine UMAT_1 From 9473e9b2ce528bfb622abbdedbee2f52f91b9d39 Mon Sep 17 00:00:00 2001 From: dongli314 Date: Sat, 4 Nov 2017 14:25:15 -0400 Subject: [PATCH 05/10] HW7 --- EN234_FEA/src/main.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/EN234_FEA/src/main.f90 b/EN234_FEA/src/main.f90 index 7a22df3..aa5c8bf 100644 --- a/EN234_FEA/src/main.f90 +++ b/EN234_FEA/src/main.f90 @@ -8,10 +8,10 @@ program en234fea character (len=80) :: VS_root_folder character (len=80) :: Eclipse_root_folder - VS_root_folder = 'C:/Users/DongLi/Source/Repos/EN234_FEA/EN234_FEA/' ! This should work with Intel Studio on the remote desktop if you follow the instructions for cloning your EN234FEA fork + VS_root_folder = 'C:/Users/dli33/Source/Repos/EN234_FEA/EN234_FEA/' ! This should work with Intel Studio on the remote desktop if you follow the instructions for cloning your EN234FEA fork Eclipse_root_folder = './' ! This should work with Eclipse - VS_root_folder = 'C:/Users/DongLi/Source/Repos/EN234_FEA/EN234_FEA/' + VS_root_folder = 'C:/Users/dli33/Source/Repos/EN234_FEA/EN234_FEA/' root_directory = VS_root_folder ! @@ -94,12 +94,12 @@ program en234fea ! HW6 Porous elasticity UMAT - infil = 'input_files/Abaqus_umat_porous_elastic.in' - outfil = 'Output_files/Abaqus_umat_porous_elastic.out' + !infil = 'input_files/Abaqus_umat_porous_elastic.in' + !outfil = 'Output_files/Abaqus_umat_porous_elastic.out' ! HW7 Hyperelastic user element -! infil = 'input_files/Abaqus_uel_hyperelastic.in' -! outfil = 'Output_files/Abaqus_uel_hyperelastic.out' + infil = 'input_files/Abaqus_uel_hyperelastic.in' + outfil = 'Output_files/Abaqus_uel_hyperelastic.out' ! Hyperelastic umat ! infil = 'input_files/Abaqus_umat_hyperelastic2.in' From 98d6f271a7f8b1d26f7c7b8d1bce4d0df5e45c66 Mon Sep 17 00:00:00 2001 From: dongli314 Date: Sat, 4 Nov 2017 14:44:38 -0400 Subject: [PATCH 06/10] HW7-UEL for hyperelastic material --- EN234_FEA/user_codes/src/abaqus_uel_2D.for | 4 +- ...tic_3d.for => abaqus_uel_hyperelastic.for} | 374 ++++++++++++++---- .../EN234_Visual_Studio.vfproj | 2 +- 3 files changed, 295 insertions(+), 85 deletions(-) rename EN234_FEA/user_codes/src/{abaqus_uel_linelastic_3d.for => abaqus_uel_hyperelastic.for} (74%) diff --git a/EN234_FEA/user_codes/src/abaqus_uel_2D.for b/EN234_FEA/user_codes/src/abaqus_uel_2D.for index 8bd2989..ff44ca8 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_2D.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_2D.for @@ -11,7 +11,7 @@ ! abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) = defines integration points for 1D line integral !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL_2D(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -259,7 +259,7 @@ return - END SUBROUTINE UEL + END SUBROUTINE UEL_2D diff --git a/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for b/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for similarity index 74% rename from EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for rename to EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for index 3a4033d..3d19aa1 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_linelastic_3d.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for @@ -3,7 +3,7 @@ ! ! This file is compatible with both EN234_FEA and ABAQUS/Standard ! -! The example implements a standard fully integrated 3D linear elastic continuum element +! The example implements a standard fully integrated 3D hyperelastic continuum element ! ! The file also contains the following subrouines: ! abq_UEL_3D_integrationpoints - defines integration ponits for 3D continuum elements @@ -13,7 +13,7 @@ ! !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL_3D(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -115,16 +115,17 @@ double precision :: strain(6) ! Strain vector contains [e11, e22, e33, 2e12, 2e13, 2e23] double precision :: stress(6) ! Stress vector contains [s11, s22, s33, s12, s13, s23] double precision :: D(6,6) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) - double precision :: B(6,60) ! strain = B*(dof_total) + double precision :: B(9,60) ! strain = B*(dof_total) double precision :: dxidx(3,3), determinant ! Jacobian inverse and determinant - double precision :: E, xnu, D44, D11, D12 ! Material properties + + double precision :: mu, Kmod, G11, G22, G33, G44 ! Material properties + double precision :: Gmat(6,6), F(3,3), Finv(3,3), strsq(3,3) + double precision :: Ctemp(3,3), Cinvtemp(3,3), H(6,9) + double precision :: Y(60,60), Sigsq(3,3), HB(6,60) + double precision :: q(9), C(6), Cbar(6), Cinv(6) + double precision :: P(6), Cstar(6), Cbarstar(6), Sigma(6) + double precision :: Jdet, Qsca, Cdet, quan - ! - ! Example ABAQUS UEL implementing 3D linear elastic elements - ! El props are: - - ! PROPS(1) Young's modulus - ! PROPS(2) Poisson's ratio if (NNODE == 4) n_points = 1 ! Linear tet if (NNODE == 10) n_points = 4 ! Quadratic tet @@ -143,23 +144,22 @@ RHS(1:MLVARX,1) = 0.d0 AMATRX(1:NDOFEL,1:NDOFEL) = 0.d0 - D = 0.d0 - E = PROPS(1) - xnu = PROPS(2) - d44 = 0.5D0*E/(1+xnu) - d11 = (1.D0-xnu)*E/( (1+xnu)*(1-2.D0*xnu) ) - d12 = xnu*E/( (1+xnu)*(1-2.D0*xnu) ) - D(1:3,1:3) = d12 - D(1,1) = d11 - D(2,2) = d11 - D(3,3) = d11 - D(4,4) = d44 - D(5,5) = d44 - D(6,6) = d44 + Gmat = 0.d0 + mu = PROPS(1) + Kmod = PROPS(2) + G11 = PROPS(3) + G22 = PROPS(4) + G33 = PROPS(5) + G44 = PROPS(6) + Gmat(1,1) = G11 + Gmat(2,2) = G22 + Gmat(3,3) = G33 + Gmat(4,4) = G44 + Gmat(5,5) = G44 + Gmat(6,6) = G44 ENERGY(1:8) = 0.d0 - ! -- Loop over integration points do kint = 1, n_points call abq_UEL_3D_shapefunctions(xi(1:3,kint),NNODE,N,dNdxi) dxdxi = matmul(coords(1:3,1:NNODE),dNdxi(1:NNODE,1:3)) @@ -170,76 +170,89 @@ B(2,2:3*NNODE-1:3) = dNdx(1:NNODE,2) B(3,3:3*NNODE:3) = dNdx(1:NNODE,3) B(4,1:3*NNODE-2:3) = dNdx(1:NNODE,2) - B(4,2:3*NNODE-1:3) = dNdx(1:NNODE,1) - B(5,1:3*NNODE-2:3) = dNdx(1:NNODE,3) - B(5,3:3*NNODE:3) = dNdx(1:NNODE,1) - B(6,2:3*NNODE-1:3) = dNdx(1:NNODE,3) - B(6,3:3*NNODE:3) = dNdx(1:NNODE,2) - - strain = matmul(B(1:6,1:3*NNODE),U(1:3*NNODE)) - - stress = matmul(D,strain) + B(5,2:3*NNODE-1:3) = dNdx(1:NNODE,1) + B(6,1:3*NNODE-2:3) = dNdx(1:NNODE,3) + B(7,3:3*NNODE:3) = dNdx(1:NNODE,1) + B(8,2:3*NNODE-1:3) = dNdx(1:NNODE,3) + B(9,3:3*NNODE:3) = dNdx(1:NNODE,2) + + F = 0.d0 + do i = 1, 3 + F(i,1:3) = matmul(U(i:3*NNODE-3+i:3),dNdx(1:NNODE,1:3)) + F(i,i) = F(i,i) + 1 + end do + call abq_UEL_invert3d(F,Finv,Jdet) + + call abq_uel_hyper_generateH(F,H) + + Ctemp = matmul(transpose(F),F) + call abq_uel_hyper_sqtovec(Ctemp,C) + Cbar(1:6) = C(1:6)/(Jdet**(2/3)) + call abq_UEL_invert3d(Ctemp,Cinvtemp,Cdet) + call abq_uel_hyper_sqtovec(Cinvtemp,Cinv) + Cstar(1:6) = C(1:6) + Cstar(4:6) = Cstar(4:6)*2 + Cbarstar(1:6) = Cstar(1:6)/(Jdet**(2/3)) + + Qsca = 0.d0 + do i = 1, 3 + Qsca = Qsca + Gmat(i,i)*(Cbarstar(i)-1)*(Cbarstar(i)-1)/4 + end do + do i = 4, 6 + Qsca = Qsca + Gmat(i,i)*Cbarstar(i)*Cbarstar(i)/4 + end do + + quan = 0.d0 + do i = 1, 3 + quan = quan + Gmat(i,i)*Cstar(i)*(Cbarstar(i)-1) + end do + do i = 4, 6 + quan = quan + Gmat(i,i)*Cstar(i)*Cbarstar(i) + end do + do i = 1, 3 + P(i) = (Gmat(i,i)*(Cbarstar(i)-1) - quan*Cinv(i)/3) + 1 /2/(Jdet**(2/3)) + end do + do i = 4, 6 + P(i) = (Gmat(i,i)*Cbarstar(i) - quan*Cinv(i)/3) + 1 /2/(Jdet**(2/3)) + end do + Sigma(1:6) = mu*exp(Qsca)*P(1:6) + 1 +Kmod*Jdet*(Jdet-1)*Cinv(1:6) + call abq_uel_hyper_q(Sigma,F,q) + + call abq_uel_hyper_TanMatrix(mu,Kmod,Qsca,Jdet, + 1 Gmat,Cstar,Cbarstar,Cinv,P,D) + + call abq_uel_hyper_Y(NNODE,Sigma,dNdx,Y) + + call abq_uel_hyper_vectosq(Sigma,Sigsq) + strsq = matmul(F,matmul(Sigsq,transpose(F)))/Jdet + call abq_uel_hyper_sqtovec(strsq,stress) + RHS(1:3*NNODE,1) = RHS(1:3*NNODE,1) - 1 - matmul(transpose(B(1:6,1:3*NNODE)),stress(1:6))* - 2 w(kint)*determinant + 1 - matmul(transpose(B(1:9,1:3*NNODE)),q(1:9))* + 2 w(kint)*determinant + HB = matmul(H,B(1:9,1:3*NNODE)) AMATRX(1:3*NNODE,1:3*NNODE) = AMATRX(1:3*NNODE,1:3*NNODE) - 1 + matmul(transpose(B(1:6,1:3*NNODE)),matmul(D,B(1:6,1:3*NNODE))) - 2 *w(kint)*determinant + 1 + (matmul(transpose(HB(1:6,1:3*NNODE)), + 2 matmul(D,HB(1:6,1:3*NNODE))) + Y(1:3*NNODE,1:3*NNODE)) + 3 *w(kint)*determinant - ENERGY(2) = ENERGY(2) - 1 + 0.5D0*dot_product(stress,strain)*w(kint)*determinant ! Store the elastic strain energy + ENERGY(2) = ENERGY(2) + mu/2*(exp(Qsca)-1) + 1 + K/2*(Jdet-1)*(Jdet-1)*w(kint)*determinant if (NSVARS>=n_points*6) then ! Store stress at each integration point (if space was allocated to do so) SVARS(6*kint-5:6*kint) = stress(1:6) endif end do - PNEWDT = 1.d0 ! This leaves the timestep unchanged (ABAQUS will use its own algorithm to determine DTIME) - ! - ! Apply distributed loads - ! - ! Distributed loads are specified in the input file using the Un option in the input file. - ! n specifies the face number, following the ABAQUS convention - ! - ! - do j = 1,NDLOAD - - call abq_facenodes_3D(NNODE,iabs(JDLTYP(j,1)), - 1 face_node_list,nfacenodes) - - do i = 1,nfacenodes - face_coords(1:3,i) = coords(1:3,face_node_list(i)) - end do - - if (nfacenodes == 3) n_points = 3 - if (nfacenodes == 6) n_points = 4 - if (nfacenodes == 4) n_points = 4 - if (nfacenodes == 8) n_points = 9 - - call abq_UEL_2D_integrationpoints(n_points, nfacenodes, xi2, w) - - do kint = 1,n_points - call abq_UEL_2D_shapefunctions(xi2(1:2,kint), - 1 nfacenodes,N2,dNdxi2) - dxdxi2 = matmul(face_coords(1:3,1:nfacenodes), - 1 dNdxi2(1:nfacenodes,1:2)) - norm(1)=(dxdxi2(2,1)*dxdxi2(3,2))-(dxdxi2(2,2)*dxdxi2(3,1)) - norm(2)=(dxdxi2(1,1)*dxdxi2(3,2))-(dxdxi2(1,2)*dxdxi2(3,1)) - norm(3)=(dxdxi2(1,1)*dxdxi2(2,2))-(dxdxi2(1,2)*dxdxi2(2,1)) - - do i = 1,nfacenodes - ipoin = 3*face_node_list(i)-2 - RHS(ipoin:ipoin+2,1) = RHS(ipoin:ipoin+2,1) - 1 - N2(1:nfacenodes)*adlmag(j,1)*norm(1:3)*w(kint) ! Note determinant is already in normal - end do - end do - end do return - END SUBROUTINE UEL_3D + END SUBROUTINE UEL subroutine abq_UEL_3D_integrationpoints(n_points, n_nodes, xi, w) @@ -569,6 +582,49 @@ end subroutine abq_UEL_3D_shapefunctions + subroutine abq_uel_hyper_generateH(F,H) + + double precision, intent(in) :: F(3,3) + double precision, intent(out) :: H(6,9) + + H = 0.d0 + H(1,1) = F(1,1) + H(4,1) = F(1,2) + H(5,1) = F(1,3) + + H(2,2) = F(2,2) + H(4,2) = F(2,1) + H(6,2) = F(2,3) + + H(3,3) = F(3,3) + H(5,3) = F(3,1) + H(6,3) = F(3,2) + + H(2,4) = F(1,2) + H(4,4) = F(1,1) + H(6,4) = F(1,3) + + H(1,5) = F(2,1) + H(4,5) = F(2,2) + H(5,5) = F(2,3) + + H(3,6) = F(1,3) + H(5,6) = F(1,1) + H(6,6) = F(1,2) + + H(1,7) = F(3,1) + H(4,7) = F(3,2) + H(5,7) = F(3,3) + + H(3,8) = F(2,3) + H(5,8) = F(2,1) + H(6,8) = F(2,2) + + H(2,9) = F(3,2) + H(4,9) = F(3,1) + H(6,9) = F(3,3) + + end subroutine abq_uel_hyper_generateH subroutine abq_UEL_invert3d(A,A_inverse,determinant) @@ -608,6 +664,7 @@ end subroutine abq_UEL_invert3d + subroutine abq_facenodes_3D(nelnodes,face,list,nfacenodes) implicit none @@ -668,5 +725,158 @@ endif end subroutine abq_facenodes_3D - - + + + subroutine abq_uel_hyper_TanMatrix(mu,Kmod,Qsca,Jdet, + 1 Gmat,Cstar,Cbarstar,Cinv,P,D) + + double precision, intent(in) :: mu, Kmod, Jdet, Qsca, P(6) + double precision, intent(in) :: Cstar(6), Cbarstar(6), Cinv(6) + double precision, intent(in) :: Gmat(6,6) + double precision, intent(out) :: D(6,6) + double precision :: Omega(6,6), M1(6,6), M2(6,6), M3(6,6) + double precision :: V1(6), V2(6), Iden(6) + + Iden = 0.d0 + Iden(1:3) = 1 + + Omega(1,1) = Cinv(1)*Cinv(1) + Omega(1,2) = Cinv(4)*Cinv(4) + Omega(1,3) = Cinv(5)*Cinv(5) + Omega(1,4) = Cinv(1)*Cinv(4) + Omega(1,5) = Cinv(1)*Cinv(5) + Omega(1,6) = Cinv(4)*Cinv(5) + + Omega(2,1) = Cinv(4)*Cinv(4) + Omega(2,2) = Cinv(2)*Cinv(2) + Omega(2,3) = Cinv(6)*Cinv(6) + Omega(2,4) = Cinv(4)*Cinv(2) + Omega(2,5) = Cinv(4)*Cinv(6) + Omega(2,6) = Cinv(2)*Cinv(6) + + Omega(3,1) = Cinv(5)*Cinv(5) + Omega(3,2) = Cinv(6)*Cinv(6) + Omega(3,3) = Cinv(3)*Cinv(3) + Omega(3,4) = Cinv(5)*Cinv(6) + Omega(3,5) = Cinv(5)*Cinv(3) + Omega(3,6) = Cinv(6)*Cinv(3) + + Omega(4,1) = Cinv(1)*Cinv(4) + Omega(4,2) = Cinv(4)*Cinv(2) + Omega(4,3) = Cinv(5)*Cinv(6) + Omega(4,4) = (Cinv(1)*Cinv(2)+Cinv(4)*Cinv(4))/2 + Omega(4,5) = (Cinv(1)*Cinv(6)+Cinv(5)*Cinv(4))/2 + Omega(4,6) = (Cinv(4)*Cinv(6)+Cinv(5)*Cinv(2))/2 + + Omega(5,1) = Cinv(1)*Cinv(5) + Omega(5,2) = Cinv(4)*Cinv(6) + Omega(5,3) = Cinv(5)*Cinv(3) + Omega(5,4) = (Cinv(1)*Cinv(6)+Cinv(5)*Cinv(4))/2 + Omega(5,5) = (Cinv(1)*Cinv(3)+Cinv(5)*Cinv(5))/2 + Omega(5,6) = (Cinv(4)*Cinv(3)+Cinv(5)*Cinv(6))/2 + + Omega(6,1) = Cinv(4)*Cinv(5) + Omega(6,2) = Cinv(2)*Cinv(6) + Omega(6,3) = Cinv(6)*Cinv(3) + Omega(6,4) = (Cinv(4)*Cinv(6)+Cinv(5)*Cinv(2))/2 + Omega(6,5) = (Cinv(4)*Cinv(3)+Cinv(5)*Cinv(6))/2 + Omega(6,6) = (Cinv(2)*Cinv(3)+Cinv(6)*Cinv(6))/2 + + Omega = -Omega + + M1 = spread(Cstar,dim=2,ncopies=6)*spread(Cinv,dim=1,ncopies=6) + M2 = spread(Cinv,dim=2,ncopies=6)* + 1 spread(matmul(Gmat,Cstar),dim=1,ncopies=6) + M3 = Gmat - (matmul(Gmat,M1)+M2)/3 + M1 = spread(Cinv,dim=2,ncopies=6)*spread(Cinv,dim=1,ncopies=6) + M3 = M3 + dot_product(Cstar,matmul(Gmat,Cstar))/9*M1 + M3 = M3/(Jdet**(4/3)) + V1 = P - Cinv/3 + M1 = spread(P,dim=2,ncopies=6)*spread(V1,dim=1,ncopies=6) + M3 = M3 + 2*M1 + V1 = matmul(Gmat,Cbarstar-Iden) + M1 = spread(Cinv,dim=2,ncopies=6)*spread(V1,dim=1,ncopies=6) + M2 = dot_product(Cstar,V1)*Omega + M3 = M3 - (M1+M2)/3/(Jdet**(2/3)) + M3 = M3*mu*exp(Qsca) + M1 = spread(Cinv,dim=2,ncopies=6)*spread(Cinv,dim=1,ncopies=6) + D = M3 + Kmod*Jdet*((2*Jdet-1)*M1+2*(Jdet-1)*Omega) + + end subroutine abq_uel_hyper_TanMatrix + + + subroutine abq_uel_hyper_q(Sigma,F,q) + double precision, intent(in) :: Sigma(6), F(3,3) + double precision, intent(out) :: q(9) + double precision :: Sigsq(3,3), qsq(3,3) + + call abq_uel_hyper_vectosq(Sigma,Sigsq) + + qsq = matmul(Sigsq,transpose(F)) + + q(1) = qsq(1,1) + q(2) = qsq(2,2) + q(3) = qsq(3,3) + q(4) = qsq(2,1) + q(5) = qsq(1,2) + q(6) = qsq(3,1) + q(7) = qsq(1,3) + q(8) = qsq(3,2) + q(9) = qsq(2,3) + + end subroutine abq_uel_hyper_q + + + subroutine abq_uel_hyper_Y(nn,Sigma,dNdx,Y) + + integer, intent(in) :: nn + double precision, intent(in) :: Sigma(6), dNdx(20,3) + double precision, intent(out) :: Y(60,60) + integer :: i, j + double precision :: Sigsq(3,3), kij + + call abq_uel_hyper_vectosq(Sigma,Sigsq) + + do i = 1, nn + do j = 1, nn + kij = dot_product(dNdx(i,1:3),matmul(Sigsq,dNdx(j,1:3))) + Y(3*i-2,3*j-2) = kij + Y(3*i-1,3*j-1) = kij + Y(3*i,3*j) = kij + end do + end do + + end subroutine abq_uel_hyper_Y + + + subroutine abq_uel_hyper_sqtovec(A,V) + + double precision, intent(in) :: A(3,3) + double precision, intent(out) :: V(6) + + V(1) = A(1,1) + V(2) = A(2,2) + V(3) = A(3,3) + V(4) = A(1,2) + V(5) = A(1,3) + V(6) = A(2,3) + + end subroutine abq_uel_hyper_sqtovec + + + subroutine abq_uel_hyper_vectosq(V,A) + + double precision, intent(in) :: V(6) + double precision, intent(out) :: A(3,3) + + A(1,1) = V(1) + A(1,2) = V(4) + A(1,3) = V(5) + A(2,1) = V(4) + A(2,2) = V(2) + A(2,3) = V(6) + A(3,1) = V(5) + A(3,2) = V(6) + A(3,3) = V(3) + + end subroutine abq_uel_hyper_vectosq \ No newline at end of file diff --git a/EN234_Visual_Studio/EN234_Visual_Studio/EN234_Visual_Studio.vfproj b/EN234_Visual_Studio/EN234_Visual_Studio/EN234_Visual_Studio.vfproj index b681ad8..ffe530f 100644 --- a/EN234_Visual_Studio/EN234_Visual_Studio/EN234_Visual_Studio.vfproj +++ b/EN234_Visual_Studio/EN234_Visual_Studio/EN234_Visual_Studio.vfproj @@ -111,7 +111,7 @@ - + From f02084079faa7b6cb430073c76d05073e887d5b3 Mon Sep 17 00:00:00 2001 From: dongli314 Date: Fri, 10 Nov 2017 17:12:44 -0500 Subject: [PATCH 07/10] HW8-phasefield --- .../input_files/Abaqus_uel_phasefield_1el.in | 2 +- .../Abaqus_uel_phasefield_coarse.in | 4 +- EN234_FEA/src/main.f90 | 16 +- EN234_FEA/user_codes/src/abaqus_uel_2D.for | 238 +++++++++--------- .../src/abaqus_uel_hyperelastic.for | 4 +- 5 files changed, 125 insertions(+), 139 deletions(-) diff --git a/EN234_FEA/input_files/Abaqus_uel_phasefield_1el.in b/EN234_FEA/input_files/Abaqus_uel_phasefield_1el.in index ccb48ed..974af43 100644 --- a/EN234_FEA/input_files/Abaqus_uel_phasefield_1el.in +++ b/EN234_FEA/input_files/Abaqus_uel_phasefield_1el.in @@ -48,7 +48,7 @@ % Define element properties - the values are passed to user subroutine elstif in the order they are listed here PROPERTIES 100.d0, 0.3d0 % E, nu - 0.0d0, 1.d0, 0.001d0, 1.d0 % Omega, W, kappa, diffusion coefft + 0.5d0, 1.d0, 0.001d0, 1.d0 % Omega, W, kappa, diffusion coefft 0.5d0 % Theta END PROPERTIES % Define element connectivity diff --git a/EN234_FEA/input_files/Abaqus_uel_phasefield_coarse.in b/EN234_FEA/input_files/Abaqus_uel_phasefield_coarse.in index c2bc079..78f7c56 100644 --- a/EN234_FEA/input_files/Abaqus_uel_phasefield_coarse.in +++ b/EN234_FEA/input_files/Abaqus_uel_phasefield_coarse.in @@ -3789,8 +3789,8 @@ % The HISTORY key defines a time history that can be applied to DOFs or distributed loads HISTORY, dof_history - 0.d0, -0.005d0 % Each line gives a time value and then a function value - 10.d0, -0.005d0 + 0.d0, -0.000d0 % Each line gives a time value and then a function value + 10.d0, -0.001d0 END HISTORY HISTORY, dload_history diff --git a/EN234_FEA/src/main.f90 b/EN234_FEA/src/main.f90 index aa5c8bf..91e5e75 100644 --- a/EN234_FEA/src/main.f90 +++ b/EN234_FEA/src/main.f90 @@ -98,8 +98,8 @@ program en234fea !outfil = 'Output_files/Abaqus_umat_porous_elastic.out' ! HW7 Hyperelastic user element - infil = 'input_files/Abaqus_uel_hyperelastic.in' - outfil = 'Output_files/Abaqus_uel_hyperelastic.out' + !infil = 'input_files/Abaqus_uel_hyperelastic.in' + !outfil = 'Output_files/Abaqus_uel_hyperelastic.out' ! Hyperelastic umat ! infil = 'input_files/Abaqus_umat_hyperelastic2.in' @@ -107,14 +107,14 @@ program en234fea ! HW8 - phase field modeling with elasticity ! Single element test -! infil = 'input_files/Abaqus_uel_phasefield_1el.in' -! outfil = 'Output_files/Abaqus_uel_phasefield_1el.out' + !infil = 'input_files/Abaqus_uel_phasefield_1el.in' + !outfil = 'Output_files/Abaqus_uel_phasefield_1el.out' -! infil = 'input_files/Abaqus_uel_phasefield_coarse.in' -! outfil = 'Output_files/Abaqus_uel_phasefield_coarse.out' + !infil = 'input_files/Abaqus_uel_phasefield_coarse.in' + !outfil = 'Output_files/Abaqus_uel_phasefield_coarse.out' -! infil = 'input_files/Abaqus_uel_phasefield_fine.in' -! outfil = 'Output_files/Abaqus_uel_phasefield_fine.out' + infil = 'input_files/Abaqus_uel_phasefield_fine.in' + outfil = 'Output_files/Abaqus_uel_phasefield_fine.out' ! Homework 9 - McCormick model with 1 element diff --git a/EN234_FEA/user_codes/src/abaqus_uel_2D.for b/EN234_FEA/user_codes/src/abaqus_uel_2D.for index ff44ca8..cebf6d9 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_2D.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_2D.for @@ -11,7 +11,7 @@ ! abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) = defines integration points for 1D line integral !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL_2D(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -93,55 +93,41 @@ ! ! ! Local Variables - integer :: i,j,n_points,kint, nfacenodes, ipoin, ksize + integer :: i,j,n_points,kint integer :: face_node_list(3) ! List of nodes on an element face + integer :: nconv ! double precision :: xi(2,9) ! Area integration points - double precision :: xi0(2,1) double precision :: w(9) ! Area integration weights double precision :: N(9) ! 2D shape functions double precision :: dNdxi(9,2) ! 2D shape function derivatives double precision :: dNdx(9,2) ! Spatial derivatives double precision :: dxdxi(2,2) ! Derivative of spatial coords wrt normalized coords - - ! Variables below are for computing integrals over element faces - double precision :: face_coords(2,3) ! Coords of nodes on an element face - double precision :: xi1(6) ! 1D integration points - double precision :: w1(6) ! Integration weights - double precision :: N1(3) ! 1D shape functions - double precision :: dN1dxi(3) ! 1D shape function derivatives - double precision :: norm(2) ! Normal to an element face - double precision :: dxdxi1(2) ! Derivative of 1D spatial coord wrt normalized areal coord + double precision :: Nbar(9) + double precision :: dNbardxi(9,2) + double precision :: dNbardx(9,2) ! - double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12, 2e13, 2e23] - double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12, s13, s23] - double precision :: D(4,4) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) - double precision :: B(4,22) ! strain = B*(dof_total) - double precision :: ktemp(22,22) ! Temporary stiffness (for incompatible mode elements) - double precision :: rhs_temp(22) ! Temporary RHS vector (for incompatible mode elements) - double precision :: kuu(18,18) ! Upper diagonal stiffness - double precision :: kaa(4,4),kaainv(4,4) ! Lower diagonal stiffness - double precision :: kau(4,18) ! Lower quadrant of stiffness - double precision :: kua(18,4) ! Upper quadrant of stiffness - double precision :: alpha(4) ! Internal DOF for incompatible mode element + double precision :: sol(9), dsol(9) ! Strain vector contains [e11, e22, e33, 2e12] + double precision :: strain(4),stress(4) ! Stress vector contains [s11, s22, s33, s12] + double precision :: Dold(4,4),D(9,9) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) + double precision :: B(9,24),B1(9,16),B2(9,8) ! strain = B*(dof_total) double precision :: dxidx(2,2), determinant, det0 ! Jacobian inverse and determinant - double precision :: E, xnu, D44, D11, D12 ! Material properties - double precision :: ualpha(22) - - ! - ! Example ABAQUS UEL implementing 2D linear elastic elements - ! Includes option for incompatible mode elements - ! El props are: + double precision :: E, xnu, d44, d11, d12 ! Material properties + double precision :: iden(4),q(9),c + double precision :: Omega,Wpara,kappa,Dpara,theta - ! PROPS(1) Young's modulus - ! PROPS(2) Poisson's ratio - - if (NNODE == 3) n_points = 1 ! Linear triangle + if (NNODE == 3) n_points = 3 ! Linear triangle if (NNODE == 4) n_points = 4 ! Linear rectangle if (NNODE == 6) n_points = 4 ! Quadratic triangle - if (NNODE == 8) n_points = 9 ! Serendipity rectangle + if (NNODE == 8) n_points = 4 ! Serendipity rectangle if (NNODE == 9) n_points = 9 ! Quadratic rect + + if (NNODE == 3) nconv = 3 ! Linear triangle + if (NNODE == 4) nconv = 4 ! Linear rectangle + if (NNODE == 6) nconv = 3 ! Quadratic triangle + if (NNODE == 8) nconv = 4 ! Serendipity rectangle + if (NNODE == 9) nconv = 4 ! Quadratic rect ! Write your code for a 2D element below call abq_UEL_2D_integrationpoints(n_points, NNODE, xi, w) @@ -156,110 +142,111 @@ RHS(1:MLVARX,1) = 0.d0 AMATRX(1:NDOFEL,1:NDOFEL) = 0.d0 - D = 0.d0 E = PROPS(1) xnu = PROPS(2) - d44 = 0.5D0*E/(1+xnu) - d11 = (1.D0-xnu)*E/( (1+xnu)*(1-2.D0*xnu) ) - d12 = xnu*E/( (1+xnu)*(1-2.D0*xnu) ) - D(1:3,1:3) = d12 - D(1,1) = d11 - D(2,2) = d11 - D(3,3) = d11 - D(4,4) = d44 + Omega = PROPS(3) + Wpara = PROPS(4) + kappa = PROPS(5) + Dpara = PROPS(6) + theta = PROPS(7) ENERGY(1:8) = 0.d0 - xi0 = 0.d0 - call abq_UEL_2D_shapefunctions(xi0,NNODE,N,dNdxi) - dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) - call abq_UEL_invert2d(dxdxi,dxidx,det0) + D = 0.d0 + Dold = 0.d0 + d44 = 0.5D0*E/(1+xnu) + d11 = (1.D0-xnu)*E/( (1+xnu)*(1-2.D0*xnu) ) + d12 = xnu*E/( (1+xnu)*(1-2.D0*xnu) ) + Dold(1:3,1:3) = d12 + Dold(1,1) = d11 + Dold(2,2) = d11 + Dold(3,3) = d11 + Dold(4,4) = d44 + D(1:3,1:3) = Dold([1,2,4],[1,2,4]) + D(1,5) = -Omega*sum(Dold(1,1:3))/3 + D(2,5) = -Omega*sum(Dold(2,1:3))/3 + D(4,1) = -Omega*sum(Dold(1,1:3))/3 + D(4,2) = -Omega*sum(Dold(2,1:3))/3 + D(4,4) = 1.d0 + D(5,5) = 1.d0/DTIME + D(6,8) = -kappa + D(7,9) = -kappa + D(8,6) = theta*Dpara + D(9,7) = theta*Dpara ! -- Loop over integration points - ktemp = 0.d0 do kint = 1, n_points - call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) - dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) - call abq_UEL_invert2d(dxdxi,dxidx,determinant) - dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) - - B = 0.d0 - B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) - B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) - B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) - B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) - B(1,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,1) - B(1,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,1) - B(2,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,2) - B(2,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,2) - B(4,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,2) - B(4,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,1) - B(4,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,2) - B(4,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,1) - ktemp = ktemp + matmul(transpose(B(1:4,1:2*NNODE+4)), - 1 matmul(D,B(1:4,1:2*NNODE+4)))*w(kint)*determinant - end do - - kuu = ktemp(1:2*NNODE,1:2*NNODE) - kua = ktemp(1:2*NNODE,2*NNODE+1:2*NNODE+4) - kau = ktemp(2*NNODE+1:2*NNODE+4,1:2*NNODE) - kaa = ktemp(2*NNODE+1:2*NNODE+4,2*NNODE+1:2*NNODE+4) - call abq_inverse_LU(kaa,kaainv,4) - - alpha = - matmul(kaainv,matmul(kau(1:4,1:2*NNODE),U(1:2*NNODE))) - ualpha(1:2*NNODE) = U(1:2*NNODE) - ualpha(2*NNODE+1:2*NNODE+4) = alpha - strain = matmul(B(1:4,1:2*NNODE+4),ualpha(1:2*NNODE+4)) - stress = matmul(D,strain) - - ENERGY(2) = ENERGY(2) - 1 + 0.5D0*dot_product(stress,strain) ! Store the elastic strain energy - - rhs_temp = 0.d0 - do kint = 1, n_points - rhs_temp(1:2*NNODE+4) = rhs_temp(1:2*NNODE+4) - 1 - matmul(transpose(B(1:4,1:2*NNODE+4)),stress(1:4))* - 2 w(kint)*determinant - - call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) - dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) - call abq_UEL_invert2d(dxdxi,dxidx,determinant) - dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) - - B = 0.d0 - B(1,1:2*NNODE-1:2) = dNdx(1:NNODE,1) - B(2,2:2*NNODE:2) = dNdx(1:NNODE,2) - B(4,1:2*NNODE-1:2) = dNdx(1:NNODE,2) - B(4,2:2*NNODE:2) = dNdx(1:NNODE,1) - B(1,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,1) - B(1,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,1) - B(2,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,2) - B(2,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,2) - B(4,2*NNODE+1) = det0/determinant*xi(1,kint)*dxidx(1,2) - B(4,2*NNODE+2) = det0/determinant*xi(1,kint)*dxidx(1,1) - B(4,2*NNODE+3) = det0/determinant*xi(2,kint)*dxidx(2,2) - B(4,2*NNODE+4) = det0/determinant*xi(2,kint)*dxidx(2,1) - - strain = matmul(B(1:4,1:2*NNODE+4),ualpha(1:2*NNODE+4)) - stress = matmul(D,strain) - - if (NSVARS>=n_points*4) then ! Store stress at each integration point (if space was allocated to do so) - SVARS(4*kint-3:4*kint) = stress(1:4) - endif + call abq_UEL_2D_shapefunctions(xi(1:2,kint),NNODE,N,dNdxi) + dxdxi = matmul(coords(1:2,1:NNODE),dNdxi(1:NNODE,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNdx(1:NNODE,1:2) = matmul(dNdxi(1:NNODE,1:2),dxidx) + + call abq_UEL_2D_shapefunctions(xi(1:2,kint),nconv,Nbar,dNbardxi) + dxdxi = matmul(coords(1:2,1:nconv),dNbardxi(1:nconv,1:2)) + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNbardx(1:nconv,1:2) = matmul(dNbardxi(1:nconv,1:2),dxidx) + + B1 = 0.d0 + B2 = 0.d0 + B1(1,1:4*nconv-3:4) = dNdx(1:nconv,1) + B1(2,2:4*nconv-2:4) = dNdx(1:nconv,2) + B1(3,1:4*nconv-3:4) = dNdx(1:nconv,2) + B1(3,2:4*nconv-2:4) = dNdx(1:nconv,1) + B1(4,3:4*nconv-1:4) = Nbar(1:nconv) + B1(5,4:4*nconv-0:4) = Nbar(1:nconv) + B1(6,3:4*nconv-1:4) = dNbardx(1:nconv,1) + B1(7,3:4*nconv-1:4) = dNbardx(1:nconv,2) + B1(8,4:4*nconv-0:4) = dNbardx(1:nconv,1) + B1(9,4:4*nconv-0:4) = dNbardx(1:nconv,2) + + B2(1,1:2*(NNODE-nconv)-1:2) = dNdx(NNODE-nconv+1:NNODE,1) + B2(2,2:2*(NNODE-nconv)-0:2) = dNdx(NNODE-nconv+1:NNODE,2) + B2(3,1:2*(NNODE-nconv)-1:2) = dNdx(NNODE-nconv+1:NNODE,2) + B2(3,2:2*(NNODE-nconv)-0:2) = dNdx(NNODE-nconv+1:NNODE,1) + B(1:9,1:16) = B1(1:9,1:16) + B(1:9,17:24) = B2(1:9,1:8) + + sol = matmul(B(1:9,1:24),U(1:24)) + dsol = matmul(B(1:9,1:24),DU(1:24,1)) + c = sol(5) + + strain(1:2) = sol(1:2) + strain(3) = 0 + strain(4) = sol(3) + iden = 0.d0 + iden(1:3) = 1.d0 + stress = matmul(Dold(1:4,1:4),(strain(1:4)-iden*Omega*c/3)) + + q(1:2) = stress(1:2) + q(3) = stress(4) + q(4) = sol(4) -2*Wpara*c*(c-1)*(2*c-1) -Omega*sum(stress(1:3))/3 + q(5) = dsol(5)/DTIME + q(6) = -kappa*sol(8) + q(7) = -kappa*sol(9) + q(8) = Dpara*(sol(6)+(theta-1.d0)*dsol(6)) + q(9) = Dpara*(sol(7)+(theta-1.d0)*dsol(7)) + + D(4,5) = -Wpara*(12*c*c-12*c+2) + 1 +Omega*Omega*sum(Dold(1:3,1:3))/9 + + RHS(1:24,1) = RHS(1:24,1) + 1 - matmul(transpose(B(1:9,1:24)),q(1:9))*w(kint)*determinant + + AMATRX(1:24,1:24) = AMATRX(1:24,1:24) + 1 + matmul(transpose(B(1:9,1:24)), + 2 matmul(D,B(1:9,1:24)))*w(kint)*determinant + + if (NSVARS>=n_points*4) then ! Store stress at each integration point (if space was allocated to do so) + SVARS(4*kint-3:4*kint) = stress(1:4) + endif end do - - RHS(1:2*NNODE,1) = rhs_temp(1:2*NNODE)-matmul(kua(1:2*NNODE,1:4), - 1 matmul(kaainv,rhs_temp(2*NNODE+1:2*NNODE+4))) - - AMATRX(1:2*NNODE,1:2*NNODE) = kuu(1:2*NNODE,1:2*NNODE) - - 1 matmul(kua(1:2*NNODE,1:4),matmul(kaainv,kau(1:4,1:2*NNODE))) - + PNEWDT = 1.d0 ! This leaves the timestep unchanged (ABAQUS will use its own algorithm to determine DTIME) return - END SUBROUTINE UEL_2D + END SUBROUTINE UEL @@ -401,7 +388,6 @@ - subroutine abq_UEL_2D_shapefunctions(xi,n_nodes,f,df) implicit none diff --git a/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for b/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for index 3d19aa1..4c5970b 100644 --- a/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for +++ b/EN234_FEA/user_codes/src/abaqus_uel_hyperelastic.for @@ -13,7 +13,7 @@ ! !=========================== ABAQUS format user element subroutine =================== - SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + SUBROUTINE UEL_hyper(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) @@ -252,7 +252,7 @@ return - END SUBROUTINE UEL + END SUBROUTINE UEL_hyper subroutine abq_UEL_3D_integrationpoints(n_points, n_nodes, xi, w) From 3cbcfe344637cb6e2ac084d93d9a72a7bcccd8a8 Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 17 Nov 2017 22:15:54 -0500 Subject: [PATCH 08/10] Add files via upload HW9-McCormick --- Abaqus_vumat_McCormick.for | 255 +++++++++++++++++++++++++++++++++++++ Abaqus_vumat_McCormick.in | 146 +++++++++++++++++++++ abaqus_vumat.f90 | 124 ++++++++++++++++++ main.f90 | 158 +++++++++++++++++++++++ 4 files changed, 683 insertions(+) create mode 100644 Abaqus_vumat_McCormick.for create mode 100644 Abaqus_vumat_McCormick.in create mode 100644 abaqus_vumat.f90 create mode 100644 main.f90 diff --git a/Abaqus_vumat_McCormick.for b/Abaqus_vumat_McCormick.for new file mode 100644 index 0000000..840717c --- /dev/null +++ b/Abaqus_vumat_McCormick.for @@ -0,0 +1,255 @@ +! +! ABAQUS format user material subroutine for explicit dynamics +! +! + + subroutine vumat(nblock, ndir, nshr, nstatev, nfieldv, nprops, + 1 lanneal, stepTime, totalTime, dt, cmname, coordMp, charLength, + 2 props, density, strainInc, relSpinInc, + 3 tempOld, stretchOld, defgradOld, fieldOld, + 4 stressOld, stateOld, enerInternOld, enerInelasOld, + 5 tempNew, stretchNew, defgradNew, fieldNew, + 6 stressNew, stateNew, enerInternNew, enerInelasNew ) +! +! include 'vaba_param.inc' +! + implicit double precision (a-h,o-z) + + dimension props(nprops), density(nblock), coordMp(nblock,*), + 1 charLength(nblock), strainInc(nblock,ndir+nshr), + 2 relSpinInc(nblock,nshr), tempOld(nblock), + 3 stretchOld(nblock,ndir+nshr), + 4 defgradOld(nblock,ndir+nshr+nshr), + 5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), + 6 stateOld(nblock,nstatev), enerInternOld(nblock), + 7 enerInelasOld(nblock), tempNew(nblock), + 8 stretchNew(nblock,ndir+nshr), + 9 defgradNew(nblock,ndir+nshr+nshr), + 1 fieldNew(nblock,nfieldv), + 2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), + 3 enerInternNew(nblock), enerInelasNew(nblock) + + character*80 cmname +! +! Local variables +! + integer k,ntens,nit,maxit + + double precision dedev(ndir+nshr),sdevstar(ndir+nshr) + double precision sestar,skkstar + double precision E,xnu,Y,e0,m,n,edot0 + double precision f,f1,dfde,dfde1 + double precision eplas,deplas,ta,dta +! +! Conventions for storing tensors: +! +! Deformation gradients are provided as components in the global basis +! (ABAQUS also allows the user to define a fixed local coord system defined with *ORIENTATION) +! +! Stresses and stretches are stored as components in a basis that 'rotates with the material' +! These are defined as follows: +! Let e_i be the initial basis vectors (the global ijk directions, or user-specified basis vectors) +! Let R be the rotation tensor, defined through the polar decomposition of deformation gradient F=RU=VR +! Then define rotated basis vectors m_i = R e_i. +! The co-rotational components of a tensor are defined as A_ij = m_i . A . m_j +! The components A_ij can also be interpreted as the global components of the tensor R^T A R +! +! +! The components of symmetric tensors (stress,strainInc and stretch) are stored as vectors in the following order +! For 3D problems (ndir=3,nshr=3) [11,22,33,12,23,13] +! For 2D problems (ndir=3,nshr=1) [11,22,33,12] +! The components of unsymmetric tensors (defgrad and relSpinInc) are stored as vectors in the following order +! For 3D problems (ndir=3,nshr=3) [11,22,33,12,23,31,21,32,13] +! For 2D problems (ndir=3,nshr=1) [11,22,33,12,21] +! +! The stresses are Cauchy (true) stress +! +! nblock No. integration points in data block (data are provided in 'blocks' of integration points) +! EN234FEA always has only 1 int pt in each block +! ndir No. direct tensor components +! nshr No. shear tensor components +! nstatev No. user defined state variables (declared in input file) +! nprops No. material properties +! lanneal Annealing flag - if set to 1, then stresses are zeroed, and state vars should be reset to initial state +! stepTime Elapsed time in this step at end of increment +! totalTime Total elapsed time at end of time increment +! dt time step +! cmname Material name +! coordMp(n,i) ith coord of nth integration point +! charLength(i) Characteristic length of element (in EN234FEA this is (element volume)^1/3) +! props(k) kth material property +! density density +! strainInc(n,i) ith strain increment component for nth int pt. Strain components are in a +! basis that rotates with material, stored in order [de11, de22, de33, de12, de23, de13] +! NOTE THAT SHEAR STRAIN COMPONENTS ARE NOT DOUBLED IN THE VECTOR +! relSpinInc(n,i) ith component of relative spin increment tensor. The global components of relSpinInc are +! dW_ij - dR_ij where dW is the skew part of displacement gradient increment, and +! dR_ij is the rotation increment. +! tempOld Temperature at start of increment +! stretchOld(n,i) ith component of right stretch V in co-rotational basis (equivalent to U in global basis) at start of increment +! defgradOld(n,i) ith component of deformation gradient in fixed global basis at start of increment +! fieldOld(n,i) ith user-defined field variable at nth integration point at start of increment +! stressOld(n,i) ith component of Cauchy stress in co-rotational basis (equivalent to R^T sigma R in fixed basis) +! stateOld(n,i) ith user defined material state variable at start of increment +! enerInternOld(n) specific internal energy at nth integration point, start of increment +! enerInelasOld(n) specific irreversible dissipation at nth integration point, start of increment +! tempNew Temperature at end of increment +! stretchNew(n,i) ith component of right stretch V in co-rotational basis (equivalent to U in global basis) at start of increment +! defgradNew(n,i) ith component of deformation gradient in fixed global basis at end of increment +! fieldNew(n,i) ith user-defined field variable at nth integration point at end of increment +! stressNew(n,i) ith component of Cauchy stress in co-rotational basis (equivale nt to R^T sigma R in fixed basis) +! stateNew(n,i) ith user defined material state variable at end of increment +! enerInternNew(n) specific internal energy at nth integration point, end of increment +! enerInelasNew(n) specific irreversible dissipation at nth integration point, end of increment +! +! Coded for 3D problems; small stretch assumption (works approximately for finite rotations) +! + E = props(1) + xnu = props(2) + Y = props(3) + e0 = props(4) + m = props(5) + edot0 = props(6) + Ssen = props(7) + Hsen = props(8) + td = props(9) + Omega = props(10) + alpha = props(11) + + ntens = ndir+nshr + + do k = 1,nblock + + eplas = stateOld(k,1) + ta = stateOld(k,2) + deplas = stateOld(k,3) + + dedev(1:ndir) = strainInc(k,1:ndir) + 1 - sum(strainInc(k,1:ndir))/3.d0 + dedev(ndir+1:ntens) = strainInc(k,ndir+1:ntens) ! No factor of 2 in dedev + skkstar = sum(stressOld(k,1:ndir)) + 1 + E*sum(strainInc(k,1:ndir))/(1.d0-2.d0*xnu) + sdevstar(1:ndir) = stressOld(k,1:ndir) + 1 - sum(stressOld(k,1:ndir))/3.d0 + sdevstar(ndir+1:ntens) = stressOld(k,ndir+1:ntens) + sdevstar(1:ntens) = sdevstar(1:ntens) + 1 + E*dedev(1:ntens)/(1.d0+xnu) + + sestar = dsqrt(1.5d0)* + 1 dsqrt(dot_product(sdevstar(1:ndir),sdevstar(1:ndir)) + + 2 2.d0*dot_product(sdevstar(ndir+1:ntens), + 3 sdevstar(ndir+1:ntens)) ) + +! Elastic increment (either no stress, or no plastic strain rate) + + if (sestar/Y<1.d-09.or.edot0==0.d0) then + stressNew(k,1:ntens) = sdevstar(1:ntens) + stressNew(k,1:ndir) = stressNew(k,1:ndir) + skkstar/3.d0 + stateNew(k,1) = eplas + stateNew(k,2) = ta + stateNew(k,3) = 0.d0 + cycle + endif + + + err = 1.d0 + maxit = 100 + nit = 1 + tol = 1.d-6*Y + if (deplas==0.d0) deplas = 1.d-09/Y + + do while (err>tol) + +! Calculate the equation that determines the plastic strain increment + call abaqus_vumat_Mcfun(Y,eplas,deplas,e0,m,dt,ta, + 1 Omega,td,alpha,Ssen,Hsen, + 2 edot0,sestar,E,xnu,f,dta) + +! Calculate the derivative of the equation wrt deplas analytically +! call abaqus_vumat_Mcdfde(Y,eplas,deplas,e0,m,dt,ta, +! 1 Omega,td,alpha,Ssen,Hsen, +! 2 edot0,sestar,E,xnu,dfde) + +! Calculate the derivative of the equation wrt deplas numerically + deplas1 = deplas*1.0001 + call abaqus_vumat_Mcfun(Y,eplas,deplas1,e0,m,dt,ta, + 1 Omega,td,alpha,Ssen,Hsen, + 2 edot0,sestar,E,xnu,f1,dta1) + dfde1 = (f1 - f)/(0.0001*deplas) + + deplas_new = deplas - f/dfde1 + if (deplas_new<0.d0) then + deplas = deplas/10.d0 + else + deplas = deplas_new + endif + + nit = nit + 1 + err = dabs(f) + if (nit>maxit) then + write(6,*) ' Newton iterations in UMAT failed to converge ' + stop + endif + + end do + + stressNew(k,1:ntens) = + 1 ( 1.d0 - 1.5d0*deplas*E/(1.d0+xnu)/sestar ) + 2 *sdevstar(1:ntens) + stressNew(k,1:ndir) = + 1 stressNew(k,1:ndir) + skkstar/3.d0 + + stateNew(k,1) = eplas + deplas + stateNew(k,2) = ta + dta + stateNew(k,3) = deplas + + end do + + End Subroutine vumat + + + subroutine abaqus_vumat_Mcfun(Y,eplas,deplas,e0,m,dt,ta, + 1 Omega,td,alpha,Ssen,Hsen, + 2 edot0,sestar,E,xnu,f,dta) + double precision, intent(in) :: Y,eplas,deplas,e0,m,dt,ta + double precision, intent(in) :: Omega,td,alpha,Ssen,Hsen + double precision, intent(in) :: edot0,sestar,E,xnu + double precision, intent(out) :: f,dta + double precision :: sig0,c1,c2,q + + + sig0 = Y*(1.d0 + (eplas + deplas)/e0)**m + q = deplas/Omega + dta = (dt/q+(exp(-q)-1)*ta)/(1+exp(-q)/q) + c1 = 1 - exp(-( ( (ta+dta)/td )**alpha )) + c2 = sig0/Ssen + Hsen*c1 + log(deplas/dt/edot0) + f = sestar/Ssen - 1.5d0*E*deplas/(1.d0+xnu)/Ssen - c2 + + end subroutine abaqus_vumat_Mcfun + + + subroutine abaqus_vumat_Mcdfde(Y,eplas,deplas,e0,m,dt,ta, + 1 Omega,td,alpha,Ssen,Hsen, + 2 edot0,sestar,E,xnu,dfde) + double precision, intent(in) :: Y,eplas,deplas,e0,m,dt,ta + double precision, intent(in) :: Omega,td,alpha,Ssen,Hsen + double precision, intent(in) :: edot0,sestar,E,xnu + double precision, intent(out) :: dfde + double precision :: dta,c,q,q1,q2,q3,q4,a,b,da,db,dtde + + q1 = 1.5d0*E/(1.d0+xnu)/Ssen + q2 = Y/Ssen*((1.d0+(eplas+deplas)/e0)**m) + 1 *(1+m*(eplas+deplas)/(eplas+deplas+e0)) + q = deplas/Omega + a = dt/q+(exp(-q)-1)*ta + b = 1+exp(-q)/q + dta = a/b + c = 1 - exp(-((ta+dta)/td)**alpha) + da = -(dt/q/q + exp(-q)*ta)/Omega + db = -exp(-q)*(1+q)/q/q/Omega + dtde = (da*b-a*db)/b/b + q3 = alpha*Hsen*(1-c)/td*(((eplas+deplas)/td)**(alpha-1))*dtde + q4 = 1.0d0/deplas + dfde = -(q1+q2+q3+q4) + + end subroutine abaqus_vumat_Mcdfde \ No newline at end of file diff --git a/Abaqus_vumat_McCormick.in b/Abaqus_vumat_McCormick.in new file mode 100644 index 0000000..07293cb --- /dev/null +++ b/Abaqus_vumat_McCormick.in @@ -0,0 +1,146 @@ +% +% Demonstration input file for simple general purpose FEA code EN234FEA +% A.F. Bower, August 2017 +% HW9 +% +% Single element test of the McCormick dynamic strain ageing plasticity model +% The material is assumed to be coded as an ABAQUS VUMAT +% +% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MESH DEFINITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% + MESH + +% The NODE command defines properties of the nodes. +% The parameters are # of coords, # of DOF, and an optional integer identifier + NODES +% The parameters are # of coords, # of DOF, and an optional integer identifier + PARAMETERS, 3, 3, 1 +% Specify which nodal DOF are displacements. In the example, DOF 1 is the x displacement, 2 is the y displacement, 3 is the z displacement + DISPLACEMENT DOF, 1, 2, 3 +% Enter x,y,z coords of nodes. The node number is optional, and is ignored in the code. + COORDINATES + 1, 0.d0, 0.d0, 0.d0 + 2, 1.d0, 0.d0, 0.d0 + 3, 1.d0, 1.d0, 0.d0 + 4, 0.d0, 1.d0, 0.d0 + 5, 0.d0, 0.d0, 1.d0 + 6, 1.d0, 0.d0, 1.d0 + 7, 1.d0, 1.d0, 1.d0 + 8, 0.d0, 1.d0, 1.d0 + END COORDINATES + END NODES +% +% The MATERIAL command creates a new material. The material properties can be assigned to ABAQUS style continuum elements to test an ABAQUS UMAT or VUMAT + MATERIAL, mccormick_vumat + STATE VARIABLES, 3 % Number of material state variables (if the key is omitted the number of state vars defaults to zero) + PROPERTIES + 70.d03, 0.3d0 % E (MPa), nu; + 70.d0, 0.001d0, 0.3d0 % Y (MPa), e0, m + 1.d-08, 2.23d0, 27.9d0 % edot0 (s^-1), S (MPa), H (MPa) + 0.02d0, 1.8d-05, 0.336d0 % td (s), Omega, alpha + END PROPERTIES + END MATERIAL +% +% The ELEMENT command defines properties of elements + ELEMENTS, INTERNAL +% The TYPE key selects an ABAQUS format continuum element +% The following elements are supported: C3D4, C3D8, C3D10, C3D20 + TYPE, C3D8 + PROPERTIES, mccormick_vumat + DENSITY, 0.005d0 % We are using mass scaling here the physical density is 2.7*10^-9 in SI units + INITIAL STATE VARIABLES + 0.d0, 10.0d0, 0.d0 % The state vars are accumulated plastic strain, ta, increment in plastic strain + END INITIAL STATE VARIABLES + +% Define element connectivity +% The element number (first number in the list) is optional, and is ignored in the code + CONNECTIVITY, zone1 + 1, 1, 2, 3, 4, 5, 6, 7, 8 + END CONNECTIVITY + +% The PROPERTIES, PARAMETERS, CONNECTIVITY keywords can be repeated here to define more set of elements with different properties + + END ELEMENTS + END MESH + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BOUNDARY CONDITIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% The BOUNDARY conditions key starts definition of BCs + BOUNDARY CONDITIONS + +% The HISTORY key defines a time history that can be applied to DOFs or distributed loads + + HISTORY, strainrate_history + 0.d0, 0.02d0 + 0.5d0, 0.02d0 + 0.50000001d0,0.004d0 + 10.d0, 0.004d0 + END HISTORY + +% The NODESET key defines a list of nodes + NODESET, left + 1, 4, 5, 8 + END NODESET + NODESET, right + 2,3,6,7 + END NODESET + +% The ELEMENTSET key defines a list of elements + ELEMENTSET, end_element + 2 + END ELEMENTSET + +% The DEGREE OF FREEDOM key assigns values to nodal DOFs +% The syntax is node set name, DOF number, and either a value or a history name. +% + DEGREES OF FREEDOM + left, 1, VALUE, 0.d0 + right, 1, HISTORY, strainrate_history, RATE % The RATE prescribes velocity instead of displacement + END DEGREES OF FREEDOM + + END BOUNDARY CONDITIONS + + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% The EXPLICIT DYNAMIC STEP key initializes a dynamic load step + + EXPLICIT DYNAMIC STEP + +% The TIME STEP key defines values of parameters controlling time stepping. +% The parameters are time step, no. steps, no. steps between state prints, no. steps between user prints +% All parameters are required + + TIME STEP, 0.00005 + NUMBER OF STEPS, 20000 + STATE PRINT STEPS, 500 + USER PRINT STEPS, 50 + + +% This prints the DOF values and projected nodal state for all solid elements to a tecplot readable file +% Nodal variables are printed as +% X, Y, (Z), Ux, Uy, (Uz), Projected states. +% Standard continuum elements allow you to print stresses (Sij), strains (Eij), and any material state variables Un if the material has state variables + PRINT STATE, Output_files\contourplots.dat + DEGREES OF FREEDOM + FIELD VARIABLES, S11,S22,S33,S12,S13,S23 + DISPLACED MESH + DISPLACEMENT SCALE FACTOR, 1.d0 + END PRINT STATE + + USER PRINT FILES + Output_files\McCormick_stressvstrain.dat + END USER PRINT FILES + USER PRINT PARAMETERS + 9 % Specify the HW problem number + END USER PRINT PARAMETERS + + + END EXPLICIT DYNAMIC STEP + + + STOP diff --git a/abaqus_vumat.f90 b/abaqus_vumat.f90 new file mode 100644 index 0000000..0ad698c --- /dev/null +++ b/abaqus_vumat.f90 @@ -0,0 +1,124 @@ +! +! ABAQUS format user material subroutine for explicit dynamics +! +! + + subroutine vumat_original(nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal, & + stepTime, totalTime, dt, cmname, coordMp, charLength, & + props, density, strainInc, relSpinInc, & + tempOld, stretchOld, defgradOld, fieldOld, & + stressOld, stateOld, enerInternOld, enerInelasOld, & + tempNew, stretchNew, defgradNew, fieldNew, & + stressNew, stateNew, enerInternNew, enerInelasNew ) +! +! include 'vaba_param.inc' +! + implicit double precision (a-h,o-z) + + dimension props(nprops), density(nblock), coordMp(nblock,*), & + charLength(nblock), strainInc(nblock,ndir+nshr), & + relSpinInc(nblock,nshr), tempOld(nblock), & + stretchOld(nblock,ndir+nshr), & + defgradOld(nblock,ndir+nshr+nshr), & + fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr), & + stateOld(nblock,nstatev), enerInternOld(nblock), & + enerInelasOld(nblock), tempNew(nblock), & + stretchNew(nblock,ndir+nshr), & + defgradNew(nblock,ndir+nshr+nshr), & + fieldNew(nblock,nfieldv), & + stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), & + enerInternNew(nblock), enerInelasNew(nblock) + + character*80 cmname +! +! Local variables +! + integer iblock + + double precision D(6,6) + double precision E,xnu,d11,d12,d44 +! +! Conventions for storing tensors: +! +! Deformation gradients are provided as components in the global basis +! (ABAQUS also allows the user to define a fixed local coord system defined with *ORIENTATION) +! +! Stresses and stretches are stored as components in a basis that 'rotates with the material' +! These are defined as follows: +! Let e_i be the initial basis vectors (the global ijk directions, or user-specified basis vectors) +! Let R be the rotation tensor, defined through the polar decomposition of deformation gradient F=RU=VR +! Then define rotated basis vectors m_i = R e_i. +! The co-rotational components of a tensor are defined as A_ij = m_i . A . m_j +! The components A_ij can also be interpreted as the global components of the tensor R^T A R +! +! +! The components of symmetric tensors (stress,strainInc and stretch) are stored as vectors in the following order +! For 3D problems (ndir=3,nshr=3) [11,22,33,12,23,13] +! For 2D problems (ndir=3,nshr=1) [11,22,33,12] +! The components of unsymmetric tensors (defgrad and relSpinInc) are stored as vectors in the following order +! For 3D problems (ndir=3,nshr=3) [11,22,33,12,23,31,21,32,13] +! For 2D problems (ndir=3,nshr=1) [11,22,33,12,21] +! +! The stresses are Cauchy (true) stress +! +! nblock No. integration points in data block (data are provided in 'blocks' of integration points) +! EN234FEA always has only 1 int pt in each block +! ndir No. direct tensor components +! nshr No. shear tensor components +! nstatev No. user defined state variables (declared in input file) +! nprops No. material properties +! lanneal Annealing flag - if set to 1, then stresses are zeroed, and state vars should be reset to initial state +! stepTime Elapsed time in this step at end of increment +! totalTime Total elapsed time at end of time increment +! dt time step +! cmname Material name +! coordMp(n,i) ith coord of nth integration point +! charLength(i) Characteristic length of element (in EN234FEA this is (element volume)^1/3) +! props(k) kth material property +! density density +! strainInc(n,i) ith strain increment component for nth int pt. Strain components are in a +! basis that rotates with material, stored in order [de11, de22, de33, de12, de23, de13] +! NOTE THAT SHEAR STRAIN COMPONENTS ARE NOT DOUBLED IN THE VECTOR +! relSpinInc(n,i) ith component of relative spin increment tensor. The global components of relSpinInc are +! dW_ij - dR_ij where dW is the skew part of displacement gradient increment, and +! dR_ij is the rotation increment. +! tempOld Temperature at start of increment +! stretchOld(n,i) ith component of right stretch V in co-rotational basis (equivalent to U in global basis) at start of increment +! defgradOld(n,i) ith component of deformation gradient in fixed global basis at start of increment +! fieldOld(n,i) ith user-defined field variable at nth integration point at start of increment +! stressOld(n,i) ith component of Cauchy stress in co-rotational basis (equivalent to R^T sigma R in fixed basis) +! stateOld(n,i) ith user defined material state variable at start of increment +! enerInternOld(n) specific internal energy at nth integration point, start of increment +! enerInelasOld(n) specific irreversible dissipation at nth integration point, start of increment +! tempNew Temperature at end of increment +! stretchNew(n,i) ith component of right stretch V in co-rotational basis (equivalent to U in global basis) at start of increment +! defgradNew(n,i) ith component of deformation gradient in fixed global basis at end of increment +! fieldNew(n,i) ith user-defined field variable at nth integration point at end of increment +! stressNew(n,i) ith component of Cauchy stress in co-rotational basis (equivale nt to R^T sigma R in fixed basis) +! stateNew(n,i) ith user defined material state variable at end of increment +! enerInternNew(n) specific internal energy at nth integration point, end of increment +! enerInelasNew(n) specific irreversible dissipation at nth integration point, end of increment +! +! Coded for 3D problems; small stretch assumption (works approximately for finite rotations) +! + D = 0.d0 + E = PROPS(1) + xnu = PROPS(2) + d44 = 0.5D0*E/(1+xnu) + d11 = (1.D0-xnu)*E/( (1+xnu)*(1-2.D0*xnu) ) + d12 = xnu*E/( (1+xnu)*(1-2.D0*xnu) ) + D(1:3,1:3) = d12 + D(1,1) = d11 + D(2,2) = d11 + D(3,3) = d11 + D(4,4) = 2.d0*d44 ! Important - shear strain components are not doubled in a VUMAT + D(5,5) = 2.d0*d44 + D(6,6) = 2.d0*d44 + + do iblock = 1,nblock + stressNew(iblock,1:6) = stressOld(iblock,1:6) + matmul(D,strainInc(iblock,1:6)) + end do + + +! +End Subroutine vumat_original diff --git a/main.f90 b/main.f90 new file mode 100644 index 0000000..5e0e013 --- /dev/null +++ b/main.f90 @@ -0,0 +1,158 @@ +program en234fea + use Types + use ParamIO + use Globals + use Controlparameters + implicit none + + character (len=80) :: VS_root_folder + character (len=80) :: Eclipse_root_folder + + VS_root_folder = 'C:/Users/dli33/Source/Repos/EN234_FEA/EN234_FEA/' ! This should work with Intel Studio on the remote desktop if you follow the instructions for cloning your EN234FEA fork + Eclipse_root_folder = './' ! This should work with Eclipse + + VS_root_folder = 'C:/Users/dli33/Source/Repos/EN234_FEA/EN234_FEA/' + root_directory = VS_root_folder + +! +! Homework Assignments 2017 +! You can test the ABAQUS UEL, UMAT and VUMAT codes that you develop for EN234 +! by uncommenting the relevant input files listed below. +! +! You will also need to write some code for each assignment: +! A user material for a static analysis should be in a subroutine called UMAT( ) +! A user material for a dynamic analysis should be in a subroutine called VUMAT( ) +! A user element for a static analysis should be in a subroutine called UEL( ) +! A user element for a dynamic analysis should be in a subroutine called VUEL( ) +! +! ABAQUS has a number of other user subroutines, but those are not duplicated in EN234FEA +! + +! Demo codes - these provide examples of coding and testing ABAQUS user elements in EN234FEA +! +! Small strain linear elasticity - the UEL is in Abaqus_uel_3d.for + !infil = 'input_files/Abaqus_uel_linear_elastic_3d.in' + !outfil = 'Output_files/Abaqus_uel_linear_elastic_3d.out' + +! Linear elastic plate with a central hole using an ABAQUS UEL +! infil = 'input_files/Abaqus_uel_holeplate_3d.in' +! outfil = 'Output_files/Abaqus_uel_holeplate_3d.out' + +! Simple 1 element demonstration of an ABAQUS VUEL +! The source code for the user element is in abaqus_vuel.for +! infil = 'input_files/Abaqus_vuel_linear_elastic_3d.in' +! outfil = 'Output_files/Abaqus_vuel_linear_elastic_3d.out' +! + +! Runs an explicit dynamic simulation of a 3D plate with a central hole with and ABAQUS VUEL +! This simulation will take a few minutes to run (running in release mode will speed it up) +! +! infil = 'input_files/Abaqus_uel_holeplate_3d.in' +! outfil = 'output_files/Abaqus_uel_holeplate_3d.out' + +! Tests an ABAQUS format UMAT subroutine (in abaqus_umat_elastic.for) with two 8 noded quadrilateral elements +! infil = 'input_files/Abaqus_umat_linear_elastic_3d.in' +! outfil = 'Output_files/Abaqus_umat_linear_elastic_3d.out' + +! Tests the UMAT on the hole in a plate problem. +! infil = 'input_files/Abaqus_umat_holeplate_3d.in' +! outfil = 'Output_files/Abaqus_umat_holeplate_3d.out' + +! Tests the VUMAT on a hole in a plate problem +! infil = 'input_files/Abaqus_vumat_linear_elastic_3d.in' +! outfil = 'Output_files/Abaqus_vumat_linear_elastic_3d.out' + + +! Homework 3: develop and test an ABAQUS user element implementing 2D linear elasticity with full integration + +! Simple test of a 2D plane element +! infil = 'input_files/Abaqus_uel_linear_elastic_2d.in' +! outfil = 'Output_files/Abaqus_uel_linear_elastic_2d.out' + +! Solve hole-in-a-plate problem with 4 noded quadrilateral elements +! infil = 'input_files/Abaqus_uel_holeplate_2d_quad4.in' +! outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad4.out' + +! Solve hole-in-a-plate problem with 8 noded quads +! infil = 'input_files/Abaqus_uel_holeplate_2d_quad8.in' +! outfil = 'Output_files/Abaqus_uel_holeplate_2d_quad8.out' + +! Solve hole-in-a-plate problem with 3 noded triangles +! infil = 'input_files/Abaqus_uel_holeplate_2d_tri3.in' +! outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri3.out' + +! infil = 'input_files/Abaqus_uel_holeplate_2d_tri6.in' +! outfil = 'Output_files/Abaqus_uel_holeplate_2d_tri6.out' + +! HW5 Cantilever beam to test incompatible mode elements + +! infil = 'input_files/Abaqus_uel_cantilever.in' +! outfil = 'Output_files/Abaqus_uel_cantilever.out' + +! HW6 Porous elasticity UMAT + +! infil = 'input_files/Abaqus_umat_porous_elastic.in' +! outfil = 'Output_files/Abaqus_umat_porous_elastic.out' + +! HW7 Hyperelastic user element +! infil = 'input_files/Abaqus_uel_hyperelastic.in' +! outfil = 'Output_files/Abaqus_uel_hyperelastic.out' + +! Hyperelastic umat +! infil = 'input_files/Abaqus_umat_hyperelastic2.in' +! outfil = 'Output_files/Abaqus_umat_hyperelastic.out' + +! HW8 - phase field modeling with elasticity +! Single element test +! infil = 'input_files/Abaqus_uel_phasefield_1el.in' +! outfil = 'Output_files/Abaqus_uel_phasefield_1el.out' + +! infil = 'input_files/Abaqus_uel_phasefield_coarse.in' +! outfil = 'Output_files/Abaqus_uel_phasefield_coarse.out' + +! infil = 'input_files/Abaqus_uel_phasefield_fine.in' +! outfil = 'Output_files/Abaqus_uel_phasefield_fine.out' + + +! Homework 9 - McCormick model with 1 element + infil = 'input_files/Abaqus_vumat_McCormick.in' + outfil = 'Output_files/Abaqus_vumat_McCormick.out' + + + +! Homework 10 - Continuum beam element solution to end loaded cantilever beam +! infil = 'input_files/Abaqus_uel_continuum_beam.in' +! outfil = 'Output_files/Abaqus_uel_continuum_beam.out' + + infil = trim(root_directory)//trim(infil) + outfil = trim(root_directory)//trim(outfil) + open (unit = IOR, file = trim(infil), status = 'old', ERR=500) + open (UNIT = IOW, FILE = trim(outfil), STATUS = 'unknown', ERR=500) + + call read_input_file + + if (printinitialmesh) call print_initial_mesh + + if (checkstiffness) call check_stiffness(checkstiffness_elementno) + if (checktangent) call check_tangent(checktangent_materialno) + + if (staticstep) then + call compute_static_step + if (checkstiffness) call check_stiffness(checkstiffness_elementno) + if (checktangent) call check_tangent(checktangent_materialno) + endif + + if (explicitdynamicstep) call explicit_dynamic_step + + write(6,*) ' Program completed successfully ' + + stop + + 500 write(6,*) ' Error opening input or output file ' + + + + + + +end program en234fea From 6cf1322e5b0d6091b5aaa002f7ffedf9376c388a Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 1 Dec 2017 01:25:56 -0500 Subject: [PATCH 09/10] HW10-Continuum beam element abaqus_uel_2D.for is the UEL subroutine for continuum beam element --- Abaqus_uel_continuum_beam.in | 165 ++++++++ abaqus_uel_2D.for | 740 +++++++++++++++++++++++++++++++++++ main.f90 | 8 +- 3 files changed, 909 insertions(+), 4 deletions(-) create mode 100644 Abaqus_uel_continuum_beam.in create mode 100644 abaqus_uel_2D.for diff --git a/Abaqus_uel_continuum_beam.in b/Abaqus_uel_continuum_beam.in new file mode 100644 index 0000000..e959877 --- /dev/null +++ b/Abaqus_uel_continuum_beam.in @@ -0,0 +1,165 @@ +% +% Demonstration input file for simple general purpose FEA code EN234FEA +% A.F. Bower, August 2017 +% HW10 +% +% Input file to test a continuum beam element +% The file solves a cantilever beam 10 units long that is clamped at x1=0 and has a point force at x1=10 +% Element is assumed to be coded as an ABAQUS UEL +% +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% MESH DEFINITION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% + MESH + +% The NODE command defines properties of the nodes. +% The parameters are # of coords, # of DOF, and an optional integer identifier + NODES +% The parameters are # of coords, # of DOF, and an optional integer identifier + PARAMETERS, 2, 3, 1 +% Specify which nodal DOF are displacements. In the example, DOF 1 is the x displacement, 2 is the y displacement, 3 is the z displacement + DISPLACEMENT DOF, 1, 2 +% Enter x,y,z coords of nodes. The node number is optional, and is ignored in the code. + COORDINATES +% Coords + 1, 0., 0. + 2, 1., 0. + 3, 2., 0. + 4, 3., 0. + 5, 4., 0. + 6, 5., 0. + 7, 6., 0. + 8, 7., 0. + 9, 8., 0. + 10, 9., 0. + 11, 10., 0. + END COORDINATES + END NODES +% +% The ELEMENT command defines properties of elements +% The parameters are no. nodes on the element, total no. state variables (4*# integration points), integer identifier + + ELEMENTS, USER + PARAMETERS, 2, 3, U1 +% Define element properties - the values are passed to user subroutine elstif in the order they are listed here +% For the example provided, the params are beam thickness h, width w, Youngs Modulus, Poissons ratio + PROPERTIES + 1.2d0, 1.d0, 100.d0, 0.3d0 + END PROPERTIES +% Define element connectivity +% The element number (first number in the list) is optional, and is ignored in the code + CONNECTIVITY, zone1 + 1, 1, 2 + 2, 2, 3 + 3, 3, 4 + 4, 4, 5 + 5, 5, 6 + 6, 6, 7 + 7, 7, 8 + 8, 8, 9 + 9, 9, 10 + 10, 10, 11 + END CONNECTIVITY + +% The PROPERTIES, PARAMETERS, CONNECTIVITY keywords can be repeated here to define more set of elements with different properties + + END ELEMENTS + END MESH + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% BOUNDARY CONDITIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% The BOUNDARY conditions key starts definition of BCs + BOUNDARY CONDITIONS + +% The HISTORY key defines a time history that can be applied to DOFs or distributed loads + HISTORY, force_history + 0.d0, 0.d0 % Each line gives a time value and then a function value + 1.d0, 0.01d0 + END HISTORY + +% The NODESET key defines a list of nodes + NODESET, left + 1 + END NODESET + NODESET, right + 11 + END NODESET + + +% The ELEMENTSET key defines a list of elements + ELEMENTSET, end_element + 1 + END ELEMENTSET + +% The DEGREE OF FREEDOM key assigns values to nodal DOFs +% The syntax is node set name, DOF number, VALUE/HISTORY/SUBROUTINE, value/history name/subroutine parameter list name. +% + DEGREES OF FREEDOM + left, 1, VALUE, 0.d0 + left, 2, VALUE, 0.d0 + left, 3, VALUE, 0.d0 + END DEGREES OF FREEDOM + + FORCES + right, 2, HISTORY, force_history + END FORCES + + + END BOUNDARY CONDITIONS + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Mesh printing, error checking %%%%%%%%%%%%%%%%%%%% + +% Print the initial mesh to a file named initial_mesh.dat + +% PRINT INITIAL MESH, Output_files\initial_mesh.dat + +% TIME, VALUE, 0.d0 % Use this to specify the initial time +% TIME, INCREMENT, 0.01d0 % Use this to specify a time increment (often needed for check stiffness) + +% The CHECK STIFFNESS key tests the element subroutine to ensure that +% the residual force vector is consistent with the stiffness + CHECK STIFFNESS, U1 + +% STOP + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % The STATIC STEP key initializes a static load step + + STATIC STEP + + + INITIAL TIME STEP, 1.d0 + MAX TIME STEP, 1.d0 + MIN TIME STEP, 0.001d0 + MAX NUMBER OF STEPS, 1 + STOP TIME, 10.d0 + USER PRINT STEP INTERVAL, 1 + + +% The SOLVER key controls the equation solver and Newton-Raphson iterations +% The options are FACTOR for direct solver, CONJUGATE GRADIENT for cg solver +% Factor will work on anything but might be slow for large equation systems. +% Conjugate gradient works well for elasticity problems but (with the diagonal preconditioner used here) is not so good for unsymmetric matrices +% LINEAR for linear equations, NONLINEAR for nonlinear equations +% For nonlinear solver, must specify convergence tolerance and max # iterations +% UNSYMMETRIC (optional - only for unsymmetric stiffness) + + + SOLVER, DIRECT, NONLINEAR, 1.d-05,15 + + USER PRINT FILES + Output_files\beam_displacements.dat + Output_files\beam_forces.dat + END USER PRINT FILES + USER PRINT PARAMETERS + 10 % Specify the HW problem number + END USER PRINT PARAMETERS + + + + + END STATIC STEP + + + STOP diff --git a/abaqus_uel_2D.for b/abaqus_uel_2D.for new file mode 100644 index 0000000..bc9e71a --- /dev/null +++ b/abaqus_uel_2D.for @@ -0,0 +1,740 @@ +! +! ABAQUS format UEL subroutine +! +! This file is compatible with both EN234_FEA and ABAQUS/Standard +! +! The example implements a standard fully integrated 3D linear elastic continuum element +! +! The file also contains the following subrouines: +! abq_UEL_2D_integrationpoints - defines integration points for 2D continuum elements +! abq_UEL_2D_shapefunctions - defines shape functions for 2D continuum elements +! abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) = defines integration points for 1D line integral +!=========================== ABAQUS format user element subroutine =================== + + SUBROUTINE UEL(RHS,AMATRX,SVARS,ENERGY,NDOFEL,NRHS,NSVARS, + 1 PROPS,NPROPS,COORDS,MCRD,NNODE,U,DU,V,A,JTYPE,TIME,DTIME, + 2 KSTEP,KINC,JELEM,PARAMS,NDLOAD,JDLTYP,ADLMAG,PREDEF,NPREDF, + 3 LFLAGS,MLVARX,DDLMAG,MDLOAD,PNEWDT,JPROPS,NJPROP,PERIOD) + ! + INCLUDE 'ABA_PARAM.INC' + ! + ! + DIMENSION RHS(MLVARX,*),AMATRX(NDOFEL,NDOFEL),PROPS(*), + 1 SVARS(*),ENERGY(8),COORDS(MCRD,NNODE),U(NDOFEL), + 2 DU(MLVARX,*),V(NDOFEL),A(NDOFEL),TIME(2),PARAMS(*), + 3 JDLTYP(MDLOAD,*),ADLMAG(MDLOAD,*),DDLMAG(MDLOAD,*), + 4 PREDEF(2,NPREDF,NNODE),LFLAGS(*),JPROPS(*) + + ! + ! Variables that must be computed in this routine + ! RHS(i) Right hand side vector. In EN234_FEA the dimensions are always RHS(MLVARX,1) + ! AMATRX(i,j) Stiffness matrix d RHS(i)/ d DU(j) + ! SVARS(1:NSVARS) Element state variables. Must be updated in this routine + ! ENERGY(1:8) + ! Energy(1) Kinetic Energy + ! Energy(2) Elastic Strain Energy + ! Energy(3) Creep Dissipation + ! Energy(4) Plastic Dissipation + ! Energy(5) Viscous Dissipation + ! Energy(6) Artificial strain energy + ! Energy(7) Electrostatic energy + ! Energy(8) Incremental work done by loads applied to the element + ! PNEWDT Allows user to control ABAQUS time increments. + ! If PNEWDT<1 then time step is abandoned and computation is restarted with + ! a time increment equal to PNEWDT*DTIME + ! If PNEWDT>1 ABAQUS may increase the time increment by a factor PNEWDT + ! + ! Variables provided for information + ! NDOFEL Total # DOF for the element + ! NRHS Dimension variable + ! NSVARS Total # element state variables + ! PROPS(1:NPROPS) User-specified properties of the element + ! NPROPS No. properties + ! JPROPS(1:NJPROPS) Integer valued user specified properties for the element + ! NJPROPS No. integer valued properties + ! COORDS(i,N) ith coordinate of Nth node on element + ! MCRD Maximum of (# coords,minimum of (3,#DOF)) on any node + ! U Vector of DOF at the end of the increment + ! DU Vector of DOF increments + ! V Vector of velocities (defined only for implicit dynamics) + ! A Vector of accelerations (defined only for implicit dynamics) + ! JTYPE Integer identifying element type (the number n in the Un specification in the input file) + ! TIME(1:2) TIME(1) Current value of step time + ! TIME(2) Total time + ! DTIME Time increment + ! KSTEP Current step number (always 1 in EN234_FEA) + ! KINC Increment number + ! JELEM User assigned element number in ABAQUS (internally assigned in EN234_FEA) + ! PARAMS(1:3) Time increment parameters alpha, beta, gamma for implicit dynamics + ! NDLOAD Number of user-defined distributed loads defined for this element + ! JDLTYP(1:NDLOAD) Integers n defining distributed load types defined as Un or (if negative) UnNU in input file + ! ADLMAG(1:NDLOAD) Distributed load magnitudes + ! DDLMAG(1:NDLOAD) Increment in distributed load magnitudes + ! PREDEF(1:2,1:NPREDF,1:NNODE) Predefined fields. + ! PREDEF(1,...) Value of predefined field + ! PREDEF(2,...) Increment in predefined field + ! PREDEF(1:2,1,k) Value of temperature/temperature increment at kth node + ! PREDEF(1:2,2:NPREDF,k) Value of user defined field/field increment at kth node (not used in EN234FEA) + ! NPREDF Number of predefined fields (1 for en234FEA) + ! LFLAGS Control variable + ! LFLAGS(1) Defines procedure type + ! LFLAGS(2) 0 => small displacement analysis 1 => Large displacement (NLGEOM option) + ! LFLAGS(3) 1 => Subroutine must return both RHS and AMATRX (always true in EN234FEA) + ! 2 => Subroutine must return stiffness AMATRX = -dF/du + ! 3 => Subroutine must return daming matrix AMATRX = -dF/dudot + ! 4 => Subroutine must return mass matrix AMATRX = -dF/duddot + ! 5 => Define the RHS only + ! 6 => Define the mass matrix for the initial acceleration calculation + ! 100 => Define perturbation quantities for output + ! LFLAGS(4) 0 => General step 1 => linear perturbation step + ! LFLAGS(5) 0 => current approximation to solution based on Newton correction; 1 => based on extrapolation + ! MLVARX Dimension variable (equal to NDOFEL in EN234FEA) + ! PERIOD Time period of the current step + ! + ! + ! Local Variables + integer :: i,j,n_points,kint, nfacenodes, ipoin, ksize + integer :: face_node_list(3) ! List of nodes on an element face + ! + double precision :: xi(2,9) ! Area integration points + double precision :: w(9) ! Area integration weights + double precision :: N(9) ! 2D shape functions + double precision :: dNdxi(9,2) ! 2D shape function derivatives + double precision :: dNdx(9,2) ! Spatial derivatives + double precision :: dxdxi(2,2), dxdxi1(2,2) ! Derivative of spatial coords wrt normalized coords + ! + double precision :: strain(4) ! Strain vector contains [e11, e22, e33, 2e12, 2e13, 2e23] + double precision :: stress(4) ! Stress vector contains [s11, s22, s33, s12, s13, s23] + double precision :: D(2,2) ! stress = D*(strain) (NOTE FACTOR OF 2 in shear strain) + double precision :: B(4,22) ! strain = B*(dof_total) + double precision :: ktemp(22,22) ! Temporary stiffness (for incompatible mode elements) + double precision :: rhs_temp(22) ! Temporary RHS vector (for incompatible mode elements) + double precision :: alpha(4) ! Internal DOF for incompatible mode element + double precision :: dxidx(2,2), determinant, det0 ! Jacobian inverse and determinant + double precision :: E, xnu, D44, D11, D12 ! Material properties + double precision :: thk, wid, lbar, dv(2), sn1, sn2, cs1, cs2 + double precision :: scoords(2,4), Tmat(8,6), Rmat(2,3) + double precision :: strainhat(2), stresshat(2) + double precision :: RBT(2,6), stavar(3) + + ! + ! Example ABAQUS UEL implementing 2D linear elastic elements + ! Includes option for incompatible mode elements + ! El props are: + + ! PROPS(1) Young's modulus + ! PROPS(2) Poisson's ratio + + + if (NNODE == 2) n_points = 5 ! 1-D element + if (NNODE == 3) n_points = 1 ! Linear triangle + if (NNODE == 4) n_points = 4 ! Linear rectangle + if (NNODE == 6) n_points = 4 ! Quadratic triangle + if (NNODE == 8) n_points = 9 ! Serendipity rectangle + if (NNODE == 9) n_points = 9 ! Quadratic rect + + ! Write your code for a 2D element below + call abq_UEL_1D_integrationpoints(n_points, NNODE, xi, w) + + if (MLVARX<2*NNODE) then + write(6,*) ' Error in abaqus UEL ' + write(6,*) ' Variable MLVARX must exceed 2*NNODE' + write(6,*) ' MLVARX = ',MLVARX,' NNODE = ',NNODE + stop + endif + + RHS(1:MLVARX,1) = 0.d0 + AMATRX(1:NDOFEL,1:NDOFEL) = 0.d0 + stavar = 0.d0 + + D = 0.d0 + thk = PROPS(1) + wid = PROPS(2) + E = PROPS(3) + xnu = PROPS(4) + D(1,1) = E + D(2,2) = E/2/(1.d0+xnu) + + ENERGY(1:8) = 0.d0 + + lbar = norm2(coords(1:2,2)-coords(1:2,1)) + dv(1:2) = (coords(1:2,2)-coords(1:2,1))/lbar + sn1 = dv(2)*cos(U(3)) + dv(1)*sin(U(3)) + cs1 = dv(1)*cos(U(3)) - dv(2)*sin(U(3)) + sn2 = dv(2)*cos(U(6)) + dv(1)*sin(U(6)) + cs2 = dv(1)*cos(U(6)) - dv(2)*sin(U(6)) + scoords(1,1) = coords(1,1) + thk/2*sn1 + scoords(2,1) = coords(2,1) - thk/2*cs1 + scoords(1,2) = coords(1,2) + thk/2*sn2 + scoords(2,2) = coords(2,2) - thk/2*cs2 + scoords(1,3) = coords(1,2) - thk/2*sn2 + scoords(2,3) = coords(2,2) + thk/2*cs2 + scoords(1,4) = coords(1,1) - thk/2*sn1 + scoords(2,4) = coords(2,1) + thk/2*cs1 + + Tmat = 0.d0 + Tmat(1,1) = 1.d0 + Tmat(1,3) = coords(2,1) - scoords(2,1) + Tmat(2,2) = 1.d0 + Tmat(2,3) = -coords(1,1) + scoords(1,1) + Tmat(3,4) = 1.d0 + Tmat(3,6) = coords(2,2) - scoords(2,2) + Tmat(4,5) = 1.d0 + Tmat(4,6) = coords(1,2) - scoords(1,2) + Tmat(5,4) = 1.d0 + Tmat(5,6) = coords(2,2) - scoords(2,3) + Tmat(6,5) = 1.d0 + Tmat(6,6) = coords(1,2) - scoords(1,3) + Tmat(7,1) = 1.d0 + Tmat(7,3) = coords(2,1) - scoords(2,4) + Tmat(8,2) = 1.d0 + Tmat(8,3) = -coords(1,1) + scoords(1,4) + + ! -- Loop over integration points + do kint = 1, n_points + call abq_UEL_2D_shapefunctions(xi(1:2,kint),4,N,dNdxi) + dxdxi = matmul(scoords(1:2,1:4),dNdxi(1:4,1:2)) + dxdxi1(1:2,1) = dxdxi(1:2,1)/norm2(dxdxi(1:2,1)) + dxdxi1(1:2,2) = dxdxi(1:2,2)/norm2(dxdxi(1:2,2)) + + Rmat(1,1) = dxdxi1(1,1)*dxdxi1(1,1) + Rmat(1,2) = dxdxi1(1,2)*dxdxi1(1,2) + Rmat(1,3) = dxdxi1(1,1)*dxdxi1(1,2) + Rmat(2,1) = 2*dxdxi1(1,1)*dxdxi1(2,1) + Rmat(2,2) = 2*dxdxi1(1,2)*dxdxi1(2,2) + Rmat(2,3) = dxdxi1(1,1)*dxdxi1(2,2) + dxdxi1(1,2)*dxdxi1(2,1) + + call abq_UEL_invert2d(dxdxi,dxidx,determinant) + dNdx(1:4,1:2) = matmul(dNdxi(1:4,1:2),dxidx) + + B = 0.d0 + B(1,1:2*4-1:2) = dNdx(1:4,1) + B(2,2:2*4:2) = dNdx(1:4,2) + B(3,1:2*4-1:2) = dNdx(1:4,2) + B(3,2:2*4:2) = dNdx(1:4,1) + + strain(1:3) = matmul(B(1:3,1:8),matmul(Tmat(1:8,1:6),U(1:6))) + strainhat = matmul(Rmat(1:2,1:3),strain(1:3)) + stresshat = matmul(D,strainhat) + + RBT = matmul(Rmat(1:2,1:3),matmul(B(1:3,1:8),Tmat(1:8,1:6))) + RHS(1:3*NNODE,1) = RHS(1:3*NNODE,1) + 1 - wid*lbar*thk/2*matmul(transpose(RBT),stresshat(1:2))*w(kint) + + AMATRX(1:3*NNODE,1:3*NNODE) = AMATRX(1:3*NNODE,1:3*NNODE) + 1 + wid*lbar*thk/2*matmul(transpose(RBT),matmul(D,RBT))*w(kint) + + stavar(1) = stavar(1) + wid*thk/2*stresshat(1)*w(kint) + stavar(2) = stavar(2) + wid*thk/2*stresshat(2)*w(kint) + stavar(3) = stavar(3) + wid*thk*thk/4 + 1 *stresshat(1)*xi(2,kint)*w(kint) + + end do + + SVARS(1:3) = stavar(1:3) + + PNEWDT = 1.d0 ! This leaves the timestep unchanged (ABAQUS will use its own algorithm to determine DTIME) + + return + + END SUBROUTINE UEL + + + + subroutine abq_UEL_2D_integrationpoints(n_points, n_nodes, xi, w) + + implicit none + integer, intent(in) :: n_points + integer, intent(in) :: n_nodes + + double precision, intent(out) :: xi(2,*) + double precision, intent(out) :: w(*) + + integer :: i,j,k,n + + double precision :: cn,w1,w2,w11,w12,w22 + + ! Defines integration points and weights for 2D continuum elements + + if ( n_points==1 ) then + if ( n_nodes==4 .or. n_nodes==9 ) then ! --- 4 or 9 noded quad + xi(1, 1) = 0.D0 + xi(2, 1) = 0.D0 + w(1) = 4.D0 + else if ( n_nodes==3 .or. n_nodes==6 ) then ! --- 3 or 6 noded triangle + xi(1, 1) = 1.D0/3.D0 + xi(2, 1) = 1.D0/3.D0 + w(1) = 1.D0/2.D0 + end if + else if ( n_points==3 ) then + xi(1, 1) = 0.5D0 + xi(2, 1) = 0.5D0 + w(1) = 1.D0/6.D0 + xi(1, 2) = 0.D0 + xi(2, 2) = 0.5D0 + w(2) = w(1) + xi(1, 3) = 0.5D0 + xi(2, 3) = 0.D0 + w(3) = w(1) + else if ( n_points==4 ) then + if ( n_nodes==4 .or. n_nodes==8 .or. n_nodes==9 ) then + ! 2X2 GAUSS INTEGRATION POINTS FOR QUADRILATERAL + ! 43 + ! 12 + cn = 0.5773502691896260D0 + xi(1, 1) = -cn + xi(1, 2) = cn + xi(1, 3) = cn + xi(1, 4) = -cn + xi(2, 1) = -cn + xi(2, 2) = -cn + xi(2, 3) = cn + xi(2, 4) = cn + w(1) = 1.D0 + w(2) = 1.D0 + w(3) = 1.D0 + w(4) = 1.D0 + else if ( n_nodes==3 .or. n_nodes==6 ) then + ! xi integration points for triangle + xi(1, 1) = 1.D0/3.D0 + xi(2, 1) = xi(1, 1) + w(1) = -27.D0/96.D0 + xi(1, 2) = 0.6D0 + xi(2, 2) = 0.2D0 + w(2) = 25.D0/96.D0 + xi(1, 3) = 0.2D0 + xi(2, 3) = 0.6D0 + w(3) = w(2) + xi(1, 4) = 0.2D0 + xi(2, 4) = 0.2D0 + w(4) = w(2) + end if + + else if ( n_points==7 ) then + ! Quintic integration for triangle + xi(1,1) = 1.d0/3.d0 + xi(2,1) = xi(1,1) + w(1) = 0.1125d0 + xi(1,2) = 0.0597158717d0 + xi(2,2) = 0.4701420641d0 + w(2) = 0.0661970763d0 + xi(1,3) = xi(2,2) + xi(2,3) = xi(1,2) + w(3) = w(2) + xi(1,4) = xi(2,2) + xi(2,4) = xi(2,2) + w(4) = w(2) + xi(1,5) = 0.7974269853d0 + xi(2,5) = 0.1012865073d0 + w(5) = 0.0629695902d0 + xi(1,6) = xi(2,5) + xi(2,6) = xi(1,5) + w(6) = w(5) + xi(1,7) = xi(2,5) + xi(2,7) = xi(2,5) + w(7) = w(5) + else if ( n_points==9 ) then + ! 3X3 GAUSS INTEGRATION POINTS + ! 789 + ! 456 + ! 123 + cn = 0.7745966692414830D0 + xi(1, 1) = -cn + xi(1, 2) = 0.D0 + xi(1, 3) = cn + xi(1, 4) = -cn + xi(1, 5) = 0.D0 + xi(1, 6) = cn + xi(1, 7) = -cn + xi(1, 8) = 0.D0 + xi(1, 9) = cn + xi(2, 1) = -cn + xi(2, 2) = -cn + xi(2, 3) = -cn + xi(2, 4) = 0.D0 + xi(2, 5) = 0.D0 + xi(2, 6) = 0.D0 + xi(2, 7) = cn + xi(2, 8) = cn + xi(2, 9) = cn + w1 = 0.5555555555555560D0 + w2 = 0.8888888888888890D0 + w11 = w1*w1 + w12 = w1*w2 + w22 = w2*w2 + w(1) = w11 + w(2) = w12 + w(3) = w11 + w(4) = w12 + w(5) = w22 + w(6) = w12 + w(7) = w11 + w(8) = w12 + w(9) = w11 + end if + + return + + end subroutine abq_UEL_2D_integrationpoints + + subroutine abq_UEL_invert2d(A,A_inverse,determinant) + + double precision, intent(in) :: A(2,2) + double precision, intent(out) :: A_inverse(2,2) + double precision, intent(out) :: determinant + +! Compute inverse and determinant of 2x2 matrix + + determinant = A(1,1)*A(2,2) - A(1,2)*A(2,1) + + IF (determinant==0.d0) THEN + write(6,*) ' Error in subroutine abq_UEL_inver2d' + write(6,*) ' A 2x2 matrix has a zero determinant' + stop + endif + + A_inverse(1,1) = A(2,2) + A_inverse(1,2) = -A(1,2) + A_inverse(2,1) = -A(2,1) + A_inverse(2,2) = A(1,1) + + A_inverse = A_inverse / determinant + + + end subroutine abq_UEL_invert2d + + + subroutine abq_UEL_2D_shapefunctions(xi,n_nodes,f,df) + + implicit none + integer, intent(in) :: n_nodes + + double precision, intent(in) :: xi(2) + double precision, intent(out) :: f(*) + double precision, intent(out) :: df(9,2) + double precision g1, g2, g3, dg1, dg2, dg3 + double precision h1, h2, h3, dh1, dh2, dh3 + double precision z,dzdp, dzdq + + if ( n_nodes==3 ) then ! SHAPE FUNCTIONS FOR 3 NODED TRIANGLE + f(1) = xi(1) + f(2) = xi(2) + f(3) = 1.D0 - xi(1) - xi(2) + df(1, 1) = 1.D0 + df(1, 2) = 0.D0 + df(2, 1) = 0.D0 + df(2, 2) = 1.D0 + df(3, 1) = -1.D0 + df(3, 2) = -1.D0 + else if ( n_nodes==4 ) then + ! SHAPE FUNCTIONS FOR 4 NODED QUADRILATERAL + ! 43 + ! 12 + g1 = 0.5D0*(1.D0 - xi(1)) + g2 = 0.5D0*(1.D0 + xi(1)) + h1 = 0.5D0*(1.D0 - xi(2)) + h2 = 0.5D0*(1.D0 + xi(2)) + f(1) = g1*h1 + f(2) = g2*h1 + f(3) = g2*h2 + f(4) = g1*h2 + dg1 = -0.5D0 + dg2 = 0.5D0 + dh1 = -0.5D0 + dh2 = 0.5D0 + df(1, 1) = dg1*h1 + df(2, 1) = dg2*h1 + df(3, 1) = dg2*h2 + df(4, 1) = dg1*h2 + df(1, 2) = g1*dh1 + df(2, 2) = g2*dh1 + df(3, 2) = g2*dh2 + df(4, 2) = g1*dh2 + + else if ( n_nodes==6 ) then + + ! SHAPE FUNCTIONS FOR 6 NODED TRIANGLE + ! 3 + + ! 6 5 + + ! 1 4 2 + + ! P = L1 + ! Q = L2 + ! Z = 1 - P - Q = L3 + + z = 1.D0 - xi(1) - xi(2) + f(1) = (2.D0*xi(1) - 1.D0)*xi(1) + f(2) = (2.D0*xi(2) - 1.D0)*xi(2) + f(3) = (2.D0*z - 1.D0)*z + f(4) = 4.D0*xi(1)*xi(2) + f(5) = 4.D0*xi(2)*z + f(6) = 4.D0*xi(1)*z + dzdp = -1.D0 + dzdq = -1.D0 + df(1, 1) = 4.D0*xi(1) - 1.D0 + df(2, 1) = 0.D0 + df(3, 1) = 4.D0*z*dzdp - dzdp + df(4, 1) = 4.D0*xi(2) + df(5, 1) = 4.D0*xi(2)*dzdp + df(6, 1) = 4.D0*z + 4.D0*xi(1)*dzdp + df(1, 2) = 0.D0 + df(2, 2) = 4.D0*xi(2) - 1.D0 + df(3, 2) = 4.D0*z*dzdq - dzdq + df(4, 2) = 4.D0*xi(1) + df(5, 2) = 4.D0*z + 4.D0*xi(2)*dzdq + df(6, 2) = 4.D0*xi(1)*dzdq + + else if ( n_nodes==8 ) then + ! SHAPE FUNCTIONS FOR 8 NODED SERENDIPITY ELEMENT + f(1) = -0.25*(1.-xi(1))*(1.-xi(2))*(1.+xi(1)+xi(2)); + f(2) = 0.25*(1.+xi(1))*(1.-xi(2))*(xi(1)-xi(2)-1.); + f(3) = 0.25*(1.+xi(1))*(1.+xi(2))*(xi(1)+xi(2)-1.); + f(4) = 0.25*(1.-xi(1))*(1.+xi(2))*(xi(2)-xi(1)-1.); + f(5) = 0.5*(1.-xi(1)*xi(1))*(1.-xi(2)); + f(6) = 0.5*(1.+xi(1))*(1.-xi(2)*xi(2)); + f(7) = 0.5*(1.-xi(1)*xi(1))*(1.+xi(2)); + f(8) = 0.5*(1.-xi(1))*(1.-xi(2)*xi(2)); + df(1,1) = 0.25*(1.-xi(2))*(2.*xi(1)+xi(2)); + df(1,2) = 0.25*(1.-xi(1))*(xi(1)+2.*xi(2)); + df(2,1) = 0.25*(1.-xi(2))*(2.*xi(1)-xi(2)); + df(2,2) = 0.25*(1.+xi(1))*(2.*xi(2)-xi(1)); + df(3,1) = 0.25*(1.+xi(2))*(2.*xi(1)+xi(2)); + df(3,2) = 0.25*(1.+xi(1))*(2.*xi(2)+xi(1)); + df(4,1) = 0.25*(1.+xi(2))*(2.*xi(1)-xi(2)); + df(4,2) = 0.25*(1.-xi(1))*(2.*xi(2)-xi(1)); + df(5,1) = -xi(1)*(1.-xi(2)); + df(5,2) = -0.5*(1.-xi(1)*xi(1)); + df(6,1) = 0.5*(1.-xi(2)*xi(2)); + df(6,2) = -(1.+xi(1))*xi(2); + df(7,1) = -xi(1)*(1.+xi(2)); + df(7,2) = 0.5*(1.-xi(1)*xi(1)); + df(8,1) = -0.5*(1.-xi(2)*xi(2)); + df(8,2) = -(1.-xi(1))*xi(2); + else if ( n_nodes==9 ) then + ! SHAPE FUNCTIONS FOR 9 NODED LAGRANGIAN ELEMENT + ! 789 + ! 456 + ! 123 + g1 = -.5D0*xi(1)*(1.D0 - xi(1)) + g2 = (1.D0 - xi(1))*(1.D0 + xi(1)) + g3 = .5D0*xi(1)*(1.D0 + xi(1)) + h1 = -.5D0*xi(2)*(1.D0 - xi(2)) + h2 = (1.D0 - xi(2))*(1.D0 + xi(2)) + h3 = .5D0*xi(2)*(1.D0 + xi(2)) + dg1 = xi(1) - 0.5d0 + dg2 = -2.d0*xi(1) + dg3 = xi(1) + 0.5d0 + dh1 = xi(2)-0.5d0 + dh2 = -2.d0*xi(2) + dh3 = xi(2) + 0.5d0 + f(1) = g1*h1 + f(2) = g2*h1 + f(3) = g3*h1 + f(4) = g1*h2 + f(5) = g2*h2 + f(6) = g3*h2 + f(7) = g1*h3 + f(8) = g2*h3 + f(9) = g3*h3 + df(1,1) = dg1*h1 + df(1,2) = g1*dh1 + df(2,1) = dg2*h1 + df(2,2) = g2*dh1 + df(3,1) = dg3*h1 + df(3,2) = g3*dh1 + df(4,1) = dg1*h2 + df(4,2) = g1*dh2 + df(5,1) = dg2*h2 + df(5,2) = g2*dh2 + df(6,1) = dg3*h2 + df(6,2) = g3*dh2 + df(7,1) = dg1*h3 + df(7,2) = g1*dh3 + df(8,1) = dg2*h3 + df(8,2) = g2*dh3 + df(9,1) = dg3*h3 + df(9,2) = g3*dh3 + end if + + end subroutine abq_UEL_2D_shapefunctions + + + subroutine abq_UEL_1D_integrationpoints(n_points, n_nodes, xi, w) + + + implicit none + integer, intent(in) :: n_points + integer, intent(in) :: n_nodes + + double precision, intent(out) :: xi(*) + double precision, intent(out) :: w(*) + + integer :: i,j,k,n + + double precision x1D(4), w1D(4) + + + + select case ( n_points ) + case (2) + xi(1) = .5773502691896257D+00 + xi(2) = -.5773502691896257D+00 + w(1) = .1000000000000000D+01 + w(2) = .1000000000000000D+01 + return + case (3) + xi(1) = 0.7745966692414834D+00 + xi(2) = .0000000000000000D+00 + xi(3) = -.7745966692414834D+00 + w(1) = .5555555555555556D+00 + w(2) = .8888888888888888D+00 + w(3) = .5555555555555556D+00 + return + case (4) + xi(1) = .8611363115940526D+00 + xi(2) = .3399810435848563D+00 + xi(3) = -.3399810435848563D+00 + xi(4) = -.8611363115940526D+00 + w(1) = .3478548451374538D+00 + w(2) = .6521451548625461D+00 + w(3) = .6521451548625461D+00 + w(4) = .3478548451374538D+00 + return + case (5) + xi(1) = .9061798459386640D+00 + xi(2) = .5384693101056831D+00 + xi(3) = .0000000000000000D+00 + xi(4) = -.5384693101056831D+00 + xi(5) = -.9061798459386640D+00 + w(1) = .2369268850561891D+00 + w(2) = .4786286704993665D+00 + w(3) = .5688888888888889D+00 + w(4) = .4786286704993665D+00 + w(5) = .2369268850561891D+00 + return + case (6) + xi(1) = .9324695142031521D+00 + xi(2) = .6612093864662645D+00 + xi(3) = .2386191860831969D+00 + xi(4) = -.2386191860831969D+00 + xi(5) = -.6612093864662645D+00 + xi(6) = -.9324695142031521D+00 + w(1) = .1713244923791703D+00 + w(2) = .3607615730481386D+00 + w(3) = .4679139345726910D+00 + w(4) = .4679139345726910D+00 + w(5) = .3607615730481386D+00 + w(6) = .1713244923791703D+00 + return + case DEFAULT + write(6,*)'Error in subroutine abq_UEL_1D_integrationpoints' + write(6,*) ' Invalid number of integration points for 1D' + write(6,*) ' n_points must be between 1 and 6' + stop + end select + + + + + + + + end subroutine ABQ_UEL_1D_integrationpoints + + + + subroutine abq_facenodes_2D(nelnodes,face,list,nfacenodes) + + implicit none + + integer, intent (in) :: nelnodes + integer, intent (in) :: face + integer, intent (out) :: list(*) + integer, intent (out) :: nfacenodes + ! + ! Subroutine to return list of nodes on an element face for standard 2D solid elements + ! + integer :: i3(3) + integer :: i4(4) + + i3(1:3) = [2,3,1] + i4(1:4) = [2,3,4,1] + + if (nelnodes == 3) then + nfacenodes = 2 + list(1) = face + list(2) = i3(face) + else if (nelnodes == 4) then + nfacenodes = 2 + list(1) = face + list(2) = i4(face) + else if (nelnodes == 6) then + nfacenodes = 3 + list(1) = face + list(2) = i3(face) + list(3) = face+3 + else if (nelnodes == 8) then + nfacenodes = 3 + list(1) = face + list(2) = i4(face) + list(3) = face+4 + else if (nelnodes == 9) then + nfacenodes = 3 + if (face==1) list(1:3) = (/1,3,2/) + if (face==2) list(1:3) = (/3,9,6/) + if (face==3) list(1:3) = (/9,7,8/) + if (face==4) list(1:3) = (/7,1,4/) + endif + + end subroutine abq_facenodes_2d + + subroutine abq_inverse_LU(Ain,A_inverse,n) ! Compute the inverse of an arbitrary matrix by LU decomposition + + implicit none + + integer, intent(in) :: n + + double precision, intent(in) :: Ain(n,n) + double precision, intent(out) :: A_inverse(n,n) + + double precision :: A(n,n), L(n,n), U(n,n), b(n), d(n), x(n) + double precision :: coeff + integer :: i, j, k + + A(1:n,1:n) = Ain(1:n,1:n) + L=0.d0 + U=0.d0 + b=0.d0 + + do k=1, n-1 + do i=k+1,n + coeff=a(i,k)/a(k,k) + L(i,k) = coeff + A(i,k+1:n) = A(i,k+1:n)-coeff*A(k,k+1:n) + end do + end do + + forall (i=1:n) L(i,i) = 1.d0 + forall (j=1:n) U(1:j,j) = A(1:j,j) + do k=1,n + b(k)=1.d0 + d(1) = b(1) + do i=2,n + d(i)=b(i) + d(i) = d(i) - dot_product(L(i,1:i-1),d(1:i-1)) + end do + x(n)=d(n)/U(n,n) + do i = n-1,1,-1 + x(i) = d(i) + x(i)=x(i)-dot_product(U(i,i+1:n),x(i+1:n)) + x(i) = x(i)/U(i,i) + end do + A_inverse(1:n,k) = x(1:n) + b(k)=0.d0 + end do + + end subroutine abq_inverse_LU + + diff --git a/main.f90 b/main.f90 index 5e0e013..1b383b6 100644 --- a/main.f90 +++ b/main.f90 @@ -115,14 +115,14 @@ program en234fea ! Homework 9 - McCormick model with 1 element - infil = 'input_files/Abaqus_vumat_McCormick.in' - outfil = 'Output_files/Abaqus_vumat_McCormick.out' +! infil = 'input_files/Abaqus_vumat_McCormick.in' +! outfil = 'Output_files/Abaqus_vumat_McCormick.out' ! Homework 10 - Continuum beam element solution to end loaded cantilever beam -! infil = 'input_files/Abaqus_uel_continuum_beam.in' -! outfil = 'Output_files/Abaqus_uel_continuum_beam.out' + infil = 'input_files/Abaqus_uel_continuum_beam.in' + outfil = 'Output_files/Abaqus_uel_continuum_beam.out' infil = trim(root_directory)//trim(infil) outfil = trim(root_directory)//trim(outfil) From 2899f65b9bb310d0ecd1cf66f96d7bcf2120e44f Mon Sep 17 00:00:00 2001 From: dongli314 <32434921+dongli314@users.noreply.github.com> Date: Fri, 1 Dec 2017 19:26:39 -0500 Subject: [PATCH 10/10] HW10-Continuum beam element UEL subroutine for continuum beam element --- abaqus_uel_2D.for | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/abaqus_uel_2D.for b/abaqus_uel_2D.for index bc9e71a..15b08ed 100644 --- a/abaqus_uel_2D.for +++ b/abaqus_uel_2D.for @@ -96,7 +96,7 @@ integer :: i,j,n_points,kint, nfacenodes, ipoin, ksize integer :: face_node_list(3) ! List of nodes on an element face ! - double precision :: xi(2,9) ! Area integration points + double precision :: xi(2,9), xi1(9) ! Area integration points double precision :: w(9) ! Area integration weights double precision :: N(9) ! 2D shape functions double precision :: dNdxi(9,2) ! 2D shape function derivatives @@ -134,7 +134,9 @@ if (NNODE == 9) n_points = 9 ! Quadratic rect ! Write your code for a 2D element below - call abq_UEL_1D_integrationpoints(n_points, NNODE, xi, w) + call abq_UEL_1D_integrationpoints(n_points, NNODE, xi1, w) + xi = 0.d0 + xi(2,1:n_points) = xi1(1:n_points) if (MLVARX<2*NNODE) then write(6,*) ' Error in abaqus UEL '