5.5 Reactive Damage Mechanics

Bond-Breaking Reaction

The reactive damage mechanics framework was first described in ¹. It is based on constrained reactive mixtures of solids (Section Constrained Reactive Mixture of Solids) and used to model damage in an elastic solid as a reaction that transforms intact (elastic) bonds into broken bonds,

\[ \begin{equation} \mathcal{E}^{i}\to\mathcal{E}^{b}\,.\label{eq:dmg-reaction} \end{equation} \]

Here, \(\mathcal{E}^{\alpha}\) is the material associated with bonds \(\alpha\) (\(\alpha=i\) for intact bonds and \(\alpha=b\) for broken bonds). The material is modeled as a constrained mixture of these two constituents \(\alpha\). Whereas intact bonds may store free energy, broken bond sustain none. This framework assumes that isothermal conditions prevail. Thus, any heat generated by the dissipative damage reaction must be radiated from the continuum to preserve a constant temperature. In an isothermal framework, the free energy density is also equal to the strain energy density.

The referential mass density of the solid mixture is \(\rho_{r}\) (mass of solid per volume in its referential, stress-free configuration), which remains constant throughout an analysis. The material associated with intact bonds has an apparent mass density \(\rho_{r}^{i}\) while that associated with broken bonds is \(\rho_{r}^{b},\)such that the mixture mass balance is satisfied by

\[ \begin{equation} \rho_{r}=\rho_{r}^{i}+\rho_{r}^{b}\,.\label{eq:dmg-mass-balance} \end{equation} \]

Strain Energy Density and Stress

Let the specific free energy stored in intact bonds be represented by \(\psi\left(\mathbf{F}\right)\); that of broken bonds is zero. Therefore, the free energy density of the mixture is

\[ \begin{equation} \Psi_{r}\left(\mathbf{F}\right)=\rho_{r}^{i}\psi\left(\mathbf{F}\right)\,.\label{eq:dmg-FED} \end{equation} \]

We may define the mass fraction \(w^{\alpha}\) of bond species \(\alpha\) as

\[ \begin{equation} w^{\alpha}=\frac{\rho_{r}^{\alpha}}{\rho_{r}}\,.\label{eq:dmg-mass-fraction} \end{equation} \]

Now, the mixture mass balance in eq.\eqref{eq:dmg-mass-balance} may be rewritten as \(\sum_{\alpha}w^{\alpha}=1\), or more specifically,

\[ \begin{equation} w^{i}+w^{b}=1\,.\label{eq:dmg-massfraction-balance} \end{equation} \]

We may also rewrite the mixture free energy density in eq.\eqref{eq:dmg-FED} as

\[ \begin{equation} \Psi_{r}\left(\mathbf{F}\right)=w^{i}\rho_{r}\psi\left(\mathbf{F}\right)=\left(1-w^{b}\right)\rho_{r}\psi\left(\mathbf{F}\right)\,,\label{eq:dmg-FED-redux} \end{equation} \]

where we have made use of eq.\eqref{eq:dmg-massfraction-balance}. The corresponding Cauchy stress may be evaluated using the standard hyperelasticity formula,

\[ \begin{equation} \boldsymbol{\sigma}=J^{-1}\frac{\partial\Psi_{r}}{\partial\mathbf{F}}\cdot\mathbf{F}^{T}=\left(1-w^{b}\right)\frac{\rho_{r}}{J}\frac{\partial\psi}{\partial\mathbf{F}}\cdot\mathbf{F}^{T}.\label{eq:dmg-stress} \end{equation} \]

These relation show that the free energy density and stress of a damaged material are scaled by the mass fraction \(w^{i}=1-w^{b}\) of remaining intact bonds. Comparing these formulas to those of classical damage mechanics ²³⁴⁵⁶⁷, it becomes immediately apparent that the classical damage variable \(D\) appearing in those theories is equivalent to the mass fraction \(w^{b}\) of broken bonds,

\[ \begin{equation} D\equiv w^{b}.\label{eq:dmg-variable} \end{equation} \]

