11 L09 — Damage Mechanics
Continuum Damage, Fracture, and Coupled Problems
11.1 Motivation: Degradation of Stiffness
Unlike plasticity (permanent strain), damage describes the progressive degradation of material stiffness due to microcrack growth and void coalescence.
Observable signatures:
- Reduction of the elastic modulus
- Stress softening (load-bearing capacity decreases with deformation)
- Final fracture / failure
Continuum Damage Mechanics (CDM) describes this via internal damage variables without resolving individual cracks.
11.2 Isotropic Damage Model
A scalar damage variable \(D \in [0, 1]\):
- \(D = 0\): virgin (undamaged) material
- \(D = 1\): fully damaged (no load-carrying capacity)
Effective stress concept (Lemaitre): \[ \tilde{\boldsymbol{\sigma}} = \frac{\boldsymbol{\sigma}}{1-D} \]
Constitutive law with damage: \[ \boldsymbol{\sigma} = (1-D)\,\mathbb{C}^e:\boldsymbol{\varepsilon} \]
11.3 Strain Equivalence Principle
Strain equivalence (Lemaitre): the strain in the damaged material under nominal stress \(\boldsymbol{\sigma}\) equals the strain in the undamaged material under effective stress \(\tilde{\boldsymbol{\sigma}}\).
This allows re-use of undamaged material models by simply replacing \(\boldsymbol{\sigma} \to \tilde{\boldsymbol{\sigma}}\).
Energy equivalence (Cordebois & Sidoroff): alternative that gives different predictions for shear.
11.4 Damage Evolution Law
Thermodynamic driving force (energy release rate \(Y\)): \[ Y = -\frac{\partial\Psi}{\partial D} = \frac{1}{2}\boldsymbol{\varepsilon}:\mathbb{C}^e:\boldsymbol{\varepsilon} \]
Damage criterion and evolution: \[ g(Y, D) = Y - r(D) \leq 0, \qquad \dot{D} = \dot{\mu}\frac{\partial g}{\partial Y}, \]
where \(r(D)\) is the damage threshold function (analogous to yield stress hardening).
11.5 Coupling Damage with Plasticity
Combined elasto-plastic-damage model free energy: \[ \Psi = (1-D)\hat{\Psi}^e(\boldsymbol{\varepsilon}^e) + \Psi^p(\bar{\varepsilon}^p) \]
The stress: \[ \boldsymbol{\sigma} = (1-D)\mathbb{C}^e:\boldsymbol{\varepsilon}^e \]
Both yield surface and damage criterion must be checked at each step. Order of priority depends on the coupling model chosen.
11.6 Softening and Regularization
Softening (negative tangent modulus) triggers:
- Loss of ellipticity of governing equations
- Mesh-dependent FEM solutions (localization into a band of one element width)
Regularization methods:
- Nonlocal damage: replace \(Y\) by a weighted spatial average \(\bar{Y}\)
- Gradient damage: add \(-c\,\Delta D\) to the damage driving force
- Phase-field fracture: variational approach with intrinsic length scale
11.7 Phase-Field Approach to Fracture
Introduce a smooth phase-field \(\phi \in [0,1]\) interpolating between intact (\(\phi=0\)) and fractured (\(\phi=1\)):
\[ \mathcal{E}(\mathbf{u}, \phi) = \int_\Omega \left[ g(\phi)\Psi^e_+(\boldsymbol{\varepsilon}) + \Psi^e_-(\boldsymbol{\varepsilon}) + \frac{G_c}{2}\left(\frac{\phi^2}{l_0} + l_0|\nabla\phi|^2\right) \right]dV \]
where \(g(\phi) = (1-\phi)^2\) is the degradation function and \(l_0\) the regularization length.
Advantages: no explicit crack tracking, natural branching/merging.
11.8 Thermo-Mechanical Coupling
Temperature \(\theta\) affects yield stress and elastic moduli: \[ \sigma_y = \sigma_y(\bar{\varepsilon}^p, \theta), \qquad \boldsymbol{\varepsilon}^\text{th} = \alpha(\theta - \theta_\text{ref})\mathbf{I}. \]
Energy equation with heat generated by plastic dissipation (Taylor-Quinney coefficient \(\beta \approx 0.9\)): \[ \rho c_p\dot{\theta} = \beta\,\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}^p + \nabla\cdot(k\nabla\theta) \]
Strong coupling requires staggered or monolithic FEM solvers.
11.9 Pressure-Dependent Yield Criteria (Recap)
Materials like soils, rocks, and concrete exhibit plastic behavior that depends strongly on hydrostatic pressure, unlike metals (which follow von Mises, pressure-independent).
Why Pressure Dependence Matters: - Frictional effects: higher normal stress increases shear strength - Dilatational effects: volume changes coupled to plastic flow - Confinement: material becomes stronger under compression
All pressure-dependent models consist of: 1. Additive strain decomposition: \(\boldsymbol{\varepsilon} = \boldsymbol{\varepsilon}^e + \boldsymbol{\varepsilon}^p\) 2. Hypoelastic law: \(\dot{\boldsymbol{\sigma}} = \mathbb{C}^e : \dot{\boldsymbol{\varepsilon}}^e\) 3. Yield function: \(f(\boldsymbol{\sigma}, q) \leq 0\) (depends on hydrostatic pressure \(p\)) 4. Flow rule: \(\dot{\boldsymbol{\varepsilon}}^p = \dot{\gamma} \frac{\partial g(\boldsymbol{\sigma}, q)}{\partial\boldsymbol{\sigma}}\) (often non-associative: \(g \neq f\)) 5. Hardening law: \(\dot{q} = \dot{\gamma} h(\boldsymbol{\sigma}, q)\) (e.g., cohesion increases with plastic strain) 6. Kuhn-Tucker conditions: \(\dot{\gamma} \geq 0\), \(f \leq 0\), \(\dot{\gamma}f = 0\)
11.10 Mohr-Coulomb Criterion
The Mohr-Coulomb criterion is the classical model for frictional materials (soils, rocks, concrete).
Physical Interpretation: Yielding occurs when shear stress on any plane reaches: \[\tau = c - \sigma_n\tan\phi\]
where: - \(\tau\) = shear stress on the failing plane - \(\sigma_n\) = normal stress on the failing plane - \(c\) = cohesion (material strength when \(\sigma_n = 0\)) - \(\phi\) = internal friction angle
This gives a linear failure envelope in \((\sigma_n, \tau)\) space.
Yield Function in Terms of Principal Stresses:
For principal stresses \(\sigma_1 \geq \sigma_2 \geq \sigma_3\): \[\Phi(\boldsymbol{\sigma}, c) = (\sigma_1 - \sigma_3) + (\sigma_1 + \sigma_3)\sin\phi - 2c\cos\phi = 0\]
Using invariants (with \(p\) = mean pressure, \(\theta\) = Lode angle, \(J_2, J_3\) = stress invariants): \[\Phi = \left(\cos\theta - \frac{\sin\theta\sin\phi}{\sqrt{3}}\right)\sqrt{J_2} + p\sin\phi - c\cos\phi = 0\]
Key Properties: - Non-smooth yield surface: hexagonal pyramid in principal stress space - Apex at \(p = c\cot\phi\) (tensile limit) - Generally non-associative flow rule (dilation angle \(\psi \neq \phi\)) - Hardening: cohesion becomes a function of accumulated plastic strain: \(c = c(\bar{\varepsilon}^p)\)
11.11 Tresca Criterion
The Tresca yield criterion (pressure-insensitive version of M-C) states that plasticity begins when the maximum shear stress reaches a critical value:
\[\tau_{\max} = \frac{1}{2}(\sigma_{\max} - \sigma_{\min})\]
Yield Function: \[\Phi(\boldsymbol{\sigma}) = (\sigma_1 - \sigma_3) - \sigma_y(\bar{\varepsilon}^p) = 0\]
where \(\sigma_y(\bar{\varepsilon}^p)\) is the uniaxial yield stress (function of hardening variable \(\bar{\varepsilon}^p\)).
Yield Surface: A hexagonal prism in principal stress space (axis = hydrostatic line).
Multisurface Representation: Six yield surfaces based on principal stress pairs: \[\Phi_1 = \sigma_1 - \sigma_3 - \sigma_y, \quad \Phi_2 = \sigma_2 - \sigma_3 - \sigma_y, \quad \Phi_3 = \sigma_1 - \sigma_2 - \sigma_y, \ldots\]
Only one or two surfaces are active at a given stress state; the algorithm must identify the active subset.
Return Mapping Procedure:
For a given principal stress ordering (\(\sigma_1 \geq \sigma_2 \geq \sigma_3\)), three cases arise:
- Main plane flow (\(\Phi_1\) only): one multiplier \(\Delta\gamma\)
- Updates: \(s_1 = s_{trial,1} - 2G\Delta\gamma\), \(s_3 = s_{trial,3} + 2G\Delta\gamma\), \(s_2\) unchanged
- Consistency: \(s_{trial,1} - s_{trial,3} - 4G\Delta\gamma - \sigma_y(\bar{\varepsilon}_{p,n} + \Delta\gamma) = 0\)
- Corner flow (two surfaces active, e.g., \(\Phi_1\) and \(\Phi_6\)): two multipliers \(\Delta\gamma_a, \Delta\gamma_b\)
- More complex stress updates and consistency conditions (system of 2 equations)
- Another corner (\(\Phi_1\) and \(\Phi_2\)): similar to case 2 but different geometry
The algorithm must determine which case applies and solve accordingly (often with iteration).
11.12 Drucker-Prager Criterion
The Drucker-Prager (DP) criterion is a smooth approximation to Mohr-Coulomb, incorporating pressure into a von Mises-like framework.
Yield Function: \[f_\text{DP}(\boldsymbol{\sigma}, c) = \sqrt{J_2(\mathbf{s}(\boldsymbol{\sigma}))} + \eta\, p(\boldsymbol{\sigma}) - \xi\, c = 0\]
where: - \(\sqrt{J_2}\) = equivalent deviatoric stress (like von Mises) - \(p\) = hydrostatic pressure - \(\eta, \xi\) = material parameters chosen to fit Mohr-Coulomb - \(c = c(\bar{\varepsilon}^p)\) = cohesion (hardening variable)
Yield Surface: A circular cone in principal stress space (smooth, isotropic about hydrostatic axis).
Advantages over M-C: - Smooth (no corners) → simpler numerical implementation - No need for multisurface logic - Single elastic region identification
Flow Rule (Non-Associative):
The flow potential is: \[\Psi(\boldsymbol{\sigma}, c) = \sqrt{J_2(\mathbf{s})} + \bar{\eta} p\]
Flow vector: \[\mathbf{N} = \frac{\partial\Psi}{\partial\boldsymbol{\sigma}} = \frac{\mathbf{s}}{2\sqrt{J_2}} + \frac{\bar{\eta}}{3}\mathbf{I}\]
where the second term represents volumetric (dilatational) flow. The dilatancy parameter \(\bar{\eta}\) can differ from \(\eta\) in the yield function.
Return Mapping Algorithm:
Two cases exist due to cone geometry:
- Return to smooth cone portion:
- Update deviatoric stress: \(\mathbf{s}_{n+1} = \left(1 - \frac{G\Delta\gamma}{\sqrt{J_2(s_{trial})}}\right) \mathbf{s}_{trial}\)
- Update mean pressure: \(p_{n+1} = p_{trial} - K\bar{\eta}\Delta\gamma\)
- Consistency condition gives scalar equation for \(\Delta\gamma\)
- Return to apex (when stress tries to enter the “forbidden” region):
- Deviatoric stress collapses: \(\mathbf{s}_{n+1} = 0\)
- Mean pressure determined from apex condition: \(p = \frac{\xi}{\eta}c\)
- Different scalar equation for volumetric plastic strain increment
The algorithm first attempts the smooth cone; if invalid, it applies the apex return.
11.13 Consistent Tangent Moduli
For Newton-Raphson convergence in FEM, we need the algorithmic (consistent) tangent stiffness \(\mathbb{C}_{\text{ep}}\).
For Drucker-Prager (Smooth Cone Return): \[\mathbb{C}_{\text{ep}} = 2G\left(1 - \frac{\Delta\gamma}{\sqrt{2}\|\varepsilon^e_{trial,\text{dev}}\|}\right)\mathbb{I}_{\text{dev}} + 2G\left(\frac{\Delta\gamma}{\sqrt{2}\|\varepsilon^e_{trial,\text{dev}}\|} - GA\right)\mathbf{D}\otimes\mathbf{D}\] \[\quad - \sqrt{2}GAK(\eta\mathbf{D}\otimes\mathbf{I} + \bar{\eta}\mathbf{I}\otimes\mathbf{D}) + K(1 - K\eta\bar{\eta}A)\mathbf{I}\otimes\mathbf{I}\]
where \(A = \frac{1}{G + K\eta\bar{\eta} + \xi^2 H_\text{iso}}\) and \(H_\text{iso}\) is the hardening slope.
For Drucker-Prager (Apex Return): \[\mathbb{C}_{\text{ep}} = K\left(1 - \frac{K}{K + \alpha^2 H_\text{iso}}\right)\mathbf{I}\otimes\mathbf{I}\]
where \(\alpha = \xi/\bar{\eta}\).
The algorithm uses the tangent consistent with the return mapping (smooth or apex) that was active in the previous iteration.
11.14 Mohr-Coulomb and Drucker-Prager: Key Relationships
Material parameters can be related via:
| Parameter | Expression |
|---|---|
| \(\eta\) (DP yield) | \(\sin\phi / \sqrt{9 + 3\sin^2\phi}\) |
| \(\xi\) (DP yield) | \(3c\cos\phi / \sqrt{9 + 3\sin^2\phi}\) |
| \(\bar{\eta}\) (DP flow) | Often chosen as \(\sin\psi\) (dilatancy angle) |
These relationships allow fitting Drucker-Prager parameters to match Mohr-Coulomb behavior while maintaining smoothness for numerical robustness.