This code is a continuum damage mechanics (CDM) material model for use with the Abaqus finite element code. This is a research code which aims to provide an accurate representation of mesoscale damage modes in fiber-reinforced polymer materials.
The CDM material model is implemented as an Abaqus VUMAT for the simulation of matrix cracks formed under tensile, compressive, and shear loading conditions and fiber fracture under tensile and compressive loading. A smeared crack approach is used to accurately represent the kinematics of matrix cracks rotating with the deformed material.
This software may be used, reproduced, and provided to others only as permitted under the terms of the agreement under which it was acquired from the U.S. Government. Neither title to, nor ownership of, the software is hereby transferred. This notice shall remain on all copies of the software.
Copyright 2016 United States Government as represented by the Administrator of the National Aeronautics and Space Administration. No copyright is claimed in the United States under Title 17, U.S. Code. All Other Rights Reserved.
Publications that use this code include:
- Cheryl A. Rose, et al. "Analysis Methods for Progressive Damage of Composite Structures" NASA/TM-2013-218024, July 2013.
- Frank A. Leone Jr. "Deformation gradient tensor decomposition for representing matrix cracks in fiber-reinforced materials" Composites Part A (2015) 76:334-341.
- Frank A. Leone Jr. "Representing matrix cracks through decomposition of the deformation gradient tensor in continuum damage mechanics methods" ICCM20.
For any questions, please contact the developers:
- Frank Leone | frank.a.leone@nasa.gov | (W) 757-864-3050
- Andrew Bergan | andrew.c.bergan@nasa.gov | (W) 757-864-3744
- Carlos Dávila | carlos.g.davila@nasa.gov | (W) 757-864-9130
- Getting started
- Model features
- Elements
- Material properties
- State variables
- Advanced debugging
- Summary of tests classes
- Known issues
The user subroutine source code is located in the for directory. The main entry point is CompDam_DGD.for.
Intel Fortran Compiler version 11.1 or newer is required to compile the code. It is recommended that Abaqus 6.14-1 or newer is used with this code.
The abaqus_v6.env file must have /fpp, /Qmkl:sequential, and /free in the ifort command where the format for Windows is used. The corresponding Linux format is: -fpp, -free, and -mkl=sequential. The /fpp option enables the Fortran preprocessor, which is required for the code to compile correctly. The /free option sets the compiler to free-formatting for the source code files. The /Qmkl:sequential enables the Intel Math Kernel Library (MKL), which provides optimized and verified functions for many mathematical operations. The MKL is used in this code for calculating eigenvalues and eigenvectors.
A sample environment file is provided in the tests directory for Windows and Linux systems.
This code is an Abaqus/Explicit VUMAT. Please refer to the Abaqus documentation for the general instructions on how to submit finite element analyses using user subroutines. Please see the example input file statements for details on how to interface with this particular VUMAT subroutine.
Analyses with this code must be run in double precision. Some of the code has double precision statements and variables hard-coded, so if Abaqus/Explicit is run in single precision, compile-time errors will arise. When submitting an Abaqus/Explicit job from the command line, double precision can be specified by including the command line argument double=both.
For example, run the test model test_C3D8R_elastic_fiberTension in the tests directory with the following command:
abaqus job=test_C3D8R_elastic_fiberTension user=../for/CompDam_DGD.for double=both
Example 1, using an external material properties file:
*Section controls, name=, distortion control=YES, element deletion=YES
**
*Material, name=IM7-8552
*Density
1.57e-09,
*Depvar, delete=11
** the above delete statement is optional
19,
1, CDM_d2
2, CDM_Fb1
3, CDM_Fb2
4, CDM_Fb3
5, CDM_B
6, CDM_Lc
7, CDM_rfT
8, CDM_d1
9, CDM_FImT
10, CDM_alpha
11, CDM_STATUS
12, CDM_Plas12
13, CDM_Inel12
14, CDM_d_comp_init
15, CDM_slide1
16, CDM_slide2
17, CDM_rfC
18, CDM_d1T
19, CDM_d1C
*User material, constants=3
** 1 2 3
** feature flags, , thickness
100001, , 0.1
**
*Initial Conditions, type=SOLUTION
elset_name, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0,
0.d0, 0.d0, -999, 1, 0.d0, 0.d0, 0.d0,
0.d0, 0.d0, 0.d0, 0.d0, 0.d0
** - alpha can be pre-set by setting SV10 to the desired angle.
** Setting alpha = -999 will evaluate cracks every 10 degrees in
** the 2-3 plane to find the correct crack initiation angle.
** In each step, NLGEOM=YES must be used. This is the default setting.
Example 2, using an input deck command:
*Section controls, name=, distortion control=YES, element deletion=YES
**
*Material, name=IM7-8552
*Density
1.57e-09,
*Depvar, delete=11
** the above delete statement is optional
19,
1, CDM_d2
2, CDM_Fb1
3, CDM_Fb2
4, CDM_Fb3
5, CDM_B
6, CDM_Lc
7, CDM_rfT
8, CDM_d1
9, CDM_FImT
10, CDM_alpha
11, CDM_STATUS
12, CDM_Plas12
13, CDM_Inel12
14, CDM_d_comp_init
15, CDM_slide1
16, CDM_slide2
17, CDM_rfC
18, CDM_d1T
19, CDM_d1C
*User material, constants=40
** 1 2 3 4 5 6 7 8
** feature flags, , thickness, 4, 5, 6, 7, 8
100001, , 0.1, , , , , ,
**
** 9 10 11 12 13 14 15 16
** E1, E2, G12, nu12, nu23, YT, SL GYT,
171420.0, 9080.0, 5290.0, 0.32, 0.52, 62.3, 92.30, 0.277,
**
** 17 18 19 20 21 22 23 24
** GSL, eta_BK, YC, alpha0 E3, G13, G23, nu13,
0.788, 1.634, 199.8, 0.925, , , , ,
**
** 25 26 27 28 29 30 31 32
** alpha11, alpha22, alpha_PL, n_PL, XT, fXT, GXT, fGXT,
-5.5d-6, 2.58d-5, , , 2326.2, 0.2, 133.3, 0.5,
**
** 33 34 35 36 37 38 39 40
** XC, fXC, GXC, fGXC, rsvd, rsvd, rsvd, mu
1200.1, , , , , , , 0.3
**
*Initial Conditions, type=SOLUTION
elset_name, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0, 0.d0,
0.d0, 0.d0, -999, 1, 0.d0, 0.d0, 0.d0,
0.d0, 0.d0, 0.d0, 0.d0, 0.d0
** - alpha can be pre-set by setting SV10 to the desired angle.
** Setting alpha = -999 will evaluate cracks every 10 degrees in
** the 2-3 plane to find the correct crack initiation angle.
** In each step, NLGEOM=YES must be used. This is the default setting.
Test cases are available in the tests directory. The tests are useful for demonstrating the capabilities of the VUMAT as well as to verify that the code performs as intended. Try running some of the test cases to see how the code works. The test cases can be submitted as a typical Abaqus job using the Abaqus command line arguments.
CompDam_DGD can be built into a shared library file. Follow these steps:
- Place a copy of the Abaqus environment file (with the compiler flags specified) in the
fordirectory - In Linux, and when using Abaqus versions prior to 2017, rename
CompDam_DGD.fortoCompDam_DGD.f - From the
fordirectory, run:
abaqus make library=CompDam_DGD
This command will create shared libraries for the operating system it is executed on (.dll for Windows and .so for Linux).
When using a pre-compiled shared library, it is only necessary to specify the location of the shared library files in the environment file (the compiler options are not required). To run an analysis using a shared library, add usub_lib_dir = <full path to shared library file> to the Abaqus environment file in the Abaqus working directory.
The CompDam_DGD material model implements a variety of features that can be enabled or disabled by the user. An overview of these features is provided in this section. The material properties required for each feature are listed. References are provided to more detailed discussions of the theoretical framework for each feature.
The composite materials modeled with CompDam_DGD can be defined assuming either transverse isotropy or orthotropy. For a transversely isotropic material definition, the following properties must be defined: E1, E2, G12, v12, and v23. For an orthotropic material definition, the following additional properties must be defined: E2, G13, G23, and v13.
Tensile and compressive matrix damage is modeled by embedding cohesive laws to represent cracks in the material according to the deformation gradient decomposition method of Leone (2015). The mixed-mode behavior of matrix cracks is defined according to the Benzeggagh-Kenane law. The initiation of compressive matrix cracks accounts for friction on the potential crack surface according to the LaRC04 failure criterion.
The following material properties are required for the prediction of matrix damage: YT, SL, GYT, GSL, eta_BK, YC, and alpha0.
The thermal strains are calculated by multiplying the 1-, 2-, and 3-direction coefficients of thermal expansion by the current ΔT, as provided by the Abaqus solver. The thermal strains are subtracted from the current total strain.
The required material properties are the coefficients of thermal expansion in the 1 and 2 directions. It is assumed that the 2- and 3-direction coefficients of thermal expansion are equal.
Hygroscopic strains are not accounted for. If the effects of hygroscopic expansion are to be modeled, it is recommended to smear the hygroscopic and thermal expansion coefficients to approximate the response using the solver-provided ΔT.
Shear nonlinearity is modeled in the 1-2 plane according to the Ramberg-Osgood equation, with its parameters selected to fit experimental data. As applied herein, the Ramberg-Ogsood equation is written in the following form:
γ12 = [τ12 + αPLsign(τ12)|τ12|nPL]/G12
where γ12 is the shear strain and τ12 is the shear stress. Prior to the initiation of matrix damage (i.e., d2), the nonlinear response due to the above is equation is plastic, and the unloading/reloading slope is unchanged. The required material inputs are the two parameters in the above equation: αPL and nPL.
A continuum damage mechanics model similar to the work of Maimí et al. (2007) is used. The model utilizes a non-interacting maximum strain failure criterion, and bilinear softening after the initiation of failure. The area under the stress-strain curve is equal to the fracture toughness divided by the element length. The required material properties are: XT, fXT, GXT, and fGXT, where fXT and fGXT are ratios of strength and fracture toughness for bilinear softening, analogous to the n and m terms in Dávila et al. (2009).
Same model as in tension, but for compression. Assumes maximum strain failure criterion and bilinear softening. The required material properties are: XC, fXC, GXC, fGXC.
Load reversal assumptions from Maimí et al. (2007).
Friction is modeled on the damaged fraction of the cross-sectional area of DGD cracks using the approach of Alfano and Sacco (2006). The coefficient of friction μ must be defined to account for friction on the failed crack surface.
The amount of sliding which has taken place in the longitudinal and transverse directions are stored in state variables CDM_slide1 and CDM_slide2.
The strain is calculated using the deformation gradient tensor provided by the Abaqus solver. The default strain definition used is the Green-Lagrange strain:
E = (FTF - I)/2
Hooke's law is applied using the Green-Lagrange strain to calculate the 2nd Piola-Kirchhoff stress S.
CompDam_DGD has been developed and tested using the Abaqus three-dimensional, reduced-integration C3D8R solid and S4R shell elements. Limited testing has been performed using the CPS4R plane stress element, the SC8R continuum shell element, and the fully-integrated C3D8 solid element.
Because CompDam_DGD is a material model, it is expected to be compatible with structural elements generally. However, users are advised to perform tests with any previously untested element types before proceeding to use CompDam_DGD in larger structural models.
A set of material properties must be defined for the material of interest. This section describes how to specify the material properties.
Material properties can be defined in the input deck or in a separate .props file. Definition of the material properties in a .props file is more convenient and generally preferred since it isolates the material properties from the structural model definition.
Using a .props file is a versatile means of defining material properties. The subroutine looks for a file named as jobName_materialName where the job name is the name of the Abaqus job (default is input deck name) and the material is name assigned as *Material, Name=materialName in the input deck. If no file is found, then the subroutine looks for materialName.props. The .props file must be located in the Abaqus working directory.
The .props should contain one property per line with the format [NAME] = [VALUE] where the name is symbolic name for the property and the value is assigned value for the property. Blank lines or commented text (denoted by //) is ignored. See the table of material properties for a complete list of material property symbolic names, acceptable values, and recommended test methods for characterizing the properties.
When a .props is used to define the material properties, the feature flags and thickness still must be defined in the input deck.
Material properties can be defined in the input deck. Any optional material property can be left blank and the corresponding feature(s) will be disabled. The ordering of the material properties for the input deck definition is given in the first (#) column of the table of material properties.
| # | Symbol | Name | Description | Units | Admissible values | Recommeded Test |
|---|---|---|---|---|---|---|
| 9 | E1 | E1 | Longitudinal Young's modulus | F/L2 | 0 < E1 < ∞ | ASTM D3039 |
| 10 | E2 | E2 | Transverse Young's modulus | F/L2 | 0 < E2 < ∞ | ASTM D3039 |
| 11 | G12 | G12 | In-plane Shear modulus | F/L2 | 0 < G12 < ∞ | ASTM D3039 |
| 12 | ν12 | nu12 | Major Poisson's ratio | - | 0 ≤ ν12 ≤ 1 | ASTM D3039 |
| 13 | ν23 | nu23 | Minor Poisson's ratio | - | 0 ≤ ν23 ≤ 1 | |
| === | ||||||
| 14 | YT | YT | Transverse tensile strength, in-situ | F/L2 | 0 < YT < ∞ | ASTM D3039 |
| 15 | SL | SL | Shear strength, in-situ | F/L2 | 0 < SL < ∞ | |
| 16 | GIc | GYT | Mode I fracture toughness | F/L | 0 < GIc < ∞ | ASTM D5528 |
| 17 | GIIc | GSL | Mode II fracture toughness | F/L | 0 < GIIc < ∞ | ASTM D7905 |
| 18 | η | eta_BK | BK exponent for mode-mixity | - | 1 ≤ η < ∞ | |
| 19 | YC | YC | Transverse compressive strength | F/L2 | 0 < YC < ∞ | ASTM D3410 |
| 20 | α0 | alpha0 | Fracture plane angle for pure trans. comp. | Radians | 0 ≤ α0 ≤ π/2 | |
| ------ | ||||||
| 21 | E3 | E3 | 3-Direction Young's modulus | F/L2 | 0 < E3 < ∞ | |
| 22 | G13 | G13 | Shear modulus in 1-3 plane | F/L2 | 0 < G13 < ∞ | |
| 23 | G23 | G23 | Shear modulus in 1-2 plane | F/L2 | 0 < G23 < ∞ | |
| 24 | ν13 | nu13 | Poisson's ratio in 2-3 plane | - | 0 ≤ ν13 ≤ 1 | |
| ------ | ||||||
| 25 | α11 | alpha11 | Coefficient of long. thermal expansion | /° | -1 ≤ α11 ≤ 1 | |
| 26 | α22 | alpha11 | Coefficient of tran. thermal expansion | /° | -1 ≤ α22 ≤ 1 | |
| ------ | ||||||
| 27 | αPL | alpha_PL | Nonlinear shear parameter | (F/L2)1-nPL | 0 ≤ αPL < ∞ | |
| 28 | nPL | n_PL | Nonlinear shear parameter | - | 0 ≤ nPL < ∞ | |
| ------ | ||||||
| 29 | XT | XT | Long. tensile strength | F/L2 | 0 < XT < ∞ | ASTM D3039 |
| 30 | fXT | fXT | Long. tensile strength ratio | - | 0 ≤ fXT ≤ 1 | |
| 31 | GXT | GXT | Long. tensile fracture toughness | F/L | 0 < GXT < ∞ | |
| 32 | fGXT | fGXT | Long. tensile fracture toughness ratio | - | 0 ≤ fGXT ≤ 1 | |
| ------ | ||||||
| 33 | XC | XC | Longitudinal compressive strength | F/L2 | 0 < XC < ∞ | ASTM D3410 |
| 34 | fXC | fXC | Long. compression strength ratio | - | 0 ≤ fXC ≤ 1 | |
| 35 | GXC | GXC | Long. compression fracture toughness | F/L | 0 < GXC < ∞ | |
| 36 | fGXC | fGXC | Long. compression fracture toughness ratio | - | 0 ≤ fGXC ≤ 1 | |
| ------ | ||||||
| 40 | μ | mu | Coefficient of friction | - | 0 ≤ μ ≤ 1 |
Notes:
- The first five inputs (above the ===) are required
- Properties for each feature are grouped and separated by ------
- The number in the first column corresponds to the property order when defined in the input deck
- Properties not listed in the table above (see next section):
- feature flags
- reserved
- thickness
- reserved
- reserved
- reserved
- reserved
- reserved
- ∞ is calculated with the Fortran intrinsic
Hugefor double precision
The feature flags and thickness are defined in the input deck on the material property data lines. These properties must be defined in the input deck whether the other material properties are defined via the .props file or via the input deck. While feature flags and thickness are not material properties per se, they are used in controlling the behavior of the material model.
Model features can be enabled or disabled by two methods. The first method is specifying only the material properties required for the features you would like to enable. CompDam_DGD disables any feature for which all of the required material properties have not been assigned. If an incomplete set of material properties are defined for a feature, a warning is issued.
The second method is by specifying the status of each feature directly as a material property in the input deck. Each feature of the subroutine is controlled by a position in an integer, where 0 is disabled and 1 is enabled. The positions correspond to the features as follows:
- Position 1: Matrix damage
- Position 2: Shear nonlinearity
- Position 3: Fiber tensile damage
- Position 4: Fiber compression damage
- Position 5: reserved
- Position 6: Friction
For example, 101000 indicates that the model will run with matrix damage and fiber tension damage enabled; 110001 indicates that the model will run with matrix damage, in-plane shear nonlinearity, and friction.
Length along the thickness-direction associated with the current integration point
The table below lists all of the state variables in the model. The model requires a minimum of 18 state variables. Additional state variables are defined depending on which (if any) fiber compression features are enabled. When fiber compression (feature flag position 4) is used, 19 state variables are required.
| # | Name | Description |
|---|---|---|
| 1 | CDM_d2 |
d2, Matrix cohesive damage variable |
| 2 | CDM_Fb1 |
Cohesive Opening Displacement 1 - fiber-direction shear |
| 3 | CDM_Fb2 |
Cohesive Opening Displacement 2 - normal |
| 4 | CDM_Fb3 |
Cohesive Opening Displacement 3 - transverse shear |
| 5 | CDM_B |
Mode Mixity (GII / (GI + GII)) |
| 6 | CDM_Lc |
Characteristic Element Length |
| 7 | CDM_rfT |
Fiber tensile threshold (maximum strain) |
| 8 | CDM_d1 |
d1, Fiber damage variable |
| 9 | CDM_FImT |
Matrix failure criterion (del/del_0) |
| 10 | CDM_alpha |
alpha, the cohesive surface normal [degrees, integer] |
| 11 | CDM_STATUS |
STATUS (for element deletion) |
| 12 | CDM_Plas12 |
1-2 plastic strain |
| 13 | CDM_Inel12 |
1-2 inelastic strain |
| 14 | CDM_d_comp_init |
Normal cohesive displacement at initiation |
| 15 | CDM_slide1 |
Cohesive sliding displacement along 1-direction |
| 16 | CDM_slide2 |
Cohesive sliding displacement along 2-direction |
| 17 | CDM_rfC |
Fiber compression failure threshold |
| 18 | CDM_d1T |
Fiber tension damage variable |
| --- | ------------------ | ----------------------------------------------------------------------- |
| 19 | CDM_d1C |
Fiber compression damage variable |
All state variables should be initialized using the *Initial conditions command. All state variables should be initialized as zero, except CDM_alpha and CDM_STATUS.
The initial condition for CDM_alpha can be used to specify a predefined angle for the cohesive surface normal. To specify a predefined CDM_alpha, set the initial condition for CDM_alpha to an integer (degrees). The range of valid values for CDM_alpha depends on the aspect ratio of the element, but values in the range of 0 to 90 degrees are always valid. Setting CDM_alpha to -999 will make the subroutine evaluate cracks every 10 degrees in the 2-3 plane to find the correct crack initiation angle. Note that CDM_alpha is measured from the 2-axis rotating about the 1-direction.
Since CDM_STATUS is used for element deletion, always initialize CDM_STATUS to 1.
Using an interactive debugger helps to identify issues in the Fortran code. Abaqus knowledge base article QA00000007986 describes the details involved. The following is a quick-start guide for direct application to CompDam.
Several statements for debugging need to be uncommented in the environment file. Follow these steps:
- Copy your system environment file (
SIMULIA\Abaqus\6.14-1\SMA\site) to your local working directory. For the example below, copy the environment file to thetestsdirectory. - Edit the local environment file: uncomment lines that end with
# <-- Debugging,# <-- Debug symbols, and# <-- Optimization Debugging
Run the job with the -debug and -explicit arguments. For example:
abaqus -j test_C3D8R_fiberTension -user ../for/CompDam_DGD.for -double both -debug -explicit
This command should open the Visual Studio debugging software automatically. Open the source file(s) to debug. At a minimum, open the file with the subroutine entry point for/CompDam_DGD.for. Set a break point by clicking in the shaded column on the left edge of the viewer. The break point will halt execution. Press F5 to start the solver. When the break point is reached, a yellow arrow will appear and code execution will pause. Press F5 to continue to the next break point, press F11 to execute the next line of code following execution into function calls (Step Into), or press F10 to execute the next line of code but not follow execution into function calls (Step Over).
To stop execution, close the visual studio window. Choose stop debugging and do not save your changes.
More tips on debugging Fortran programs from Intel.
This section includes a brief summary of each test implemented in the tests folder. The input deck file names briefly describe the test. All of the input decks start with test_<elementType>_ and end with a few words describing the test. A more detailed description for each is given below:
- elastic_fiberTension: Demonstrates the elastic response in the 1-direction under prescribed extension. The 1-direction stress-strain curve has the modulus E1.
- elastic_matrixTension: Demonstrates the elastic response in the 2-direction under prescribed extension. The 2-direction stress-strain curve has the modulus E2.
- fiberCompression_CDM: Demonstrates the constitutive response in the 1-direction under prescribed shortening. The 1-direction stress-strain curve is trilinear.
- fiberLoadReversal: Demonstrates the constitutive response in the 1-direction under prescribed extension and shortening reversals. The 1-direction stress-strain curve shows the intended behavior under load reversal.
- fiberTension: Demonstrates the constitutive response in the 1-direction under prescribed extension. The 1-direction stress-strain curve is trilinear.
- matrixTension: Demonstrates the constitutive response in the 2-direction under prescribed extension. The 2-direction stress-strain curve is bilinear.
- nonlinearShear12: Demonstrates the constitutive response under prescribed simple shear. Shows the response under load reversal.
- simpleShear12: Demonstrates the constitutive response under prescribed simple shear. The shear stress-strain curve is bilinear.
- simpleShear12friction: Demonstrates the constitutive response under prescribed simple shear with friction enabled. The element is loaded under transverse compression and then sheared. Shows the friction-induced stresses.
We invite your contributions to CompDam_DGD! Please submit contributions (including a test case) with pull requests so that we can reproduce the behavior of interest. Commits history should be clean. Please contact the developers if you would like to make a major contribution to this repository.