To further clarify this equivalence, we may let \(\Psi_{0}\equiv\rho_{r}\psi\) represent the free energy density of an intact elastic solid, such that eq.\eqref{eq:dmg-FED-redux} may be rewritten as \(\Psi_{r}=\left(1-D\right)\Psi_{0}\). Similarly, eq.\eqref{eq:dmg-stress} may be rewritten as \(\boldsymbol{\sigma}=\left(1-D\right)\boldsymbol{\sigma}_{0}\), where \(\boldsymbol{\sigma}_{0}\) is the stress in the intact elastic solid, derived from the hyperelasticity relation \(\boldsymbol{\sigma}_{0}=J^{-1}\left(\partial\Psi_{0}/\partial\mathbf{F}\right)\cdot\mathbf{F}^{T}\).

For nearly-incompressible hyperelastic materials (Section Nearly-Incompressible Hyperelasticity), the strain energy density of the intact material has the form of eq.(2.6-53), thus \(\Psi_{0}\left(\mathbf{C}\right)=\tilde{\Psi}_{0}\left(\tilde{\mathbf{C}}\right)+U\left(J\right)\). In this case, we assume that the damage only affects the distortional part of the strain energy density \(\tilde{\Psi}_{0}\left(\tilde{\mathbf{C}}\right)\), consistent with the general framework advocated in ⁸. Thus, for uncoupled damage, we assume that \(\Psi_{r}\left(\mathbf{C}\right)=\left(1-D\right)\tilde{\Psi}_{0}\left(\tilde{\mathbf{C}}\right)+U\left(J\right)\). The resulting damage stress similarly takes the form \(\boldsymbol{\sigma}=\left(1-D\right)\dev\tilde{\boldsymbol{\sigma}}_{0}+p\mathbf{I}\), consistent with eq.(2.6-61), where \(\tilde{\boldsymbol{\sigma}}_{0}\) is evaluated from \(\tilde{\Psi}_{0}\left(\tilde{\mathbf{C}}\right)\) as given in eq.(2.6-62) and \(p\) is evaluated from \(U\left(J\right)\) as given in eq.(2.6-55).

When investigating the damage mechanics of tension-bearing fibrous materials, described in Section Tension-Bearing Fiber Materials, it is important to use the unconstrained version of the fiber and damage mechanics models, even when the fibers are embedded in a ground matrix with a nearly-incompressible response (uncoupled formulation). This is a necessary requirement since uncoupled fiber formulations are now understood to be non-physical. Nevertheless, for historical reasons, FEBio allows users to use uncoupled fiber formulations in an uncoupled damage material.

Damage Criterion

At each material point \(\mathbf{X}\) in the continuum, damage occurs when a scalar damage (or failure) measure \(\Xi\left(\mathbf{F}\right)\) achieves a critical value \(\Xi_{m}\) over the loading history,

\[ \begin{equation} \Xi_{m}\left(\mathbf{X}\right)=\max_{-\infty<s\le t}\Xi\left(\mathbf{F}\left(\mathbf{X},s\right)\right)\,.\label{eq:dmg-critical-measure} \end{equation} \]

The scalar damage measure \(\Xi\left(\mathbf{F}\right)\) must be invariant to orthogonal transformations \(\mathbf{Q}\) that preserve material symmetry, or else the damage formulation would not be observer-independent. For example, for isotropic materials, \(\Xi\left(\mathbf{F}\right)\) must be an isotropic function of the deformation, in which case it should be expressed as \(\Xi\left(\mathbf{U}\right)\) or \(\Xi\left(\mathbf{E}\right)\), where \(\mathbf{U}\) is the right stretch tensor in the polar decomposition \(\mathbf{F}=\mathbf{R}\cdot\mathbf{U}\) of the deformation gradient, and

\[ \begin{equation} \mathbf{E}=\frac{1}{2}\left(\mathbf{F}^{T}\cdot\mathbf{F}-\mathbf{I}\right)=\frac{1}{2}\left(\mathbf{U}^{2}-\mathbf{I}\right)\label{eq:Lagrange-strain-tensor} \end{equation} \]

is the Green-Lagrange strain tensor. It follows that

\[ \left(\frac{\partial\mathbf{E}}{\partial\mathbf{U}}\right)_{ijmn}=\frac{1}{2}\left(\frac{1}{2}\left(\delta_{im}U_{nj}+\delta_{in}U_{mj}\right)+\frac{1}{2}\left(U_{im}\delta_{jn}+U_{in}\delta_{jm}\right)\right)\,. \]

