Damage Evolution and Element Removal for Fiber-Reinforced Composites

The damage evolution capability for fiber-reinforced materials in Abaqus:

assumes that damage is characterized by progressive degradation of material stiffness, leading to material failure;
requires linearly elastic behavior of the undamaged material (see Linear Elastic Behavior);
can be used with the Hashin damage initiation criterion for unidirectional fiber-reinforced composites to take into account four different failure modes: fiber tension, fiber compression, matrix tension, and matrix compression, with four damage variables to describe damage for each failure mode;
can be used with the LaRC05 damage initiation criterion for unidirectional fiber-reinforced composites in Abaqus/Standard when used with enriched elements to model discontinuities (such as cracks) in an extended finite element method (XFEM) analysis (see Modeling Discontinuities as an Enriched Feature Using the Extended Finite Element Method);
can be used with the ply fabric damage initiation criterion for bidirectional fabric-reinforced composites to take into account five different failure modes: four modes of fiber failure and one of matrix shear failure, with five damage variables to describe damage for each failure mode;
is based on energy dissipation during the damage process;
offers options for what occurs upon failure, including the removal of elements from the mesh; and
can be used with a viscous regularization of the constitutive equations to improve the convergence rate in the softening regime if used with the Hashin damage initiation criterion.

This page discusses:

Damage Evolution for the Hashin Model
Damage Evolution for the Ply Fabric Model
Maximum Degradation and Choice of Element Removal
Viscous Regularization
Material Damping
Elements
Output
References

Damage Evolution for the Hashin Model

The previous section (Damage Initiation for Fiber-Reinforced Composites) discussed the damage initiation in plane stress fiber-reinforced composites. This section will discuss the post-damage initiation behavior for cases in which a damage evolution model has been specified. Prior to damage initiation the material is linearly elastic, with the stiffness matrix of a plane stress orthotropic material. Thereafter, the response of the material is computed from

σ = C_{d} ε^{e l},

where $ε^{e l}$ is the elastic strain and $C_{d}$ is the damaged elasticity matrix, which has the form

C_{d} = \frac{1}{D} [\begin{matrix} (1 - d_{f}) E_{1} & (1 - d_{f}) (1 - d_{m}) ν_{21} E_{1} & 0 \\ (1 - d_{f}) (1 - d_{m}) ν_{12} E_{2} & (1 - d_{m}) E_{2} & 0 \\ 0 & 0 & (1 - d_{s}) G D \end{matrix}],

where $D = 1 - (1 - d_{f}) (1 - d_{m}) ν_{12} ν_{21}$ , $d_{f}$ reflects the current state of fiber damage, $d_{m}$ reflects the current state of matrix damage, $d_{s}$ reflects the current state of shear damage, $E_{1}$ is the Young's modulus in the fiber direction, $E_{2}$ is the Young's modulus in the matrix direction, $G$ is the shear modulus, and $ν_{12}$ and $ν_{21}$ are Poisson's ratios.

The damage variables $d_{f}$ , $d_{m}$ , and $d_{s}$ are derived from damage variables $d_{f}^{t}$ , $d_{f}^{c}$ , $d_{m}^{t}$ , and $d_{m}^{c}$ , corresponding to the four failure modes previously discussed, as follows:

\begin{matrix} d_{f} & = {\begin{matrix} d_{f}^{t} & if {\hat{σ}}_{11} \geq 0, \\ d_{f}^{c} & if {\hat{σ}}_{11} < 0, \end{matrix} \\ d_{m} & = {\begin{matrix} d_{m}^{t} & if {\hat{σ}}_{22} \geq 0, \\ d_{m}^{c} & if {\hat{σ}}_{22} < 0, \end{matrix} \\ d_{s} & ​ = 1 - (1 - d_{f}^{t}) (1 - d_{f}^{c}) (1 - d_{m}^{t}) (1 - d_{m}^{c}), \end{matrix}

${\hat{σ}}_{11}$ and ${\hat{σ}}_{22}$ are components of the effective stress tensor. The effective stress tensor is primarily used to evaluate damage initiation criteria; see Damage Initiation for Fiber-Reinforced Composites for a description of how the effective stress tensor is computed.

Evolution of Damage Variables for Each Mode

To alleviate mesh dependency during material softening, Abaqus introduces a characteristic length into the formulation, so that the constitutive law is expressed as a stress-displacement relation. The damage variable evolves such that the stress-displacement behaves as shown in Figure 1 in each of the four failure modes. The positive slope of the stress-displacement curve prior to damage initiation corresponds to linear elastic material behavior; the negative slope after damage initiation is achieved by evolution of the respective damage variables according to the equations shown below.

Equivalent displacement and stress for each of the four damage modes are defined as follows:

Fiber tension $({\hat{σ}}_{11} \geq 0)$ :

δ_{e q}^{f t} = L^{c} \sqrt{{⟨ ε_{11} ⟩}^{2} + α ε_{12}^{2}},

σ_{e q}^{f t} = \frac{⟨ σ_{11} ⟩ ⟨ ε_{11} ⟩ + α τ_{12} ε_{12}}{δ_{e q}^{f t} / L^{c}},

Fiber compression $({\hat{σ}}_{11} < 0)$ :

δ_{e q}^{f c} = L^{c} ⟨ - ε_{11} ⟩,

σ_{e q}^{f c} = \frac{⟨ - σ_{11} ⟩ ⟨ - ε_{11} ⟩}{δ_{e q}^{f c} / L^{c}} .

Matrix tension $({\hat{σ}}_{22} \geq 0)$ :

δ_{e q}^{m t} = L^{c} \sqrt{{⟨ ε_{22} ⟩}^{2} + ε_{12}^{2}},

σ_{e q}^{m t} = \frac{⟨ σ_{22} ⟩ ⟨ ε_{22} ⟩ + τ_{12} ε_{12}}{δ_{e q}^{m t} / L^{c}} .

Matrix compression $({\hat{σ}}_{22} < 0)$ :

δ_{e q}^{m c} = L^{c} \sqrt{{⟨ - ε_{22} ⟩}^{2} + ε_{12}^{2}},

σ_{e q}^{m c} = \frac{⟨ - σ_{22} ⟩ ⟨ - ε_{22} ⟩ + τ_{12} ε_{12}}{δ_{e q}^{m c} / L^{c}} .

The characteristic length, $L^{c}$ , is based on the element geometry and formulation: it is a typical length of a line across an element for a first-order element; it is half of the same typical length for a second-order element. For membranes and shells it is a characteristic length in the reference surface, computed as the square root of the area. Alternatively, this characteristic length could be defined as a function of the element topology and material orientation in user subroutine VUCHARLENGTH (see Defining the Characteristic Element Length at a Material Point in Abaqus/Explicit). The symbol $⟨ ⟩$ in the equations above represents the Macaulay bracket operator, which is defined for every $α \in ℜ$ as $⟨ α ⟩ = (α + | α |) / 2$ .

After damage initiation (i.e., $δ_{e q} \geq δ_{e q}^{0}$ ) for the behavior shown in Figure 1, the damage variable for a particular mode is given by the following expression

d = \frac{δ_{e q}^{f} (δ_{e q} - δ_{e q}^{0})}{δ_{e q} (δ_{e q}^{f} - δ_{e q}^{0})},

where $δ_{e q}^{0}$ is the initial equivalent displacement at which the initiation criterion for that mode was met and $δ_{e q}^{f}$ is the displacement at which the material is completely damaged in this failure mode. The above relation is presented graphically in Figure 2.

The values of $δ_{e q}^{0}$ for the various modes depend on the elastic stiffness and the strength parameters specified as part of the damage initiation definition (see Damage Initiation for Fiber-Reinforced Composites). For each failure mode you must specify the energy dissipated due to failure, $G^{c}$ , which corresponds to the area of the triangle OAC in Figure 3.

The values of $δ_{e q}^{f}$ for the various modes depend on the respective $G^{c}$ values.

Unloading from a partially damaged state, such as point B in Figure 3, occurs along a linear path toward the origin in the plot of equivalent stress versus equivalent displacement; this same path is followed back to point B upon reloading as shown in the figure.

Damage Evolution for the Ply Fabric Model

The ply fabric damage initiation criterion for bidirectional fabric-reinforced composites is described in Damage Initiation for Fiber-Reinforced Composites. This section discusses the post-damage initiation behavior for cases in which a damage evolution model is specified.

The elastic stress-strain relations are given by orthotropic damaged bilamina elasticity (see Defining Orthotropic Elasticity in Plane Stress with Different Moduli in Tension and Compression). Thereafter, the damaged elastic response of the plane stress orthotropic material in the local orthogonal coordinate is computed from

σ = C_{d} ε^{e l},

where $ε^{e l}$ is the elastic strain, and $C_{d}$ is the damaged elasticity matrix. The damaged elasticity matrix has the form

C_{d} = \frac{1}{D} [\begin{matrix} (1 - d_{1}) E_{1} & (1 - d_{1}) (1 - d_{2}) ν_{21} E_{1} & 0 \\ (1 - d_{1}) (1 - d_{2}) ν_{12} E_{2} & (1 - d_{2}) E_{2} & 0 \\ 0 & 0 & (1 - d_{12}) G_{12} D \end{matrix}],

where $D = 1 - (1 - d_{1}) (1 - d_{2}) ν_{12} ν_{21}$ , $d_{1}$ reflects the current state of fiber damage in the local 1-direction, $d_{2}$ reflects the current state of fiber damage in the local 2-direction, $d_{12}$ reflects the current state of matrix shear damage, $E_{1}$ is the Young's modulus in the local 1-direction, $E_{2}$ is the Young's modulus in the local 2-direction, $G_{12}$ is the shear modulus, and $ν_{12}$ and $ν_{21}$ are Poisson's ratios.

The model differentiates between tensile and compressive fiber failure modes by activating the corresponding damage variable depending on the stress state in the fiber directions. Thus:

d_{1} = d_{1 +} \frac{< σ_{11} >}{σ_{11}} + d_{1 -} \frac{< - σ_{11} >}{σ_{11}},

and

d_{2} = d_{2 +} \frac{< σ_{22} >}{σ_{22}} + d_{2 -} \frac{< - σ_{22} >}{σ_{22}} .

The symbol $⟨ ⟩$ in the equations above represents the Macaulay bracket operator. For the four damage variables $d_{α}$ , $(α \in {1 +, 1 -, 2 +, 2 -})$ are associated with fiber failure. An index $α$ is introduced to simplify notation and is used in subsequent discussions.

To incorporate different initial (undamaged) stiffness in tension and compression, the values of the elastic constants $E_{1}$ , $E_{2}$ , and $ν_{12}$ are interpolated between their tensile values ( $E_{1 +}$ , $E_{2 +}$ , and $ν_{12 +}$ ) or compressive values ( $E_{1 -}$ , $E_{2 -}$ , and $ν_{12 -}$ ) depending on the state of the material. You must define these elastic constants using bilamina elasticity (see Defining Orthotropic Elasticity in Plane Stress with Different Moduli in Tension and Compression).

Evolution of Damage Variables for Each Mode

For the four damage variables $d_{α}$ , $(α \in {1 +, 1 -, 2 +, 2 -})$ associated with fiber failure, the evolution of these damage variables is a function of the damage thresholds and the fracture energy per unit area under uniaxial tensile/compressive loading, $G_{f α}$ . The evolution of the damage variables is given by the equation:

d_{α} = 1 - \frac{1}{r_{α}} e^{- (A_{α} (r_{α} - 1))}, {\dot{d}}_{α} \geq 0,

where

A_{α} = \frac{2 g_{0}^{α} L^{c}}{G_{f α} - g_{0}^{α} L^{c}} .

Here, $L^{c}$ is the characteristic length of the element, $G_{f α}$ is the fracture energy per unit area under uniaxial tensile/compressive loading, and $g_{0}^{α}$ is the elastic energy density (that is. per unit volume) at the point of damage initiation:

g_{0}^{α} = \frac{X_{α}^{2}}{{2 E}_{α}},

where $X_{α}$ are the tensile/compressive strengths for uniaxial loading along the fiber directions. The formulation of the damage evolution law ensures that the damage variables are monotonically increasing quantities. It also ensures that the correct amount of energy is dissipated when the lamina is subjected to uniaxial loading conditions along the fiber directions.

The characteristic length, $L^{c}$ , is based on the element geometry and formulation. It is a typical length of a line across an element for a first-order element; it is half of the same typical length for a second-order element. For membrane and shell elements it is a characteristic length in the reference surface, computed as the square root of the area. Alternatively, you can define this characteristic length as a function of the element topology and material orientation in user subroutine VUCHARLENGTH (see Defining the Characteristic Element Length at a Material Point in Abaqus/Explicit).

It is assumed that the shear damage variable $d_{12}$ increases with the logarithm of $r_{12}$ until a maximum value $d_{12}^{\max}$ is reached. Thus:

d_{12} = \min (α_{12} \ln (r_{12}), d_{12}^{\max}),

where $α_{12} > 0$ , and $d_{12}^{\max} \leq 1$ are material properties.

Calibration Procedure

You can measure the elastic constants and the fiber tensile/compressive strengths, $X_{α}$ , from standard coupon tests in uniaxial tension/compression loading of $0^{°} / 90^{°}$ laminates. The calibration of damage evolution in the fiber failure modes is based on the fracture energy per unit area of the material, $G_{f α}$ , which can be measured experimentally.

The shear response is usually calibrated with a cyclic tensile test on a $\pm 45^{°}$ laminate, where the strains along the fiber directions can be neglected. Figure 4 shows the typical shear response of a fabric-reinforced composite. The unloading/reloading paths correspond to an idealization of the actual response, which usually exhibits hysteretic behavior. Figure 4 serves as a starting point for the discussion of a general calibration procedure for the parameters that enter the damage and plasticity equations.

The level of damage can be measured from the ratio of the unloading stiffness to the initial (undamaged) elastic stiffness. You can compute pairs of stress-damage values, $(σ_{12}, d_{12})$ , for each unloading curve. This data can be represented in the space of $d_{12}$ versus $\ln ({\hat{σ}}_{12})$ , where ${\hat{σ}}_{12} = σ_{12} / (1 - d_{12})$ . A linear fit of the data provides the values of $α_{12}$ (slope of the line) and $S$ (intersection with the horizontal axis), as shown in Figure 5. Sometimes the damage data show indication of a saturation value, which would be used to determine $d_{12}^{\max} \leq 1$ . Otherwise, you should use a value of $d_{12}^{\max} = 1$ .

Finally, for each unloading curve in Figure 4, the plastic strain $ε_{12}^{p l}$ at the onset of unloading is determined from the value of residual deformation in the unloaded state. The values of $({\hat{σ}}_{12}, ε_{12}^{p l})$ at the onset of unloading are then used to fit the Johnson-Cook parameters of the yield function or, alternatively, they can be used directly as tabular input of the hardening curve.

Maximum Degradation and Choice of Element Removal

You have control over how Abaqus treats elements with severe damage. By default, the upper bound to all damage variables at a material point is $d_{m a x} = 1.0$ . You can reduce this upper bound as discussed in Controlling Element Deletion and Maximum Degradation for Materials with Damage Evolution.

In Abaqus/Standard a material point is considered to fail once all failure modes reach $d_{m a x}$ . In Abaqus/Explicit a material point is assumed to fail when either of the damage variables associated with fiber failure modes (tensile or compressive) reaches $d_{m a x}$ . By default, an element is deleted from a mesh upon material failure. Details for element deletion driven by material failure are described in Material Failure and Element Deletion. The status of a material point and an element can be determined by requesting output variables STATUSMP and STATUS, respectively.

Alternatively, you can specify that an element should remain in the model even after all of the damage variables reach $d_{m a x}$ (see Controlling Element Deletion and Maximum Degradation for Materials with Damage Evolution). In this case, once all of the damage variables reach the maximum value, the stiffness, $C_{d}$ , remains constant (see the expression for $C_{d}$ earlier in this section).

For ply fabric damage evolution (available only in Abaqus/Explicit), you have more control over the criteria for material failure. You can consider a material point to fail if any one of the following criteria is met:

when either fiber damage variable reaches a specified maximum value $d_{1} = d_{m a x}$ or $d_{2} = d_{m a x}$ , or when both fiber damage variables reach the maximum specified value $d_{1} = d_{m a x}$ and $d_{2} = d_{m a x}$ ;
when the shear plastic strain reaches a specified maximum value, ${\bar{ε}}^{p l} = {\bar{ε}}_{\max}^{p l}$ ; or
when a deformation-based material failure criterion depending on the values of the maximum ( ${\hat{ε}}_{\max} > 0$ ) and minimum ( ${\hat{ε}}_{\min} < 0$ ) principal logarithmic strains that the material can sustain is met.

If you do not define conditions that trigger material failure for the ply fabric model, the material point fails when either fiber damage variable reaches a specified maximum value.

Viscous Regularization

Material models exhibiting softening behavior and stiffness degradation often lead to severe convergence difficulties in implicit analysis programs, such as Abaqus/Standard. You can overcome some of these convergence difficulties by using the viscous regularization scheme, which causes the tangent stiffness matrix of the softening material to be positive for sufficiently small time increments.

In this regularization scheme a viscous damage variable is defined by the evolution equation:

{\dot{d}}_{v} = \frac{1}{η} (d - d_{v}),

where $η$ is the viscosity coefficient representing the relaxation time of the viscous system and d is the damage variable evaluated in the inviscid backbone model. The damaged response of the viscous material is given as

σ = C_{d} ϵ,

where the damaged elasticity matrix, $C_{d}$ , is computed using viscous values of damage variables for each failure mode. Using viscous regularization with a small value of the viscosity parameter (small compared to the characteristic time increment) usually helps improve the rate of convergence of the model in the softening regime, without compromising results. The basic idea is that the solution of the viscous system relaxes to that of the inviscid case as $t / η \to \infty$ , where t represents time.

Viscous regularization is also available in Abaqus/Explicit for the Hashin damage model. Viscous regularization slows down the rate of increase of damage and leads to increased fracture energy with increasing deformation rates, which can be exploited as an effective method of modeling rate-dependent material behavior.

In Abaqus/Standard the approximate amount of energy associated with viscous regularization over the whole model or over an element set is available using output variable ALLCD.

You cannot use the viscous regularization scheme with damage evolution for the ply fabric damage initiation criterion.

Defining Viscous Regularization Coefficients

You can specify different values of viscous coefficients for different failure modes.

Applying a Single Viscous Coefficient in Abaqus/Standard

Alternatively, in Abaqus/Standard you can specify the viscous coefficients as part of a section controls definition. In this case the same viscous coefficient is applied to all failure modes. For more information, see Using Viscous Regularization with Cohesive Elements, Connector Elements, and Elements That Can Be Used with the Damage Evolution Models for Ductile Metals and Fiber-Reinforced Composites in Abaqus/Standard.

Material Damping

If stiffness proportional damping is specified in combination with the damage evolution law for fiber-reinforced materials, Abaqus calculates the damping stresses using the damaged elastic stiffness.

Elements

The damage evolution law for fiber-reinforced materials must be used with elements with a plane stress formulation, which include plane stress, shell, continuum shell, and membrane elements.

Output

In addition to the standard output identifiers available in Abaqus (Abaqus/Standard Output Variable Identifiers and Abaqus/Explicit Output Variable Identifiers ), the following variables relate specifically to damage evolution in the fiber-reinforced composite damage model:

STATUS: Status of the element (the status of an element is 1.0 if the element is active, 0.0 if the element is not).
STATUSMP: Status of each material point in the element (1.0 if a material point is active, 0.0 if it is not). Abaqus/Explicit only.
DAMAGEFT: Fiber tensile damage variable for the Hashin model.
DAMAGEFC: Fiber compressive damage variable for the Hashin model.
DAMAGEMT: Matrix tensile damage variable for the Hashin model.
DAMAGEMC: Matrix compressive damage variable for the Hashin model.
DAMAGESHR: Shear damage variable for the Hashin model and ply fabric model.
DAMAGEF1T: Fiber tensile damage variable in the local 1-direction for the ply fabric model.
DAMAGEF1C: Fiber compressive damage variable in the local 1-direction for the ply fabric model.
DAMAGEF2T: Fiber tensile damage variable in the local 2-direction for the ply fabric model.
DAMAGEF2C: Fiber compressive damage variable in the local 2-direction for the ply fabric model.
EDMDDEN: Energy dissipated per unit volume in the element by damage.
ELDMD: Total energy dissipated in the element by damage.
DMENER: Energy dissipated per unit volume by damage.
ALLDMD: Energy dissipated in the whole (or partial) model by damage.
ALLPD: Energy dissipated in the whole (or partial) model by shear plasticity in the ply fabric model.
ECDDEN: Energy per unit volume in the element that is associated with viscous regularization.
ELCD: Total energy in the element that is associated with viscous regularization.
CENER: Energy per unit volume that is associated with viscous regularization.
ALLCD: The approximate amount of energy over the whole model or over an element set that is associated with viscous regularization.

References

Lapczyk, I., and J. A. Hurtado, “Progressive Damage Modeling in Fiber-Reinforced Materials,” Composites Part A: Applied Science and Manufacturing, vol. 38, no. 11, pp. 2333–2341, 2007.
Sokolinsky, V. S., K. C. Indermuehle, and J. A. Hurtado, “Numerical Simulation of the Crushing Process of a Corrugated Composite Plate,” Composites Part A: Applied Science and Manufacturing, vol. 42, no. 9, pp. 1119–1126, 2011.