For anisotropic materials where \(\mathbf{a}\) is the unit normal to a symmetry plane and \(\mathbf{A}=\mathbf{a}\otimes\mathbf{a}\), the damage measure \(\Xi\) must satisfy

\[ \begin{equation} \Xi\left(\mathbf{U},\mathbf{A}\right)=\Xi\left(\mathbf{Q}\cdot\mathbf{U}\cdot\mathbf{Q}^{T},\mathbf{Q}\cdot\mathbf{A}\cdot\mathbf{Q}^{T}\right)\,,\label{eq:dmg-criterion-invariance} \end{equation} \]

for transformations \(\mathbf{Q}\) that satisfy \(\mathbf{Q}\cdot\mathbf{A}\cdot\mathbf{Q}^{T}=\mathbf{A}\) (or \(\mathbf{Q}\cdot\mathbf{a}=\mathbf{a}\)). We may replace \(\mathbf{U}\) with \(\mathbf{E}\) in the above expression.

We assume that the amount of damage (the fraction \(D\) of broken bonds) is given by the function of state

\[ \begin{equation} D=F\left(\Xi_{m}\right)\,,\label{eq:dmg-CDF} \end{equation} \]

where \(0\le F\left(\Xi_{m}\right)\le1\). As shown in ¹, the Clausius-Duhem inequality imposes the constraint that \(F\left(\Xi_{m}\right)\) must be a monotonically increasing function of its argument. Therefore, we may understand \(F\) to represent a cumulative density function (CDF), whose derivative \(f\left(\Xi_{m}\right)=F^{\prime}\left(\Xi_{m}\right)\) is a probability distribution function (PDF) that describes the probability of bonds breaking at the specific threshold \(\Xi_{m}\).

Reaction Kinetics and Thermodynamics

The axiom of mass balance in a reactive constrained mixture reduces to

\[ \begin{equation} \dot{\rho}_{r}^{\alpha}=\hat{\rho}_{r}^{\alpha}\,,\label{eq:dmg-mass-balance-alpha} \end{equation} \]

where \(\dot{\rho}_{r}^{\alpha}\) is the material time derivative of \(\rho_{r}^{\alpha}\) and \(\hat{\rho}_{r}^{\alpha}\) is a function of state representing the referential mass supply density to constituent \(\alpha\) due to reactions with all other constituents. In the damage framework, the above relations show that

\[ \begin{equation} \begin{aligned}\rho_{r}^{i} & =\rho_{r}\left(1-F\left(\Xi_{m}\right)\right)\\ \rho_{r}^{b} & =\rho_{r}F\left(\Xi_{m}\right) \end{aligned} \,.\label{eq:dmg-mass-densities} \end{equation} \]

Substituting these expressions into eq.\eqref{eq:dmg-mass-balance-alpha} shows that the referential mass density supplies are given by

\[ \begin{equation} \begin{aligned}\hat{\rho}_{r}^{i} & =-\rho_{r}\dot{F}\left(\Xi_{m}\right)\\ \hat{\rho}_{r}^{b} & =\rho_{r}\dot{F}\left(\Xi_{m}\right) \end{aligned} \,,\label{eq:dmg-mass-supplies} \end{equation} \]

where

\[ \begin{equation} \dot{F}\left(\Xi_{m}\right)=\begin{cases} \left.f\left(\Xi\right)\dot{\Xi}\right|_{\Xi_{m}} & \text{advancing damage}\\ 0 & \text{otherwise} \end{cases}\,.\label{eq:dmg-CDF-mtd} \end{equation} \]

In these expressions, the damage is advancing when \(\Xi_{m}\) increases over two consecutive time points. In this expression for \(\dot{F}\) we need to evaluate

\[ \begin{equation} \dot{\Xi}\left(\mathbf{U}\right)=\frac{\partial\Xi}{\partial\mathbf{U}}:\dot{\mathbf{U}}=\mathbf{N}:\dot{\mathbf{U}}\,,\label{eq:dmg-Xsi-dot} \end{equation} \]

where we defined

\[ \begin{equation} \mathbf{N}\equiv\frac{\partial\Xi}{\partial\mathbf{U}}\,\label{eq:dmg-surface-normal} \end{equation} \]

to represent the tensorial normal to the damage hypersurface, which needs to be evaluated at \(\Xi_{m}\).

In this isothermal damage framework it can be shown from the energy balance that a heat supply density \(\rho_{r}r\) must radiate the bond-breaking energy out of the continuum to maintain isothermal conditions, where

\[ \begin{equation} \rho_{r}r=\hat{\rho}_{r}^{i}\psi\left(\mathbf{F}\right)=-\rho_{r}\dot{F}\left(\Xi_{m}\right)\psi\left(\mathbf{F}\right)\,.\label{eq:dmg-heat-supply} \end{equation} \]

Since \(F\) is a monotonically increasing function of \(\Xi_{m}\), its material time derivative \(\dot{F}\) is always positive when the damage is increasing, and zero otherwise as per eq.\eqref{eq:dmg-Xsi-dot}. Since the specific strain energy \(\psi\left(\mathbf{F}\right)\) is always positive, it follows that the specific heat supply \(r=-\dot{F}\left(\Xi_{m}\right)\psi\left(\mathbf{F}\right)\) in eq.\eqref{eq:dmg-heat-supply} is negative or zero, consistent with the expectation that heat needs to leave the continuum to maintain isothermal conditions.

Constitutive Models for Damage and Yield Criteria

Constitutive models for the damage or yield measure \(\Xi\left(\mathbf{U}\right)\) may be derived from energy- or stress-based potentials or (less commonly) strain measures. A summary of constitutive models for damage or yield criteria currently implemented in FEBio is presented below.

Strain Energy Density

It may be assumed that damage occurs when the strain energy density \(\Psi_{r}\) achieves a certain threshold \(\Psi_{m}\). In that case, the constitutive model for the damage measure is

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\Psi_{r}\left(\mathbf{U}\right)\,.\label{eq:dmg-model-SED} \end{equation} \]

This damage measure is valid for isotropic or anisotropic materials, since the strain energy density \(\Psi_{r}\) must satisfy the frame-invariance of eq.(5.5-11) by construction. The resulting damage surface normal is

\[ \begin{equation} \mathbf{N}=\frac{\partial\Psi_{r}}{\partial\mathbf{U}}=\frac{1}{2}\left(\mathbf{S}\cdot\mathbf{U}+\mathbf{U}\cdot\mathbf{S}\right)\,,\label{eq:dmg-model-SED-N} \end{equation} \]

where \(\mathbf{S}=\partial\Psi_{r}/\partial\mathbf{E}=J\mathbf{F}^{-1}\cdot\boldsymbol{\sigma}\cdot\mathbf{F}^{-T}\) is the second Piola-Kirchhoff stress associated with the material.

Simo Damage Criterion

Simo ⁹⁷ proposed a damage criterion related to the strain energy density,

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\sqrt{2\Psi_{r}\left(\mathbf{U}\right)}\,.\label{eq:dmg-model-Simo} \end{equation} \]

This damage measure is valid for isotropic or anisotropic materials, since the strain energy density \(\Psi_{r}\) must satisfy the frame-invariance of eq.(5.5-11) by construction. Its resulting damage surface normal is

\[ \begin{equation} \mathbf{N}=\frac{1}{\sqrt{2\Psi_{r}}}\frac{1}{2}\left(\mathbf{S}\cdot\mathbf{U}+\mathbf{U}\cdot\mathbf{S}\right)\,,\label{eq:dmg-model-Simo-N} \end{equation} \]

where we have used the result of eq.\eqref{eq:dmg-model-SED-N}. The normal reduces to the null tensor in the limit when \(\Psi_{r}\) and \(\boldsymbol{\sigma}\) tend to zero.

Specific Strain Energy

It may be assumed that damage occurs when the specific strain energy \(\Psi_{r}/\rho_{r}\) achieves a certain threshold \(\psi_{m}\). In that case, the constitutive model for the damage measure is

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\frac{1}{\rho_{r}}\Psi_{r}\left(\mathbf{U}\right)\,.\label{eq:dmg-model-SSE} \end{equation} \]

This damage measure is valid for isotropic or anisotropic materials, since the strain energy density \(\Psi_{r}\) must satisfy the frame-invariance of eq.(5.5-11) by construction. The resulting damage surface normal is

\[ \begin{equation} \mathbf{N}=\frac{1}{\rho_{r}}\frac{\partial\Psi_{r}}{\partial\mathbf{U}}=\frac{1}{2\rho_{r}}\left(\mathbf{S}\cdot\mathbf{U}+\mathbf{U}\cdot\mathbf{S}\right)\,,\label{eq:dmg-model-SSE-N} \end{equation} \]

where we have used the result of eq.\eqref{eq:dmg-model-SED-N}.

Von Mises Stress

For this criterion it is assumed that damage or yield is initiated by increases in the von Mises (or effective) stress, \(\sigma_{Y}\). Thus,

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\sigma_{Y}\left(\mathbf{U}\right)=\sqrt{\frac{3}{2}\dev\boldsymbol{\sigma}:\dev\boldsymbol{\sigma}}\,,\label{eq:dmg-model-VMS} \end{equation} \]

where \(\dev\boldsymbol{\sigma}\) is the deviatoric part of \(\boldsymbol{\sigma}\). To evaluate the damage surface normal in this case, we must use the chain rule,

\[ \begin{equation} \mathbf{N}=\frac{\partial\Xi}{\partial\boldsymbol{\sigma}}:\frac{\partial\boldsymbol{\sigma}}{\partial\mathbf{F}}:\frac{\partial\mathbf{F}}{\partial\mathbf{U}}\,.\label{eq:dmg-normal-chain-rule} \end{equation} \]

From the hyperelasticity relation in eq.\eqref{eq:dmg-stress} it can be shown that

\[ \begin{equation} \frac{\partial\boldsymbol{\sigma}}{\partial\mathbf{F}}=\left(\boldsymbol{\mathcal{C}}+\mathbf{I}\oslash\boldsymbol{\sigma}+\boldsymbol{\sigma}\obslash\mathbf{I}-\boldsymbol{\sigma}\otimes\mathbf{I}\right)\cdot\mathbf{F}^{-T}\,,\label{eq:dmg-stress-tangent} \end{equation} \]

where \(\boldsymbol{\mathcal{C}}\) is the fourth-order spatial elasticity tensor associated with the strain energy density \(\Psi_{r}\). Then, it can be shown that

\[ \begin{equation} \mathbf{N}=\frac{1}{2}\mathbf{R}^{T}\cdot\mathbf{M}\cdot\mathbf{R}\cdot\mathbf{U}^{-1}+\frac{1}{2}\mathbf{U}^{-1}\cdot\mathbf{R}^{T}\cdot\mathbf{M}^{T}\cdot\mathbf{R}\,,\label{eq:dmg-N-redux} \end{equation} \]

where

\[ \begin{equation} \mathbf{M}=\frac{\partial\Xi}{\partial\boldsymbol{\sigma}}:\boldsymbol{\mathcal{C}}+2\frac{\partial\Xi}{\partial\boldsymbol{\sigma}}\cdot\boldsymbol{\sigma}-\left(\frac{\partial\Xi}{\partial\boldsymbol{\sigma}}:\boldsymbol{\sigma}\right)\mathbf{I}\,.\label{eq:dmg-N-M} \end{equation} \]

From the von Mises criterion in eq.\eqref{eq:dmg-model-VMS}, it can be shown that

\[ \begin{equation} \frac{\partial\Xi}{\partial\boldsymbol{\sigma}}=\frac{3}{2\sigma_{Y}}\dev\boldsymbol{\sigma}\,.\label{eq:dmg-VMS-tangent} \end{equation} \]

Substituting Eqs.\eqref{eq:dmg-N-M}-\eqref{eq:dmg-VMS-tangent} into eq.\eqref{eq:dmg-N-redux} provides the surface normal \(\mathbf{N}\).

Maximum Normal Stress

For this criterion, the damage measure takes the form

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\sigma_{1}\,,\label{eq:dmg-model-MNS} \end{equation} \]

where \(\sigma_{1}\) is the maximum principal stress (under the assumption that the three principal stresses of \(\boldsymbol{\sigma}\) are ordered such that \(\sigma_{1}\ge\sigma_{2}\ge\sigma_{3}\)). To evaluate the damage surface normal, we may use Eqs.\eqref{eq:dmg-N-redux}-\eqref{eq:dmg-N-M} where

\[ \begin{equation} \frac{\partial\Xi}{\partial\boldsymbol{\sigma}}=\frac{\partial\sigma_{1}}{\partial\boldsymbol{\sigma}}=\mathbf{n}_{1}\otimes\mathbf{n}_{1}\,.\label{eq:dmg-MNS-tangent} \end{equation} \]

Here, \(\mathbf{n}_{1}\) is a unit vector along the principal direction of maximum normal stress (the eigenvector of \(\boldsymbol{\sigma}\) corresponding to the eigenvalue \(\sigma_{1}\)).

Maximum Shear Stress

For this criterion, the damage measure takes the form

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=\frac{\sigma_{1}-\sigma_{3}}{2}\,,\label{eq:dmg-model-MSS} \end{equation} \]

where \(\sigma_{1}\) and \(\sigma_{3}\) are the maximum and minimum principal normal stresses (under the assumption that the three principal stresses of \(\boldsymbol{\sigma}\) are ordered such that \(\sigma_{1}\ge\sigma_{2}\ge\sigma_{3}\)). To evaluate the damage surface normal, we may use Eqs.\eqref{eq:dmg-N-redux}-\eqref{eq:dmg-N-M} where

\[ \begin{equation} \frac{\partial\Xi}{\partial\boldsymbol{\sigma}}=\frac{1}{2}\left(\mathbf{n}_{1}\otimes\mathbf{n}_{1}-\mathbf{n}_{3}\otimes\mathbf{n}_{3}\right)\,.\label{eq:dmg-MSS-tangent} \end{equation} \]

Here, \(\mathbf{n}_{1}\) and \(\mathbf{n}_{3}\) are unit vectors along the principal directions of maximum and minimum normal stress, respectively.

Drucker Shear Stress

This criterion is based on the yield criterion for plasticity introduced in ¹⁰. Its damage (or yield) measure takes the form

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=k=\left(J_{2}^{3}-cJ_{3}^{2}\right)^{1/6}\,,\label{eq:dmg-model-Drucker} \end{equation} \]

where \(J_{2}=\frac{1}{2}\dev\boldsymbol{\sigma}:\dev\boldsymbol{\sigma}\), \(J_{3}=\det\left(\dev\boldsymbol{\sigma}\right)\), \(k\) is the yield limit in simple shear and \(c\) is a user-specified non-dimensional material constant which must lie in the range \(-\frac{27}{8}\le c\le\frac{9}{4}\). To better understand the meaning of \(k\), consider uniaxial loading of a bar which yields at the normal stress value of \(\sigma_{y}\). In this case,

\[ \begin{equation} k=\frac{\sigma_{y}}{\sqrt{3}}\left(1-\frac{4}{27}c\right)^{1/6}\quad\frac{\sigma_{y}}{\sqrt{3}}\left(\frac{2}{3}\right)^{1/6}\le k\le\frac{\sigma_{y}}{\sqrt{3}}\left(\frac{3}{2}\right)^{1/6}\,.\label{eq:dmg-Drucker-range} \end{equation} \]

In the special case when \(c=0\) the Drucker criterion reduces to the von Mises criterion, with \(k=\sigma_{y}/\sqrt{3}\). To evaluate the damage or yield surface normal, we may use Eqs.\eqref{eq:dmg-N-redux}-\eqref{eq:dmg-N-M} where

\[ \begin{equation} \frac{\partial\Xi}{\partial\boldsymbol{\sigma}}=\frac{1}{k^{5}}\left(\frac{J_{2}^{2}}{2}\dev\boldsymbol{\sigma}-\frac{cJ_{3}^{2}}{3}\dev\left(\dev\boldsymbol{\sigma}\right)^{-1}\right)\,.\label{eq:dmg-Drucker-tangent} \end{equation} \]

Maximum Normal Lagrange Strain

The Lagrange strain tensor \(\mathbf{E}\) is related to \(\mathbf{F}\) via eq.\eqref{eq:Lagrange-strain-tensor}. Its maximum principal value is denoted by \(E_{1}\). For this criterion, we let

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=E_{1}\,.\label{eq:dmg-model-MNLS} \end{equation} \]

Then, the damage surface normal is given by

\[ \begin{equation} \mathbf{N}=\frac{\partial E_{1}}{\partial\mathbf{U}}=\sqrt{1+2E_{1}}\mathbf{n}_{1}\otimes\mathbf{n}_{1}\,,\label{eq:dmg-MNLS-N} \end{equation} \]

where \(\mathbf{n}_{1}\) is a unit vector along the principal direction of normal strain.

Note that the maximum normal Lagrange strain is a kinematic measure, thus it does not represent an intrinsic material property. Though FEBio allows users to specify this criterion, it is strongly recommended to employ stress- or strain-energy-density-based failure or yield criteria in practice.

Octahedral Lagrange Strain

The octahedral Lagrange strain is given by

\[ \begin{equation} e\left(\mathbf{U}\right)=\sqrt{\frac{2}{3}\dev\mathbf{E}:\dev\mathbf{E}}\,,\label{eq:octahedral-Lagrange-strain} \end{equation} \]

where \(\dev\mathbf{E}\) is the deviatoric part of the Lagrange strain tensor,

\[ \begin{equation} \dev\mathbf{E}=\mathbf{E}-\frac{1}{3}\tr\left(\mathbf{E}\right)\mathbf{I}\,,\label{eq:deviatoric-Lagrange-strain} \end{equation} \]

and \(\mathbf{E}\) is given in eq.\eqref{eq:Lagrange-strain-tensor}. For this damage measure we let

\[ \begin{equation} \Xi\left(\mathbf{U}\right)=e\left(\mathbf{U}\right)\,.\label{eq:dmg-model-OLS} \end{equation} \]

The damage surface normal can be evaluated from

\[ \begin{equation} \mathbf{N}=\frac{\partial e}{\partial\mathbf{E}}:\frac{\partial\mathbf{E}}{\partial\mathbf{U}}=\frac{1}{3e}\left(\dev\mathbf{E}\cdot\mathbf{U}+\mathbf{U}\cdot\dev\mathbf{E}\right)\,,\label{eq:dmg-OLS-N} \end{equation} \]

where we used

\[ \begin{equation} \frac{\partial e}{\partial\mathbf{E}}=\frac{2}{3e}\dev\mathbf{E}\label{eq:dmg-e-tangent} \end{equation} \]

and

\[ \begin{equation} \frac{\partial\mathbf{E}}{\partial\mathbf{U}}=\frac{1}{2}\frac{\partial\mathbf{U}^{2}}{\partial\mathbf{U}}=\frac{1}{2}\left(\mathbf{I}\odot\mathbf{U}+\mathbf{U}\odot\mathbf{I}\right)\,.\label{eq:dmg-E-U-tangent} \end{equation} \]

Note that the octahedral Lagrange strain is a kinematic measure, thus it does not represent an intrinsic material property. Though FEBio allows users to specify this criterion, it is strongly recommended to employ stress- or strain-energy-density-based failure or yield criteria in practice.

---

1. Nims, Robert J; Durney, Krista M; Cigan, Alexander D; Duss{\'e}aux, Antoine; Hung, Clark T; Ateshian, Gerard A. "Continuum theory of fibrous tissue damage mechanics using bond kinetics: application to cartilage tissue engineering." Interface Focus, vol. 6, pp. 20150063 (2016). ↩↩

2. Kachanov, Lazar M. "Rupture time under creep conditions." International Journal of Fracture (1958). ↩

3. Rabotnov, Yu N. "Elements of hereditary solid mechanics." MIT Publishers, Moscow (1980). ↩

4. Chaboche, Jean-Louis. "Continuous damage mechanics---a tool to describe phenomena before crack initiation." Nuclear Engineering and Design, vol. 64, pp. 233--247 (1981). ↩

5. Lemaitre, Jean. "How to use damage mechanics." Nuclear engineering and design, vol. 80, pp. 233--245 (1984). ↩

6. Lemaitre, Jean. "A continuous damage mechanics model for ductile fracture." J. Eng. Mater. Technol. (1985). ↩

7. Simo, Juan C; Ju, JW. "Strain-and stress-based continuum damage models---I. Formulation." International journal of solids and structures, vol. 23, pp. 821--840 (1987). ↩↩

8. Holzapfel, Gerhard A. "Nonlinear solid mechanics: a continuum approach for engineering." Wiley (2000). ↩

9. Simo, JC. "On a fully three-dimensional finite-strain viscoelastic damage model: formulation and computational aspects." Computer methods in applied mechanics and engineering, vol. 60, pp. 153--173 (1987). ↩

10. Drucker, Daniel Charles. "Relation of experiments to mathematical theories of plasticity." Journal of Applied Mechanics (1949). ↩