K. FeikeP. KurzejaJ. MoslerK. LangenfeldTU Dortmund University, Institute of Mechanics, Leonhard-Euler-Str. 5, D-44227 Dortmund, Germany

###### Abstract

The stress triaxiality and the Lode angle parameter are two well established stress invariants for the characterization of damage evolution. This work assesses the limits of this tuple by using it for damage predictions in a continuum damage mechanics framework. Isotropic and anisotropic formulations of two well-established models are used to avoid model-specific restrictions. The damage evolution is analyzed for different load paths, while the stress triaxiality and the Lode angle parameter are controlled. The equivalent plastic strain is moreover added as a third parameter, but still does not suffice to uniquely define the damage state. As a consequence, well-established concepts such as fracture surfaces depending on this triple have to be taken with care, if complex paths are to be investigated. These include, e.g., load paths observed during metal forming applications with varying load directions or multiple stages.

###### keywords:

ductile damage, triaxiality, lode angle, path dependence

^{†}

^{†}journal: Engineering Fracture Mechanics

## 1 Introduction

The prediction of damage in materials is still a major challenge in many engineering disciplines since the underlying mechanisms are even nowadays not fully understood. As a result, large safety factors are mandatory to statistically ensure the functionality of structures and components. The modeling of the underlying material behavior is thus a key aspect for understanding and predicting the performance of these parts, cf. [1, 2, 3].

From a microscopic view, the mechanism of ductile damage can be related to the nucleation, the growth and finally the coalescence of voids induced by plastic deformation[4]. Depending on the material of interest, voids can show different origins and shapes, e.g., decohesion at hard particles. Since voids reduce the effective (load-carrying) area, they result in a reduced effective material stiffness[5], which ultimately also reduces the load-bearing capacity of the material. A natural modeling framework is thus the effective stress concept, cf.[6, 7]. Within this concept, which is also known as principle of strain equivalence, the true stresses are introduced. They are defined with respect to the effective cross sectional area due to void nucleation, growth and coalescence. An alternative modeling concept is based on effective strains, which can be interpreted as the decomposition of the total strains into elastic and damage-related parts[8, 9]. Finally, an effective damage-induced elastic stiffness can also be derived by means of the concept of strain energy equivalence[10, 11, 12]. In contrast to the concept of effective stresses and the concept of effective strains, the variational principle of strain energy equivalence naturally enforces physically relevant invariances. For instance, it naturally leads to a symmetric Cauchy stress even for anisotropic damage degradation (certainly, provided a Boltzmann continuum is considered).

Turning the view to the macroscopic description, a common approach for damage characterization is based on stress invariants like triaxiality $\eta$ and Lode angle parameter $\bar{\theta}$

$\displaystyle\eta$ | $\displaystyle=\dfrac{\sqrt{2}}{3}\dfrac{\sigma_{1}+\sigma_{2}+\sigma_{3}}{%\sqrt{\left[\sigma_{1}-\sigma_{2}\right]^{2}+\left[\sigma_{2}-\sigma_{3}\right%]^{2}+\left[\sigma_{3}-\sigma_{1}\right]^{2}}}$ | (1) | ||

$\displaystyle\bar{\theta}$ | $\displaystyle=1-\dfrac{2}{\pi}\,\arccos\!\left({\frac{L\,\left[L-3\right]\,%\left[L+3\right]}{\left[L^{2}+3\right]^{3/2}}}\right)\qquad\text{with}\quad L=%\dfrac{2\,\sigma_{2}-\sigma_{1}-\sigma_{3}}{\sigma_{1}-\sigma_{3}}$ | (2) |

combined with an equivalent plastic strain, cf.[13, 14], via

$\displaystyle\varepsilon^{\mathrm{p,eq}}$ | $\displaystyle=\int\limits_{0}^{t}\sqrt{\dfrac{2}{3}\,\dot{\boldsymbol{%\varepsilon}}^{\text{p}}:\dot{\boldsymbol{\varepsilon}}^{\text{p}}}\,\mathrm{d%}\bar{t}\,\text{.}$ | (3) |

$L$ is the Lode parameter, $\sigma_{i}$ are the stress eigenvalues and $\dot{\boldsymbol{\varepsilon}}^{\text{p}}$ is the plastic strain rate tensor. From a mathematical point of view, a characterization by the aforementioned invariants corresponds to an isotropic approximation. Early research showed that the stress triaxiality can be related to the physical process of void nucleation, growth and coalescence, cf.[15, 16]. Likewise, the Lode angle parameter was found to characterize damage evolution at low stress triaxialities, see[14, 17]. More recently, the Lode angle parameter has been investigated for different fracture criteria in order to characterize the onset of damage[18]. The entire triple of stress triaxiality, Lode angle parameter and equivalent plastic strain can be found in the context of so-called fracture surfaces[13, 14, 19]. A fracture surface usually represents an iso-surface in the space of $\{\eta,\,\bar{\theta},\,\varepsilon^{\mathrm{p,eq}}\}$, at which macroscopic failure is observed in experiments. This concept can also be adapted for damage models[20, 21, 22].

To investigate the effect of the stress state on ductile damage, loads at constant stress triaxiality have been analyzed in[23, 24, 25] with a focus on the microscale. Further studies have extended the analyses to include or to additionally control the Lode (angle) parameter[26, 27, 28, 29]. Findings clearly indicate a dependency of the underlying microstructure (distribution of voids in a unit cell) on the macroscopic response[29]. This indicates limits of common isotropic approaches for damage characterization.

Additional influences of the complete stress state on the damage evolution are often not considered, though, for instance those accounting for anisotropic behavior. Moreover, the stress invariants stress triaxiality and Lode angle parameter often change in time for complex load paths. Therefore, such load paths – as occuring in forming processes[30, 31] – require an extended description.

The present work thus analyzes the limits of simplified isotropic damage characterizations such as that by Eqs.(1)–(3) and identifies suitable and unsuitable load conditions. It uses numerical models for this purpose, which we understand as one major scientific instrument next to experimental investigations. While experimental data will clearly be needed for a final confirmation of the modeling predictions, the present numerical approach benefits from the possibility of broadband parameter studies. This numerical investigation shall thus systematically disclose the missing relationships and guide future experimental exploration.

The leading property of interest is ductile damage evolution in the context of continuum damage mechanics. The key highlights are:

- 1.
The frequently applied recommendation ”the lower the stress triaxiality, the lower the damage accumulation” is shown to be generally invalid for complex load paths.

- 2.
Such load paths are identified and discussed by means of computational optimization.

- 3.
An extended characterization of damage by the triple of stress triaxiality, Lode (angle) parameter and equivalent plastic strain, as employed in the concept of fracture surfaces, is shown to be still insufficient.

- 4.
Finally, even a characterization of damage by a complete isotropic representation of the stress tensor (all invariants) is generally insufficient for complex load paths.

In order to elaborate these highlights, two anisotropic damage models depending on the complete stress tensor are considered. Their prediction capabilities are then evaluated, for instance, by comparison of two different load paths with identical stress triaxialities and/or Lode angle parameters. Since such numerical experiments may also depend on the underlying constitutive model, two different modeling frameworks are adopted. While the first one is based on the principle of energy equivalence[10, 11, 12], the second one is the by now classic Lemaitre model[32].

This work is structured as follows. Section 2 describes the two constitutive models based on continuum damage mechanics, which form the basis for the subsequent investigations. The numerical methodology for the discretization, control and optimization of load paths is presented in Section 3. Particularly, techniques for prescribing the stress-triaxiality and the Lode (angle) parameter are presented. The main findings are reported in Section 4 by numerical experiments that highlight the limits of different isotropic damage characterizations. Section 5 concludes with implications for the requirement of a more general damage characterization.

## 2 Continuum damage mechanics – two prototype models

Two constitutive models are briefly presented for the prediction of ductile damage, both as isotropic and anisotropic formulations. This allows a comprehensive load path analysis that is not bound to the characteristics of a single model. The first model, ECC, is a characteristic prototype of the effective configuration concept (also known by the principle of strain energy equivalence). The second model, Lemaitre model LEM, is a characteristic prototype of the effective stress concept (also known by the principle of strain equivalence). All models have been calibrated with the same experimental data and form the basis for the subsequent numerical investigations.

### 2.1 ECC-Model – Effective Configuration Concept

The first model is named after the effective configuration concept and is based on the principle of strain energy equivalence, cf.[10]. The constitutive relations are adopted from the works in [11, 12].The Helmholtz energy $\psi$ is additively decomposed into an elastic part ($\psi^{\text{e}}$) and a plastic part ($\psi^{\text{p}}$) according to

$\displaystyle\psi$ | $\displaystyle=\psi^{\text{e}}(\boldsymbol{\varepsilon},\,\boldsymbol{%\varepsilon}^{\text{p}},\,\boldsymbol{b})+\psi^{\text{p}}(\boldsymbol{a},\,%\boldsymbol{k},\,\boldsymbol{b},\,\boldsymbol{\varepsilon}^{\text{p}})\,\text{,}$ | (4) | ||

$\displaystyle\psi^{\text{e}}$ | $\displaystyle=\dfrac{\lambda}{2}\,\left[\boldsymbol{b}:\boldsymbol{\varepsilon%}^{\text{e}}_{+}+\boldsymbol{I}:\boldsymbol{\varepsilon}^{\text{e}}_{-}\right]%^{2}+\mu\,\left[\boldsymbol{b}:\left[\boldsymbol{\varepsilon}^{\text{e}}_{+}%\cdot\boldsymbol{b}\cdot\boldsymbol{\varepsilon}^{\text{e}}_{+}\right]+%\boldsymbol{\varepsilon}^{\text{e}}_{-}:\boldsymbol{\varepsilon}^{\text{e}}_{-%}\right]\,\text{,}$ | (5) | ||

$\displaystyle\psi^{\text{p}}$ | $\displaystyle=\dfrac{H_{i}}{2}\,\left[\boldsymbol{b}:\boldsymbol{k}\right]^{2}%+\dfrac{H_{\mathrm{a}}}{2}\,\boldsymbol{b}:\left[\boldsymbol{a}\cdot%\boldsymbol{b}\cdot\boldsymbol{a}\right]\text{.}$ | (6) |

The strain tensor $\boldsymbol{\varepsilon}$ is also additively decomposed into an elastic and a plastic part, i.e., $\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{\text{e}}+\boldsymbol{%\varepsilon}^{\text{p}}$, cf.[33].The model captures isotropic as well as kinematic hardening due to the state variables $\boldsymbol{k}$ and $\boldsymbol{a}$, respectively.Here, isotropic hardening is captured by a tensorial internal variable, following the principle of strain energy equivalence [12]. Both the elastic and the plastic energy contribution are affected by the damage evolution due to integrity tensor $\boldsymbol{b}$. The integrity tensor can be interpreted as an inverse damage measure, i.e., $\boldsymbol{b}=\boldsymbol{I}$ represents the virgin material, while at least one of the eigenvalues of $\boldsymbol{b}$ converge to zero for completely damaged material points. By choosing $\boldsymbol{b}=\boldsymbol{I}$, Hooke’s law is recovered. Accordingly, $\lambda$ and $\mu$ are the Lame parameters. Furthermore, $H_{\mathrm{k}}$ is the isotropic hardening modulus and $H_{\mathrm{a}}$ is the kinematic hardening modulus.

A crucial contribution to the material model is the microcrack-closure-reopening (MCR) effect, which only allows damage evolution under tensile modes. Within the present model, the MCR effect is incorporated through an additive decomposition of the elastic strain tensor into its tensile and its compression modes according to

$\displaystyle\boldsymbol{\varepsilon}^{\text{e}}=\boldsymbol{\varepsilon}^{%\text{e}}_{+}+\boldsymbol{\varepsilon}^{\text{e}}_{-}\text{,}\quad\boldsymbol{%\varepsilon}^{\text{e}}_{+}$ | $\displaystyle=\sum_{i=1}^{3}\mathcal{H}\!\left({\varepsilon^{\text{e}}_{i}}%\right)\boldsymbol{N}_{i}\otimes\boldsymbol{N}_{i}.$ | (7) |

Here, $\varepsilon^{\text{e}}_{i}$ are the eigenvalues, $\boldsymbol{N}_{i}$ are the associated eigenvectors of the elastic strain tensor and $\mathcal{H}$ is the Heaviside function.Within the numerical implementation, the discontinuous Heaviside function has been approximated in line with[12] by setting numerical parameters $g_{0}=0$, $x_{0}=0$ and $x_{R}=10^{-6}$. A straightforward calculation of the thermodynamic driving forces leads to

$\displaystyle\boldsymbol{\sigma}\phantom{{}^{\text{e}}}$ | $\displaystyle=\phantom{-}\dfrac{\partial\psi\phantom{{}^{\text{e}}}}{\partial%\boldsymbol{\varepsilon}}=\lambda\,\left[\boldsymbol{b}:\boldsymbol{%\varepsilon}^{\text{e}}_{+}+\boldsymbol{I}:\boldsymbol{\varepsilon}^{\text{e}}%_{-}\right]\,\left[\boldsymbol{b}:\dfrac{\partial\boldsymbol{\varepsilon}^{%\text{e}}_{+}}{\partial\boldsymbol{\varepsilon}^{\text{e}}}+\boldsymbol{I}:%\dfrac{\partial\boldsymbol{\varepsilon}^{\text{e}}_{-}}{\partial\boldsymbol{%\varepsilon}^{\text{e}}}\right]+2\,\mu\,\left[\left[\boldsymbol{b}\cdot%\boldsymbol{\varepsilon}^{\text{e}}_{+}\cdot\boldsymbol{b}\right]:\dfrac{%\partial\boldsymbol{\varepsilon}^{\text{e}}_{+}}{\partial\boldsymbol{%\varepsilon}^{\text{e}}}+\boldsymbol{\varepsilon}^{\text{e}}_{-}:\dfrac{%\partial\boldsymbol{\varepsilon}^{\text{e}}_{-}}{\partial\boldsymbol{%\varepsilon}^{\text{e}}}\right]\,\text{,}$ | (8) | ||

$\displaystyle\boldsymbol{\alpha}\phantom{{}^{\text{e}}}$ | $\displaystyle=-\dfrac{\partial\psi\phantom{{}^{\text{e}}}}{\partial\boldsymbol%{a}}=-H_{a}\,\boldsymbol{b}\cdot\boldsymbol{a}\cdot\boldsymbol{b}\,\text{,}$ | (9) | ||

$\displaystyle\boldsymbol{\kappa}\phantom{{}^{\text{e}}}$ | $\displaystyle=-\dfrac{\partial\psi\phantom{{}^{\text{e}}}}{\partial\boldsymbol%{k}}=-H_{i}\,\left[\boldsymbol{b}:\boldsymbol{k}\right]\,\boldsymbol{b}\,\text%{,}$ | (10) | ||

$\displaystyle\boldsymbol{\beta}^{\text{e}}$ | $\displaystyle=-\dfrac{\partial\psi^{\text{e}}}{\partial\boldsymbol{b}}=-%\lambda\,\left[\boldsymbol{b}:\boldsymbol{\varepsilon}^{\text{e}}_{+}+%\boldsymbol{I}:\boldsymbol{\varepsilon}^{\text{e}}_{-}\right]\boldsymbol{%\varepsilon}^{\text{e}}_{+}-2\,\mu\,\boldsymbol{\varepsilon}^{\text{e}}_{+}%\cdot\boldsymbol{b}\cdot\boldsymbol{\varepsilon}^{\text{e}}_{+}\,\text{,}$ | (11) | ||

$\displaystyle\boldsymbol{\beta}^{\text{p}}$ | $\displaystyle=-\dfrac{\partial\psi^{\text{p}}}{\partial\boldsymbol{b}}=-H_{a}%\,\boldsymbol{a}\cdot\boldsymbol{b}\cdot\boldsymbol{a}-H_{i}\,\left[%\boldsymbol{b}:\boldsymbol{k}\right]\,\boldsymbol{k}\,\text{,}$ | (12) | ||

$\displaystyle\boldsymbol{\beta}\phantom{{}^{\text{e}}}$ | $\displaystyle=\boldsymbol{\beta}^{\text{e}}+\boldsymbol{\beta}^{\text{p}}\,%\text{,}$ | (13) |

where $\boldsymbol{\sigma}$ is the stress tensor, $\boldsymbol{\alpha}$ is the back stress tensor, $\boldsymbol{\kappa}$ is the drag stress tensor and $\boldsymbol{\beta}$ the energy release rate. Using the thermodynamic driving forces, the dissipation inequality reads in its reduced form

$\displaystyle\mathcal{D}^{\text{red}}$ | $\displaystyle=\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}-\dot{\psi}=%\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}^{\text{p}}+\boldsymbol{%\alpha}:\dot{\boldsymbol{a}}+\boldsymbol{\kappa}:\dot{\boldsymbol{k}}+%\boldsymbol{\beta}:\dot{\boldsymbol{b}}\geq 0.$ | (14) |

By following the principle of Generalized Standard Materials [34], a plastic potential is introduced for the definition of the evolution equations. The dissipation inequality is known to be automatically fulfilled if that potential is convex, non-negative and contains the origin. The plastic potential is specified as

$\displaystyle g$ | $\displaystyle=\Phi(\boldsymbol{\sigma},\,\boldsymbol{\alpha},\,\boldsymbol{%\kappa,\,\boldsymbol{b}})+\Gamma_{\alpha}(\boldsymbol{\alpha},\,\boldsymbol{b}%)+\Gamma_{\beta}(\boldsymbol{\beta}^{\text{e}},\,\boldsymbol{b})$ | (15) |

and consists of the yield function $\Phi$ and non-associated parts $\Gamma_{\!\alpha}$ and $\Gamma_{\beta}$.The yield function is chosen as

$\displaystyle\Phi$ | $\displaystyle=\sqrt{\bar{\tau}}-\tau_{y}-\frac{1}{3}\,\boldsymbol{b}^{-1}:%\boldsymbol{\kappa}-\Delta\tau\,\left[1-\mathrm{exp}\left(-\frac{\left|%\boldsymbol{b}^{-1}:\boldsymbol{\kappa}\right|}{\kappa_{\mathrm{u}}}\right)%\right]\,\text{,}$ | (16) | ||

$\displaystyle\bar{\tau}$ | $\displaystyle=\dfrac{3}{2}\,\boldsymbol{b}^{-1}:\left[\boldsymbol{\tau}\cdot%\boldsymbol{b}^{-1}\cdot\boldsymbol{\tau}\right]-\dfrac{1}{2}\,\left[%\boldsymbol{b}^{-1}:\boldsymbol{\tau}\right]^{2}\quad\text{with }\boldsymbol{%\tau}=\boldsymbol{\sigma}-\boldsymbol{\alpha}\,\text{,}$ | (17) |

where the equivalent stress $\sqrt{\bar{\tau}}$ is of von Mises-type and depends on the relative stress tensor $\boldsymbol{\tau}$.According to Eq.(16), the yield function captures linear and exponential isotropic hardening. The respective model parameters are denoted as $\tau_{y}$, $\Delta\tau$ and $\kappa_{\mathrm{u}}$.The non-associative parts of potential(15) are given as

$\displaystyle\Gamma_{\alpha}$ | $\displaystyle=\dfrac{B_{a}}{2\,H_{a}}\,\boldsymbol{b}^{-1}:\left[\boldsymbol{%\alpha}\cdot\boldsymbol{b}^{-1}\cdot\boldsymbol{\alpha}\right]\,\text{,}$ | (18) | ||

$\displaystyle\Gamma_{\beta}$ | $\displaystyle=\dfrac{C_{i}}{2}\,\left[\boldsymbol{b}^{m}:\boldsymbol{\beta}^{%\text{e}}\right]^{2}+\dfrac{C_{a}}{2}\,\boldsymbol{b}^{m}:\left[\boldsymbol{%\beta}^{\text{e}}\cdot\boldsymbol{b}^{m}\cdot\boldsymbol{\beta}^{\text{e}}\right]$ | (19) |

and extend the model to non-linear kinematic hardening of Armstrong-Frederick type by means of model parameter $B_{\mathrm{a}}$, cf.[35].The damage evolution is controlled by potential $\Gamma_{\beta}$ and consists of an isotropic part (damage modulus $C_{i}$) and an anisotropic part (damage modulus $C_{a}$).The damage exponent $m$ allows for greater flexibility when calibrating the model to experiments. Furthermore, potential $\Gamma_{\beta}$ depends only on the elastic part of the energy release rate $\boldsymbol{\beta}$. This novel modification of the original model ensures that no damage evolution occurs under compression. The evolution equations follow as gradients of plastic potential(15) – in line with the principle of Generalized Standard Materials – and read

$\displaystyle\dot{\boldsymbol{\varepsilon}}^{\text{p}}=\dot{\lambda}\,\dfrac{%\partial g}{\partial\boldsymbol{\sigma}},\quad\dot{\boldsymbol{a}}=\dot{%\lambda}\,\dfrac{\partial g}{\partial\boldsymbol{\alpha}},\quad\dot{%\boldsymbol{k}}=\dot{\lambda}\,\dfrac{\partial g}{\partial\boldsymbol{\kappa}}%,\quad\dot{\boldsymbol{b}}=\dot{\lambda}\,\dfrac{\partial g}{\partial%\boldsymbol{\beta}}=\dot{\lambda}\,\dfrac{\partial g}{\partial\boldsymbol{%\beta}^{\text{e}}}\text{,}$ | (20) |

which concludes the constitutive model.

##### Isotropic model

The isotropic version of the anisotropic ECC-Model follows straightforward by setting anisotropic damage modulus $C_{a}=0$. This model will also be employed within the numerical experiments.

### 2.2 LEM-Model – Effective Stress Concept (Lemaitre model)

The second model is based on the effective stress concept with equivalent strains and is adopted from[32]. Gibbs energy $\mathcal{G}$ is specified as

$\displaystyle\mathcal{G}$ | $\displaystyle=\mathcal{G}^{\text{e}}+\boldsymbol{\sigma}:\boldsymbol{%\varepsilon}^{\text{p}}-\psi^{\text{p}}\,\text{,}$ | (21) | ||

$\displaystyle\mathcal{G}^{\text{e}}$ | $\displaystyle=\frac{1+2\,\nu}{2\,E}\,\left[\boldsymbol{H}:\left[\boldsymbol{%\sigma}_{+}^{\mathrm{dev}}\cdot\boldsymbol{H}\cdot\boldsymbol{\sigma}_{+}^{%\mathrm{dev}}\right]+\boldsymbol{\sigma}_{-}^{\mathrm{dev}}:\boldsymbol{\sigma%}_{-}^{\mathrm{dev}}\right]+\frac{1-2\,\nu}{2\,E}\,\left[\frac{\left\langle%\sigma^{\mathrm{hyd}}\right\rangle^{2}}{1-D^{\mathrm{hyd}}}+\left\langle-%\sigma^{\mathrm{hyd}}\right\rangle^{2}\right]\,\text{,}$ | (22) | ||

$\displaystyle\psi^{\text{p}}$ | $\displaystyle=\dfrac{H_{a}}{2}\,\boldsymbol{a}:\boldsymbol{a}+\Delta\tau\,%\left[k+\kappa_{u}\,\mathrm{exp}\left(-\dfrac{k}{\kappa_{u}}\right)\right]\,$ | (23) |

and depends on model parameters Young’s modulus $E$, Poisson’s ratio $\nu$, the isotropic hardening parameters $\Delta\tau$ and $c$, and the kinematic hardening modulus $K$. The additive decomposition of the total strains into elastic and plastic parts is again considered as $\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{\text{e}}+\boldsymbol{%\varepsilon}^{\text{p}}$. Furthermore, the stress tensor is decomposed into a deviatoric ($\boldsymbol{\sigma}^{\mathrm{dev}}$) and a hydrostatic part ($\sigma^{\mathrm{hyd}}\,\boldsymbol{I}$) and subsequently into their respective positive and negative parts in order to capture the MCR-effect. Strain-like variables $\boldsymbol{a}$ and $p$ are suitable for kinematic and isotropic hardening, respectively. The damage evolution is captured by means of second-order tensor $\boldsymbol{D}$, which enters energy(22) in terms of variables

$\displaystyle D^{\mathrm{hyd}}$ | $\displaystyle\coloneqq\frac{1}{3}\,\mathrm{tr}\left(\boldsymbol{D}\right)\text%{,}\qquad\qquad\boldsymbol{H}\coloneqq\left[\boldsymbol{I}-\boldsymbol{D}%\right]^{-\frac{1}{2}}\,\text{.}$ |

While undamaged material is represented by $\boldsymbol{D}\!=\!\boldsymbol{0}$, at least one eigenvalue of $\boldsymbol{D}$ converges to one for completely damaged material points.

The thermodynamic driving forces follow as the derivatives of $\mathcal{G}$ and read

$\displaystyle\phantom{-}\boldsymbol{\varepsilon}^{\text{e}}$ | $\displaystyle=\phantom{-}\dfrac{\partial\mathcal{G}}{\partial\boldsymbol{%\sigma}}=\frac{1+\nu}{E}\,\left[\mathrm{dev}\left(\boldsymbol{H}\cdot%\boldsymbol{\sigma}_{+}^{\mathrm{dev}}\cdot\boldsymbol{H}\right)+\boldsymbol{%\sigma}_{-}^{\mathrm{dev}}\right]+\frac{1-2\,\nu}{E}\,\left[\frac{\left\langle%\sigma^{\mathrm{hyd}}\right\rangle}{1-D^{\mathrm{hyd}}}-\left\langle-\sigma^{%\mathrm{hyd}}\right\rangle\right]\,\boldsymbol{I}\,\text{,}$ | (24) | ||

$\displaystyle\phantom{-}\boldsymbol{\alpha}$ | $\displaystyle=-\dfrac{\partial\mathcal{G}}{\partial\boldsymbol{a}}=\frac{2}{3}%\,H_{a}\,\boldsymbol{a}\,\text{,}$ | (25) | ||

$\displaystyle\phantom{-}\kappa$ | $\displaystyle=-\dfrac{\partial\mathcal{G}}{\partial k}=\Delta\tau\,\left[1-%\mathrm{exp}\left(-\dfrac{k}{\kappa_{u}}\right)\right]\,\text{,}$ | (26) | ||

$\displaystyle\phantom{-}\boldsymbol{Y}$ | $\displaystyle=-\dfrac{\partial\mathcal{G}}{\partial\boldsymbol{D}}=-\frac{1+2%\,\nu}{E}\,\left[\boldsymbol{\sigma}^{\mathrm{dev}}_{+}\cdot\boldsymbol{H}%\cdot\boldsymbol{\sigma}^{\mathrm{dev}}_{+}\right]:\dfrac{\partial\boldsymbol{%H}}{\partial\boldsymbol{D}}-\frac{1-2\,\nu}{6\,E}\,\frac{\left\langle\sigma^{%\mathrm{hyd}}\right\rangle^{2}}{\left[1-D^{\mathrm{hyd}}\right]^{2}}\,%\boldsymbol{I}\,\text{.}$ | (27) |

The transformation of Gibbs energy $\mathcal{G}$ into Helmholtz energy $\psi$ reads

$\displaystyle\psi$ | $\displaystyle=\boldsymbol{\sigma}:\boldsymbol{\varepsilon}-\mathcal{G}=%\boldsymbol{\sigma}:\boldsymbol{\varepsilon}-\mathcal{G}^{\text{e}}-%\boldsymbol{\sigma}:\boldsymbol{\varepsilon}^{\text{p}}+\psi^{\text{p}}\,$ | (28) |

and allows to derive the reduced dissipation inequality as

$\displaystyle\mathcal{D}^{\mathrm{red}}$ | $\displaystyle=\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}-\dot{\psi}=%\boldsymbol{\sigma}:\dot{\boldsymbol{\varepsilon}}^{\text{p}}-\boldsymbol{Y}:%\dot{\boldsymbol{D}}-\boldsymbol{\alpha}:\dot{\boldsymbol{a}}-\kappa\,\dot{k}%\geq 0\,\text{.}$ | (29) |

Following the principle of Generalized Standard Materials once again, plastic potential $g$ is specified. Analogously to the previous model, $g$ is assumed to be of type

$\displaystyle g$ | $\displaystyle=\Phi\!\left({\boldsymbol{\sigma},\boldsymbol{\alpha},\kappa,%\boldsymbol{D}}\right)+\Gamma_{\!\alpha}\!\left({\boldsymbol{\alpha}}\right)$ | (30) |

with von Mises-type yield function

$\displaystyle\Phi$ | $\displaystyle=\sqrt{\bar{\tau}}-\tau_{\mathrm{y}}-\kappa\quad\text{with}\,\,%\bar{\tau}=\frac{3}{2}\,\boldsymbol{\tau}:\boldsymbol{\tau}\text{,}\quad%\boldsymbol{\tau}=\boldsymbol{H}\cdot\boldsymbol{\sigma}^{\mathrm{dev}}\cdot%\boldsymbol{H}-\boldsymbol{\alpha}$ | (31) |

and initial yield limit $\tau_{y}$. The non-associative part of potential(30) reads

$\displaystyle\Gamma_{\!\alpha}$ | $\displaystyle=\frac{3\,B_{a}}{4\,H_{a}}\,\boldsymbol{\alpha}:\boldsymbol{\alpha}$ | (32) |

and extends the model to non-linear kinematic hardening due to model parameter $B_{a}$. The evolution equations then result in

$\displaystyle\dot{\boldsymbol{\varepsilon}}^{\text{p}}$ | $\displaystyle=\dot{\lambda}\,\dfrac{\partial g}{\partial\boldsymbol{\sigma}}\,%\text{,}\qquad\dot{\boldsymbol{a}}=-\dot{\lambda}\,\dfrac{\partial g}{\partial%\boldsymbol{\alpha}}\,\text{,}\qquad\dot{k}=-\dot{\lambda}\,\dfrac{\partial g}%{\partial\kappa}=\dot{\lambda}\,\text{.}$ | (33) |

Following[32], the evolution equation associated with damage tensor $\boldsymbol{D}$ does not follow from potential $g$, but is explicitly specified as

$\displaystyle\dot{\boldsymbol{D}}$ | $\displaystyle=\left[C\,\bar{Y}\right]^{m}\,\left[\dot{\boldsymbol{\varepsilon}%}^{\text{p}}_{+}-\dot{\boldsymbol{\varepsilon}}^{\text{p}}_{-}\right]\,\text{,}$ | (34) | ||

$\displaystyle\text{with }\bar{Y}$ | $\displaystyle=\frac{1+\nu}{2\,E}\,\mathrm{tr}\left(\left[\boldsymbol{H}\cdot%\boldsymbol{\sigma}^{\mathrm{dev}}_{+}\cdot\boldsymbol{H}\right]^{2}\right)+%\frac{3-6\,\nu}{2\,E}\,\frac{\left\langle{\sigma^{\mathrm{hyd}}}\right\rangle^%{2}}{\left[1-D^{\mathrm{hyd}}\right]^{2}}$ | (35) |

and depends on model parameters $C$ and $m$. This evolution is formulated in terms of a scalar-valued energy release rate $\bar{Y}$ and the additive decomposition of the plastic strain rate tensor $\dot{\boldsymbol{\varepsilon}}^{\text{p}}$ into its positive and negative parts, cf.Eq.(7).

##### Isotropic model

The isotropic version of the anisotropic ECC-model results through the modification of evolution Eq.(34) into a spherical counterpart. Here, the choice

$\displaystyle\dot{\boldsymbol{D}}$ | $\displaystyle=\left[\hat{C}\,\bar{Y}\right]^{\hat{m}}\,\dot{\varepsilon}^{%\mathrm{p,eq}}\,\boldsymbol{I}$ | (36) |

is made since it leads to similar results as the underlying anisotropic model when considering a uniaxial tensile test.

### 2.3 Model calibration

All four models — the isotropic and the anisotropic ECC model (ECC-i and ECC-a) as well as the isotropic and the anisotropic LEM model (LEM-i and LEM-a) — are calibrated to experimental data from[36] for comparability. This involves tensile experiments of case-hardened steel 16MnCrS5 up to strain amplitudes of $\varepsilon_{11}\!\approx 0.13$, which represents a typical material used in bulk forming such as forward extrusion. The strains remain spatially homogeneous in this deformation range and necking does not yet occur. Both models, the anisotropic ECC model (ECC-a) and the anisotropic LEM model (LEM-a) were calibrated by means of a standard least-squares approach. All models share the same elastic parameters $\lambda$ and $\mu$ and the same initial yield stress $\tau_{\mathrm{y}}$. Isotropic and anisotropic variants of the same model also share the same material parameters for plasticity. A visualization of the calibrations is shown in Fig.1 and the final list of model parameters is given in Tabs.1 and 2.

Parameter name | Symbol | Value | Unit |

First Lame parameter | $\lambda$ | $118870$ | $\mathrm{MPa}$ |

Second Lame parameter | $\mu$ | $79249$ | $\mathrm{MPa}$ |

Yield stress | $\tau_{\mathrm{y}}$ | $308.260$ | $\mathrm{MPa}$ |

Kinematic hardening modulus | $H_{a}$ | $7728.863$ | $\mathrm{MPa}$ |

Kinematic yield stress saturation | $B_{a}$ | $38.218$ | - |

Isotropic hardening modulus | $H_{i}$ | $1.829\cdot 10^{-4}$ | $\mathrm{MPa}$ |

Isotropic hardening increment | $\Delta\tau$ | $2.261\cdot 10^{-2}$ | $\mathrm{MPa}$ |

Isotropic hardening parameter | $\kappa_{\mathrm{u}}$ | $0.159$ | $\mathrm{MPa}$ |

Anisotropic (ECC-a model) | |||

Isotropic damage modulus | $C_{i}$ | $0.0$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Anisotropic damage modulus | $C_{a}$ | $14.503$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Damage exponent | $m$ | $11.217$ | - |

Isotropic (ECC-i model) | |||

Isotropic damage modulus | $\hat{C}_{i}$ | $14.408$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Anisotropic damage modulus | $\hat{C}_{a}$ | $0.0$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Damage exponent | $\hat{m}$ | $11.373$ | - |

Parameter name | Symbol | Value | Unit |

First Lame parameter | $\lambda$ | $118875.0$ | $\mathrm{MPa}$ |

Second Lame parameter | $\mu$ | $79250.0$ | $\mathrm{MPa}$ |

Initial yield stress | $\tau_{\mathrm{y}}$ | $308.26$ | $\mathrm{MPa}$ |

Kinematic hardening modulus | $H_{a}$ | $3774.25$ | $\mathrm{MPa}$ |

Kinematic yield stress saturation | $B_{a}$ | $175.55$ | - |

Saturation hardening stress | $\Delta\tau$ | $1.1761\cdot 10^{6}$ | $\mathrm{MPa}$ |

Isotropic hardening parameter | $\kappa_{u}$ | $301.41$ | - |

Anisotropic (LEM-a model) | |||

Damage modulus | $C$ | $1256.7$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Damage exponent | $m$ | $0.2$ | - |

Isotropic (LEM-i model) | |||

Damage modulus | $\hat{C}$ | $623.1$ | $\mathrm{M}\mathrm{Pa}^{-1}$ |

Damage exponent | $\hat{m}$ | $0.2$ | - |

All models match the experimental data well, see Fig.1 (a, b). The largest error between the experiments and the calibration is associated with the LEM-a model and occurs at a strain of $\varepsilon_{11}=0.07$. However, the respective error is below $10\text{\,}\mathrm{MPa}$ and hence, less than $2.5\%$. The behavior with respect to a shear mode is additionally shown in Fig.1 (c) in order to highlight the models’ different response – the shear mode was not part of the calibration. All models show qualitatively similar behavior with a maximum deviation of less than $8\%$. It shall be underlined that different models are used on purpose for a better and model-independent interpretation of the results. An evaluation of the individual models themselves is not intended and clearly depends on the intended application.

###### Remark 1

The constitutive models are evaluated at the material point level and calibrated to experiments with spatially constant stress and strain fields (no necking). Therefore, no regularization is required. If boundary value problems with localization are to be analyzed, a regularization is required in order to render finite element simulations well-posed. The respective energy contributions are, however, disregarded in this work, cf.[37].

## 3 Methodology of evaluation: load path parametrization and damage measures

The analysis of varying load paths is a key of the present evaluation of damage descriptions. For example, the predictions of the previously introduced models will be compared along load paths of identical stress triaxialities and/or Lode angle parameters. For that reason, the parametrization of the load paths is explained in the following Subsection3.1. Characteristic damage measures are furthermore defined in Subsection3.2 and based on elastic stiffness and elastic compliance tensors. Each of these measures allows to evaluate specific types of material degradation. For example, the Frobenius norm of the degraded elasticity tensor can be interpreted as an average stiffness measure. These characteristic measures will take the role of target functions for the load path analysis. They allow the comparison of load paths with respect to their damage impact, the identification of sensitivities and the exploration of missing influences.

### 3.1 Parametrization of load paths

The numerical studies of the next section analyze the behavior of a representative material point. The mechanical loading of these points will hence be defined by prescribing either components of the stress or the strain tensor. Since these tensors are energetically dual, symmetric second-order tensors, one has to prescribe precisely six components in total. The dual six components follow as reactions. For instance, a strain-driven uniaxial tensile test corresponds to the following strain and stress tensors — boxed components are prescribed, the other ones follow as reactions:

$\displaystyle\boldsymbol{\varepsilon}=\begin{bmatrix}\boxed{\varepsilon_{11}}&%&\varepsilon_{12}&&\varepsilon_{13}\\\text{sym}&&\varepsilon_{22}&&\varepsilon_{23}\\\text{sym}&&\text{sym}&&\varepsilon_{33}\end{bmatrix}_{\left[\boldsymbol{e}_{1%},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\right]}\qquad\boldsymbol{\sigma}=%\begin{bmatrix}\sigma_{11}&&\boxed{\sigma_{12}}&&\boxed{\sigma_{13}}\\\text{sym}&&\boxed{\sigma_{22}}&&\boxed{\sigma_{23}}\\\text{sym}&&\text{sym}&&\boxed{\sigma_{33}}\end{bmatrix}_{\left[\boldsymbol{e}%_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}\right]}$ | (37) |

In Eq.(37), the components are defined with respect to a chosen basis $\boldsymbol{e}_{i}$ (here, a cartesian setting is adopted). Loading is implemented by enforcing a strain in $11$-direction for this example of a uniaxial tensile test. Hence, the respective components of the stress tensor are unknown, but follow from the constitutive model. The other components of the stress tensor have to vanish, prescribing $\sigma_{ij}=0$ except for $\sigma_{11}$.

The evolution of prescribed components is either parametrized and discretized by piecewise linear functions of the form

$\displaystyle f_{\mathrm{pwl}}\!\left({t;\boldsymbol{P}}\right)$ | $\displaystyle=P_{i}+\frac{P_{i}-P_{j}}{t_{i}-t_{j}}\,t\quad\text{if }t\in\left%[t_{i},\,t_{j}\right]$ | (38) |

or by Bézier curves defined as

$\displaystyle f_{\mathrm{Bez}}\!\left({t;\boldsymbol{Q}}\right)$ | $\displaystyle=\sum_{i=0}^{N}{N\choose i}\,\left[1-t\right]^{\left[N-i\right]}%\,t^{i}\,Q_{i}$ | (39) |

with control points $P_{i}$ and $Q_{i}$, respectively. The time is normalized within the interval $t\in[0;1]$, where the time-scaling is irrelevant since rate-independent models are considered.

From the perspective of the numerical implementation, prescribing components of the stress tensor requires an iteration. This is implemented by a constitutive driver where the respective constraints are solved by means of a Newton-iteration. An analogous procedure is also employed for controlling the stress triaxiality and/or the Lode (angle) parameter.

### 3.2 Characteristic damage measures and target functions

Characteristic damage measurements are defined in order to compare the damage evolution along the load paths and between different models. Since the internal variables capturing damage are not identical for both basic models, the respective elastic stiffness and elastic compliance tensors are instead considered for the definition of suitable damage measures. Elastic stiffness tensor and elastic compliance tensors are inverse to each other, i.e., $\mathbb{E}^{\text{e}}:\mathbb{C}^{\text{e}}=\mathbb{I}$ with $\mathbb{I}$ being the symmetric fourth-order identity. Starting from the Helmholtz energy $\psi$, the elastic moduli of the ECC-model are obtained as

$\displaystyle\left[\mathbb{E}_{\mathrm{\text{ECC}}}^{\text{e}}\right]_{ijkl}:=%\frac{\partial^{2}\psi}{\partial\varepsilon_{ij}\partial\varepsilon{kl}}=%\lambda\,\left[b_{ij}\,b_{kl}\right]+\mu\left[b_{ik}\,b_{jl}+b_{il}\,b_{jk}\right]$ | (40) |

with $\mathbb{E}_{\mathrm{\text{ECC}}}^{\text{e}}:\mathbb{C}_{\mathrm{\text{ECC}}}^{%\text{e}}=\mathbb{I}$. Likewise, the LEM-model yields the elastic compliance

$\displaystyle\begin{split}\left[\mathbb{C}_{\mathrm{\text{LEM}}}^{\text{e}}%\right]_{ijkl}:=\frac{\partial^{2}\mathcal{G}}{\partial\sigma_{ij}\partial%\sigma_{kl}}=&\frac{1+\nu}{E}\,\left[\frac{1}{2}\,\left[H_{ik}\,H_{jl}+H_{il}%\,H_{jk}\right]-\frac{1}{3}\,\left[\delta_{ij}\,H_{km}\,H_{ml}+H_{in}\,H_{nj}%\,\delta_{kl}\right]\right.\\&\phantom{\coloneqq\frac{1+\nu}{E}\,}+\left.\dfrac{1}{9}\,\left[\,H_{op}\,H_{%op}+\frac{3}{1-D^{\mathrm{hyd}}}\right]\,\delta_{ij}\,\delta_{kl}\right]-\frac%{\nu}{E}\,\frac{1}{1-D^{\mathrm{hyd}}}\,\delta_{ij}\,\delta_{kl}\end{split}$ | (41) |

where $\delta_{ij}$ is the Kronecker-Delta and $\mathbb{E}_{\mathrm{\text{LEM}}}^{\text{e}}:\mathbb{C}_{\mathrm{\text{LEM}}}^{%\text{e}}=\mathbb{I}$.

Effective-scalar-valued variables can be introduced based on either the elastic stiffness or the elastic compliance tensor. For instance, the Frobenius-norm of the fourth-order elasticity tensor

$\displaystyle f_{\mathbb{E}}\coloneqq\sqrt{\mathbb{E}^{\text{e}}::\mathbb{E}^{%\text{e}}}$ | (42) |

with $::$ denoting a quadruple contraction, defines a suitable average value of $\mathbb{E}^{\text{e}}$. Based on Eq.(42) and denoting the elasticity tensor associated with a reference model or load path as $\bar{\mathbb{E}}$, the difference between two load paths (or model descriptions) can be defined as

$\displaystyle f_{\mathbb{E}}^{\mathrm{diff}}$ | $\displaystyle=\left[{\sqrt{\mathbb{E}^{\text{e}}::\mathbb{E}^{\text{e}}}-\sqrt%{\bar{\mathbb{E}}^{\text{e}}::\bar{\mathbb{E}}^{\text{e}}}}\right]^{2}.$ | (43) |

The undamaged configuration represents a natural reference solution, i.e., $\bar{\mathbb{E}}^{\text{e}}=\mathbb{E}^{\text{e}}(t=0)$.As an alternative to the Frobenius-norm(42), the directional elastic stiffness and compliance

$\displaystyle E_{\boldsymbol{r}}=\left[\boldsymbol{r}\otimes\boldsymbol{r}%\right]:\mathbb{E}^{\text{e}}:\left[\boldsymbol{r}\otimes\boldsymbol{r}\right]$ | (44) | ||

$\displaystyle C_{\boldsymbol{r}}=\left[\boldsymbol{r}\otimes\boldsymbol{r}%\right]:\mathbb{C}^{\text{e}}:\left[\boldsymbol{r}\otimes\boldsymbol{r}\right]$ | (45) |

can be considered respectively, where unit vector $\boldsymbol{r}$ defines the direction. A suitable measure for the damage accumulation based on Eq.(44) reads, for instance,

$\displaystyle\xi_{\mathbb{E}}$ | $\displaystyle=\frac{\min_{\boldsymbol{r}}E_{\boldsymbol{r}}\!\left({%\boldsymbol{r}}\right)}{E_{0}}.$ | (46) |

Here, $E_{0}$ corresponds to the initial state at $t=0$, which is isotropic and independent of $\boldsymbol{r}$. In contrast to measure(43), measure(46) depends on the directional-dependent lowest stiffness and thus largest damage-induced degradation. Finally, if the direction-dependent elastic modulus is to be evaluated, $C_{\boldsymbol{r}}$ can be compared to the reference resulting in

$\displaystyle\xi_{\mathbb{C}}(\boldsymbol{r})=\dfrac{C_{0}}{C_{\boldsymbol{r}}}.$ | (47) |

## 4 Numerical examples

This section highlights the application limits of isotropic damage representations by means of illustrative examples. These examples are structured by the questioning of well-established statements on damage characterization and control. This allows a connection to practically relevant load path domains and to better identify influences that have not yet been considered in full depth. First, the well-known statement

- 1.
The smaller the stress triaxiality, the smaller the damage accumulation

(see[38] and Subsection4.1)

as well as

- 1.
Damage accumulation is uniquely governed by the stress triaxiality and the Lode (angle) parameter (see[21, 39])

are analyzed. It will be shown that contradictory load paths can indeed be designed by means of optimization algorithms. Second, the concept of fracture surfaces in line with[13, 14] is investigated in Subsection4.2. This concept is based on the statment

- 1.
Damage accumulation is uniquely governed by the triple stress triaxiality, Lode (angle) parameter and equivalent plastic strain (see[20, 21, 22] and Subsection4.2).

Optimization algorithms again show that this statement does also not hold true in general – particularly for complex load paths. Since the present studies rely on numerical constitutive modeling, the two different protoype models ECC and LEM with isotropic and anisotropic variants are considered throughout the discussion.

### 4.1 The smaller the stress triaxiality, the smaller the damage accumulation?

Coupled tension-torsion experiments

Inspired by experimental axial-torsional test systems and with the goal of multi-directional load paths in mind, a tension-torsion problem is considered within a cylindrical coordinate system $\left[\boldsymbol{e}_{\mathrm{r}},\boldsymbol{e}_{\theta},\boldsymbol{e}_{z}\right]$, cf.Fig.2.

$\displaystyle\boldsymbol{\varepsilon}$ | $\displaystyle=\begin{bmatrix}\varepsilon_{rr}&&\varepsilon_{r\theta}&&%\varepsilon_{rz}\\\text{sym}&&\varepsilon_{\theta\theta}&&\boxed{\varepsilon_{\theta z}}\\\text{sym}&&\text{sym}&&\boxed{\varepsilon_{zz}}\end{bmatrix}_{\left[%\boldsymbol{e}_{\mathrm{r}},\boldsymbol{e}_{\theta},\boldsymbol{e}_{z}\right]}$ | ||

$\displaystyle\boldsymbol{\sigma}$ | $\displaystyle=\begin{bmatrix}\boxed{\sigma_{rr}}&&\boxed{\sigma_{r\theta}}&&%\boxed{\sigma_{rz}}\\\text{sym}&&\boxed{\sigma_{\theta\theta}}&&\sigma_{\theta z}\\\text{sym}&&\text{sym}&&\sigma_{zz}\end{bmatrix}_{\left[\boldsymbol{e}_{%\mathrm{r}},\boldsymbol{e}_{\theta},\boldsymbol{e}_{z}\right]}$ |

A virtual tension-torsion sequence serves as the reference load path. Within this test, torsional strain $\varepsilon_{\theta z}$ is increased first up to $\varepsilon_{\theta z}=0.03$ (while $\varepsilon_{zz}=0=\text{const.}$), and subsequently tensile strain $\varepsilon_{zz}$ is increased up to $\varepsilon_{zz}=0.06$ (while $\varepsilon_{\theta z}=0.03=\text{const.}$). The final strain components at normalized time $t=1$ are $\varepsilon_{zz}=0.06$ and $\varepsilon_{\theta z}=0.03$. The evolution of the stress-triaxiality corresponding to this reference load path is denoted by $\eta^{\text{ref}}(t)$ and the Frobenius-norm of the final elastic stiffness tensor is $f_{\mathbb{E}}^{\text{ref}}(t=1)$. Since four different models are analyzed, four different reference solutions are computed (isotropic and anisotropic ECC-model, isotropic and anisotropic LEM-model).

Next, load path variations are parametrized to search for damage-mitigating alternatives with higher stress triaxialities. These would then provide counter examples for the first statement (”The smaller the stress triaxiality, the smaller the damage accumulation”). Strain components $\varepsilon_{\theta z}$ and $\varepsilon_{zz}$ are discretized by means of fifth-order Bézier curves, see Eq.(39). While the initial conditions are $\varepsilon_{\theta z}(t=0)=\varepsilon_{zz}(t=0)=0$, the final state is defined by $\varepsilon_{\theta z}(t=1)=0.03$ and $\varepsilon_{zz}(t=1)=0.06$. The other coordinates of the Bézier representation are computed by mitigating the final damage state in the sense of an optimization. To be more precise, the following optimization problem is considered:

$\displaystyle\underset{Q_{i}}{\text{max}}\,f_{\mathbb{E}}\!\left({t=1}\right)%\quad s.t.\,\,\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)\,\forall\,t\in%\left[0,1\right]$ | (48) |

Accordingly, we search for the load path in terms of control points $Q_{i}$ that maximizes the average elastic stiffness. This path has to fulfill the constraint $\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)\,\forall\,t$. Since the resulting optimization problem is highly non-convex, a two-stage gradient-free optimization method is employed. To be more explicit, a particle swarm optimization (see[40]), followed up by a Nelder-Mead simplex algorithm (see[41]) is employed. The inequality $\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)\,\forall\,t$ is enforced by means of a penalty function.

The results in Fig.3 correspond to eight different simulations – four constitutive models (isotropic and anisotropic ECC-model, isotropic and anisotropic LEM-model) times two load paths (reference path and optimized path).

80pt][t]0.07\psfrag{A}[l][l]{\scalebox{0.7}{\text{ECC-a} ref.}}\psfrag{B}[l][l]{\scalebox{0.7}{\text{ECC-i} ref.}}\psfrag{C}[l][l]{\scalebox{0.7}{\text{LEM-a} ref.}}\psfrag{D}[l][l]{\scalebox{0.7}{\text{LEM-i} ref.}}\psfrag{E}[l][l]{\scalebox{0.7}{\text{ECC-a} opt.}}\psfrag{F}[l][l]{\scalebox{0.7}{\text{ECC-i} opt.}}\psfrag{G}[l][l]{\scalebox{0.7}{\text{LEM-a} opt.}}\psfrag{H}[l][l]{\scalebox{0.7}{\text{LEM-i} opt.}}\includegraphics[width=368.57964pt]{Chapter/0X_NumericalResults/TensionShear/Figures/Legend2.eps}

It can be seen that the optimized paths indeed lead to less damage, i.e., the average elastic stiffness is larger for such paths, see the right columns in Fig.3 (a). Since ductile damage is primarily driven by plastic deformations, one might expect that $\varepsilon^{\text{p,eq}}$ is smaller for the optimized paths. This is confirmed in Fig.3 (b). For this reason, a low equivalent plastic strain is suggested as a further damage-mitigating indicator. This hypothesis can also be drawn from the monotonically decreasing courses in Fig.3 (c). To be more specific, Figs.3 (a, b) are end-point projections of the curves presented in Fig.3 (c) onto the $\xi_{\mathbb{E}}$-axis and the $\varepsilon^{\text{p,eq}}$-axis, respectively. Although the curves associated with optimal paths are below their respective reference (i.e. less remaining stiffness at distinct values of $\varepsilon^{\text{p,eq}}$), the final remaining stiffness can be increased by reducing the path length in $\varepsilon^{\text{p,eq}}$. Therefore, the equivalent plastic strain $\varepsilon^{\text{p,eq}}$ is an additional important measure to characterize ductile damage, alongside of the stress triaxiality and the Lode angle parameter. Its influence will be examined in the next subsection.

According to Fig.4 (b), the stress triaxiality of the optimized paths is always higher than that of the reference paths, i.e., inequality $\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)$ is indeed fulfilled at any time.

80pt][t]0.07\psfrag{A}[l][l]{\scalebox{0.7}{\text{ECC-a} ref.}}\psfrag{B}[l][l]{\scalebox{0.7}{\text{ECC-i} ref.}}\psfrag{C}[l][l]{\scalebox{0.7}{\text{LEM-a} ref.}}\psfrag{D}[l][l]{\scalebox{0.7}{\text{LEM-i} ref.}}\psfrag{E}[l][l]{\scalebox{0.7}{\text{ECC-a} opt.}}\psfrag{F}[l][l]{\scalebox{0.7}{\text{ECC-i} opt.}}\psfrag{G}[l][l]{\scalebox{0.7}{\text{LEM-a} opt.}}\psfrag{H}[l][l]{\scalebox{0.7}{\text{LEM-i} opt.}}\includegraphics[width=390.25534pt]{Chapter/0X_NumericalResults/TensionShear/Figures/Legend2.eps}

As evident from Fig.4 (a), shear component $\varepsilon_{\theta z}$ and axial component $\varepsilon_{zz}$ of the strain tensor are simultaneously applied for the optimal paths – in contrast to the sequential reference path. Finally, the evolution of the Lode angle parameter is shown in Fig.4 (c). The evolution is similar to that of the stress triaxiality. In particular, the Lode angle parameter of the optimal path is always larger than its referential counterpart. The present experiment includes no means to control the Lode angle parameter. However, literature indicates that it is an important indicator for ductile damage, along with stress triaxiality. Therefore, the Lode angle parameter will be controlled in subsequent experiments.

A more detailed comparison of the material degradation at the final time $t=1$ is given in Fig.5. It shows the directional relative stiffness of the anisotropic damage models ECC-a and LEM-a relative to the initial isotropic, undamaged state. While the model LEM-a leads to an almost isotropic material degradation, the model ECC-a predicts a more pronounced anisotropy. It shall be noted again that the comparison of the models does not aim at their individual evaluation but rather at finding model-invariant conclusions on the damage evolution and the corresponding input parameters. For both constitutive models, the optimized load paths are characterized by a larger elastic stiffness compared to the reference solution. Stress triaxiality – even in combination with the equivalent plastic strain – hence remains an insufficient indicator for damage.

\psfrag{e1}[c][c]{$\boldsymbol{e}_{\mathrm{z}}$}\psfrag{e2}[c][c]{$\boldsymbol{e}_{\theta}$}\psfrag{e3}[c][c]{$\boldsymbol{e}_{\mathrm{r}}$}\includegraphics[width=260.17464pt]{Chapter/0X_NumericalResults/TensionShear/Figures/CoordinateSystemPSFRAG.eps}

\psfrag{E1}[c][c]{\makebox{\raisebox{4.0pt}{$\xi_{\mathbb{C}}$}}}\includegraphics[width=303.53267pt]{Chapter/0X_NumericalResults/TensionShear/Figures/StiffnessEkh_Leg.eps}

### 4.2 Is damage accumulation uniquely governed by the triple stress triaxiality, Lode (angle) parameter and equivalent plastic strain?

It has been shown in the previous section that the stress triaxiality as well as the equivalent plastic strain are primary influencing factors of ductile damage evolution. However, it has also been shown that these factors are not sufficient in order to characterize anisotropic material degradation under complex load paths. For this reason, the Lode angle parameter is often also considered. So-called fracture surfaces are a frequently employed damage framework that is based on the stress triaxiality, the Lode angle parameter and the equivalent plastic strain, cf.[13, 14]. More precisely, the damage evolution is then assumed to be proportional to the plastic strain rate and a scaling factor. The latter is a function in terms of the triple stress triaxiality, Lode angle parameter and the equivalent plastic strain.Fracture surfaces, such as those in[13, 14], can be interpreted as iso-surfaces of an equivalent scalar-valued damage measure. Accordingly, they can be computed for a given model by connecting triples of stress triaxiality, Lode angle parameter and the equivalent plastic strain to the same equivalent damage measure. First, such surfaces are computed in paragraph4.2.1 for proportional load paths. Subsequently, the influence of more complex load paths is studied in paragraph4.2.2 as well as in paragraph4.2.3. It will be shown in this section that even this ansatz is still not sufficient for a full characterization of ductile damage accumulation.

#### 4.2.1 Fracture surfaces for proportional loading

In order to compute iso-surfaces of damage, a suitable scalar-valued measure is required. In this section, the directional relative stiffness of $\xi_{\mathbb{E}}=0.8$ is chosen for that purpose. Clearly other measures and thresholds are also possible. Starting from the initial unloaded configuration, the stresses are increased until threshold $\xi_{\mathbb{E}}=0.8$ is reached. For the whole load path, the stress triaxiality and the Lode angle parameter are kept constant (proportional loading). $17\times 13$ tuples of stress triaxiality and Lode angle parameter are considered in order to span the fracture surface. Along the evolving load path, the stress tensor is prescribed as

$\displaystyle\boldsymbol{\sigma}=\begin{bmatrix}\boxed{\sigma_{11}}&&\boxed{0}%&&\boxed{0}\\\text{sym}&&\boxed{\sigma_{22}(\sigma_{11},\eta,\bar{\theta})}&&\boxed{0}\\\text{sym}&&\text{sym}&&\boxed{\sigma_{33}(\sigma_{11},\eta,\bar{\theta})}\end%{bmatrix}_{\left[\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}%\right]}$ | (49) |

where $\eta$ is the (prescribed) stress triaxiality and $\bar{\theta}$ the (prescribed) Lode angle parameter. Closed-form solutions for $\sigma_{22}(\sigma_{11},\eta,L\!\left({\bar{\theta}}\right))$ and $\sigma_{33}(\sigma_{11},\eta,L\!\left({\bar{\theta}}\right))$ are given inA. Thus, the only free loading parameter is $\sigma_{11}$. This parameter is increased until the damage threshold is reached.

In line with [13, 14], the ECC model shows the trend the smaller the stress triaxialities, the larger the equivalent plastic strains – both for the anisotropic and for the isotropic formulation, see Fig.7. The ECC model does not predict damage evolution for biaxial compression. Although the stress triaxiality is the major influencing factor of the stress state, the fracture surface corresponding to the ECC model also depends on the Lode angle parameter with the largest sensitivity for small stress triaxialities.

\psfrag{E}[c][c]{\scalebox{0.8}{Equivalent plastic strain $\varepsilon^{\text{p,eq}}$ [-]}}\includegraphics[width=34.69038pt]{Chapter/0X_NumericalResults/FractureSurfaces/Figures/LEG.eps}

In contrast to the ECC model, the fracture surface predicted by the LEM model is less sensitive with respect to the stress triaxialities and the Lode angle parameter, see Figs.7 (c) and (d). However, the Lode angle parameter again contributes separately to the damage evolution and the overall trend is similar: The smaller the stress triaxialities, the larger the equivalent plastic strains. Only for larger negative triaxialities a plateau can be seen – in contrast to model ECC. This plateau can be partly explained by the structure of the driving force, cf. Eq.(35). Since driving force $\bar{Y}$ is quadratic in the stress tensor, the pre-factor scaling the damage evolution is symmetric in the stress triaxiality. Clearly, the MCR-Effect dampens this effect.

#### 4.2.2 Fracture surfaces for non-proportional load paths

This paragraph analyzes whether the fracture surface is also meaningful for other loads paths (in contrast to the proportional ones before with constant stress triaxiality and Lode angle parameter). As a representative reference path, uniaxial tension is considered, i.e., $\boldsymbol{\sigma}=\sigma_{11}\,\boldsymbol{e}_{1}\otimes\boldsymbol{e}_{1}$. It corresponds to a stress triaxiality of $\eta=1/3$ and Lode angle parameter of $\bar{\theta}=1$. Inspired by uniaxial tension, the following more general stress path is considered:

$\displaystyle\boldsymbol{\sigma}(t)=\begin{bmatrix}\boxed{\sigma_{11}(t)}&&%\boxed{0}&&\boxed{0}\\\text{sym}&&\boxed{\sigma_{22}(\sigma_{11}(t),\eta(t))}&&\boxed{0}\\\text{sym}&&\text{sym}&&\boxed{0}\end{bmatrix}_{\left[\boldsymbol{e}_{1},%\boldsymbol{e}_{2},\boldsymbol{e}_{3}\right]}$ | (50) |

Clearly, by setting $\eta=1/3$ and $\bar{\theta}=1$, uniaxial tension is recovered. However, the stress triaxiality is prescribed by means of a Bézier interpolation in time. While it corresponds to uniaxial tension for the initial and the final state, i.e., $\eta(t=0)=\eta(t=1)=1/3$, it reaches its maximum deviation from this state at $t=0.5$ with a stress triaxiality of $\eta(t=0.5)=1/3-7/30=1/10$ for the first load path and $\eta(t=0.5)=1/3+7/30=17/30$ for the second load path. The Lode angle parameter is not explicitly controlled, but follows from $\sigma_{11}$ and $\eta$. Finally, it is noted that the time-scaling of $\sigma_{11}(t)$ is iteratively adapted until the state $\sigma_{11}(t=1)$ and $\eta(t=1)$ correspond to uniaxial tension and a damage threshold of $\xi_{\mathbb{E}}=0.8$.

The aforementioned two different load paths are shown in Fig.8 in dashed and dotted lines. They are represented by means of the triple equivalent plastic strain, stress triaxiality and Lode angle parameter. According to Fig.8, these paths are non-proportional and even more important, they lead to different equivalent plastic strains at the final time $t=1$ corresponding to $\xi_{\mathbb{E}}=0.8$, see the zoom-in in Fig.8 (right). The deviation of the equivalent plastic strain for these paths is above $30\%$. This deviation clearly confirms that a characterization and description of damage accumulation only by means of the three influencing variables — equivalent plastic strain, stress triaxiality and Lode angle parameter — is generally not sufficient. This holds in particular if complex, non-proportional load paths are to be analyzed. Since the path-dependent deviation of the fracture energy also relies on the underlying constitutive model, the numerical experiment has been recomputed and confirmed for all models (ECC-a, ECC-i, LEM-a, LEM-i); see the differences between the load paths in Fig.9.

Another perspective with relevant applications in, e.g., metal forming, is the prospect of damage control. In this setting, the aim is to minimize the accumulating damage during forming. Therefore and in contrast to the previous experiment the equivalent plastic strain is prescribed (together with triaxiality and Lode angle parameter), while the final damage state is free to evolve. The time-scaling of $\sigma_{11}(t)$ is iteratively adapted to synchronize $\eta(t=1)$ with the predefined equivalent plastic strain $\bar{\varepsilon}^{\mathrm{p,eq}}=0.0462$. This specific value corresponds to the respective location on the damage iso-surface (uniaxial tension, proportional load path, ECC-a, $\eta=1/3$, $\bar{\theta}=1$).

The results in Fig.10 show a dependence of the damage behavior on the path in the triple-space of stress triaxiality, Lode angle parameter and equivalent plastic strain. Choosing the virgin material ($\xi_{\mathbb{E}}=1$) as a reference point, lower stress triaxialities increase the remaining stiffness $\xi_{\mathbb{E}}$ by more than $9\,\%$, while the higher stress triaxialities decrease it by $4\,\%$.This viewpoint from damage control supports the previous findings of insufficient parameterization. Not only can different levels of equivalent plastic strain be achieved for a given damage threshold, but the inverse relationship holds as well. That is, different levels of damage can be achieved for a given equivalent plastic strain.

#### 4.2.3 Non-proportional load paths with constant stress invariants

According to the previous section, a damage characteriziation by means of the triple equivalent plastic strain, stress triaxiality and Lode angle parameter is usually not sufficient for complex load paths. Since the stress triaxiality and the Lode angle parameter are already two invariants of the stress tensor, one could even more so monitor and control all three invariants of the stress tensor (instead of just $\eta$ and $\bar{\theta}$). This is precisely the idea analyzed in this paragraph. For that purpose, two different stress tensors are considered. While $\boldsymbol{\sigma}$ corresponds to proportional loading, $\tilde{\boldsymbol{\sigma}}$ is associated with a more complex load path. Both share the same eigenvalues $\sigma_{i}\!\left({t}\right)=\tilde{\sigma}_{i}\!\left({t}\right)$ (and thus the same triaxiality and Lode parameter). But they have different directions in terms of the associated eigenvectors. $\boldsymbol{\sigma}$ and $\tilde{\boldsymbol{\sigma}}$ can be implemented by means of the spectral decomposition, i.e.,

$\displaystyle\boldsymbol{\sigma}$ | $\displaystyle=\sum_{i=1}^{3}\sigma_{i}\,\boldsymbol{N}^{i}\otimes\boldsymbol{N%}^{i}$ | (51) | ||

$\displaystyle\tilde{\boldsymbol{\sigma}}$ | $\displaystyle=\sum_{i=1}^{3}\sigma_{i}\,\left[\boldsymbol{R}\cdot\boldsymbol{N%}^{i}\right]\otimes\left[\boldsymbol{R}\cdot\boldsymbol{N}^{i}\right]\quad%\text{with }\boldsymbol{R}\!\left({\alpha,\beta,\gamma}\right)=\boldsymbol{R}_%{\bar{\bar{\boldsymbol{e}}}_{1}}\!\left({\alpha}\right)\cdot\boldsymbol{R}_{%\bar{\boldsymbol{e}}_{2}}\!\left({\beta}\right)\cdot\boldsymbol{R}_{%\boldsymbol{e}_{3}}\!\left({\gamma}\right)\,\text{.}$ | (52) |

where $\boldsymbol{N}_{i}$ are the eigenvectors of $\boldsymbol{\sigma}$ and $\boldsymbol{R}_{\boldsymbol{e}_{i}}$ are rotation matrices depending on the Euler angles $\alpha$, $\beta$ and $\gamma$. Four different mechanical tests are considered:

- 1.
Uniaxial tension: $\eta=1/3$ and $\bar{\theta}=1$

- 2.
Combined tension and shear: $\eta=1/6$ and $\bar{\theta}=0.506$

- 3.
Simple shear: $\eta=0$ and $\bar{\theta}=0$

- 4.
Combined compression and shear: $\eta=-1/6$ and $\bar{\theta}=-0.506$

While the first eigenvalue $\sigma_{1}$ is increased, the second as well as the third follow from the prescribed and constant parameters $\eta$ and $\bar{\theta}$, seeA for details. For the proportional load paths, stress tensor(51) is considered with constant eigenvectors $\boldsymbol{N}_{i}$ and $\sigma_{1}$ is increased up to a damage threshold of $\xi_{\mathbb{E}}=0.8$. In the case of the non-proportional load paths, the Euler angles $\alpha$, $\beta$ and $\gamma$ evolve as illustrated in Fig.11.In line with the example studied in the previous paragraph, the temporal discretization of $\sigma_{1}(t)$ is chosen such that all Euler angles vanish at final time $t=1$ when damage threshold $\xi_{\mathbb{E}}=0.8$ is reached, i.e., the time scaling of $\sigma_{1}$ and the Euler angles is synchronized.

The bar charts in Fig.12 summarize the results for all four models (ECC-a, ECC-i, LEM-a, LEM-i) and for all four mechanical tests (uniaxial tension, simple shear, combined tension and shear, combined compression and shear). They show a relative difference of the equivalent plastic strain at final damaged state $\xi_{\mathbb{E}}=0.8$ between the proportional and the non-proportional load paths. This implies again, that the chosen set of three independent stress invariants (or eigenvalues) is insufficient to fully characterize the damage evolution history. It can be seen that the anisotropic model ECC-a leads to the largest difference – followed by its isotropic counterpart ECC-i. The sensitivity on the load path is less pronounced for the LEM models. However, a minor influence is still evident – even for model LEM-i. While for some load paths, such invariant-based isotropic approximations might lead to reasonable predictions, they are generally inadequate for complex load paths.

\psfrag{e1}[c][c]{$\boldsymbol{e}_{x}$}\psfrag{e2}[c][c]{$\boldsymbol{e}_{y}$}\psfrag{e3}[c][c]{$\boldsymbol{e}_{z}$}\includegraphics[width=108.405pt]{Chapter/0X_NumericalResults/TensionShear/Figures/CoordinateSystemPSFRAG.eps}

## 5 Conclusion

Application limits of isotropic damage characterizations for complex load paths were investigated in this paper. Since such numerical experiments may depend on the underlying constitutive model, two well-established anisotropic damage models have been calibrated and considered individually for that purpose: The ECC-model is based on the principle of strain energy equivalence and the LEM-model is the Lemaitre-model based on effective stresses. Their prediction of ductile damage has been analyzed for different load paths.

The first numerical experiment – inspired by a coupled tension-torsion experiment – revealed that the frequently employed recommendation ”the smaller the stress triaxiality, the smaller the damage accumulation” can be wrong in general for complex load paths. In order to design paths contradicting this hypothesis, an algorithm was proposed to analyze differences in the predictions and unconsidered influences. Also the additional consideration of the equivalent plastic strain remained insufficient.

In addition to the invariant stress triaxiality, the Lode angle parameter was controlled in the next example. Since the characterization of damage by means of the stress triaxiality and the Lode angle parameter is closely related to the concept of fracture surfaces, the latter have been studied in detail. Again, load paths were designed showing that a damage characterization only by means of the extended triple of influencing variables stress triaxiality, Lode angle parameter and equivalent plastic strain is often still not possible. This limitation holds notably for both damage models and even their isotropic simplifications. Deviations due to unconsidered influences on fracture (iso-)surfaces have been detected by up to 30 % in terms of the equivalent plastic strain and up to 13 % of damage, respectively.

Finally, the full set of invariants of the stress tensor were controlled. This can be interpreted as an isotropic approximation. However, even in this more flexible case, load paths with different damage accumulation have been found, despite completely identical stress invariants along entire load paths.

The numerical studies showed that damage can indeed be influenced and improved by controlling the stress triaxiality as well as the Lode angle parameter. Also the equivalent plastic strain can be chosen as a possible control variable. However, the remaining coordinates of the load paths still have great potential for further improvement of damage accumulation and thus, control. This potential is to be further studied, e.g., in the context of forming processes. Furthermore, the numerical studies are to be further confirmed by real experiments. The load paths identified in this systematic analysis are intended to guide experimental load settings for the further exploration of damage-influencing factors.

## Appendix A Parametrization in terms of stress triaxiality and Lode parameter

The inversion of the equations for stress triaxiality(1)and Lode parameter(2), i.e.,

$\displaystyle\left[{\sigma_{2}}\right]_{1,2}$ | $\displaystyle=-\chi\,\sigma_{1}\pm\sqrt{\chi^{2}-\omega}\left|\sigma_{1}\right|$ | (53) | ||

$\displaystyle\left[{\sigma_{3}}\right]_{1,2}$ | $\displaystyle=\left[{a+b\,\chi}\right]\,\sigma_{1}\mp b\,\sqrt{\chi^{2}-\omega%}\,\left|\sigma_{1}\right|$ | (54) |

with

$\displaystyle a$ | $\displaystyle=\frac{L+1}{L-1},\quad b=\frac{2}{L-1}$ |

and

$\displaystyle\chi$ | $\displaystyle=\frac{18\,\eta^{2}\,\left[L-L^{2}-2\right]-4\,L\,\left[L-3\right%]}{18\,\eta^{2}\,\left[L^{2}-3\right]-2\,\left[L-3\right]^{2}},\quad\omega=%\frac{18\,\eta^{2}\,\left[L^{2}-L+2\right]-8\,L^{2}}{18\,\eta^{2}\,\left[L^{2}%-3\right]-2\,\left[L-3\right]^{2}}\text{.}$ |

allows to prescribe, and hence control, both the stress triaxiality as well as the Lode parameter. This allows to analyze load paths with identical stress triaxiality and Lode parameter with respect to damage evolution. Clearly, special care is required, if $L\rightarrow 1$.

## Acknowledgements

Financial support from the German Research Foundation (DFG) via SFB/TR TRR 188 (project number 278868966), project C01, is gratefully acknowledged. We also gratefully acknowledge the computing time provided on the Linux HPC cluster at TU Dortmund University (LiDO3), partially funded in the course of the Large-Scale Equiment Initiative by the German Research Foundation (DFG) (project number 271512359).

## References

- [1]R.I. Stephens, H.O. Fuchs (Eds.), Metal fatigue in engineering, 2nd Edition,A Wiley-Interscience publication, Wiley, New York Weinheim, 2001.
- [2]S.Ma, H.Yuan, Computational investigation of multi-axial damage modeling forporous sintered metals with experimental verification, Engineering FractureMechanics 149 (2015) 89–110.doi:10.1016/j.engfracmech.2015.09.049.
- [3]A.Tekkaya, P.-O. Bouchard, S.Bruschi, C.Tasan, Damage in metal forming, CIRPAnnals 69(2) (2020) 600–623.doi:10.1016/j.cirp.2020.05.005.
- [4]J.Koplik, A.Needleman, Void growth and coalescence in porous plastic solids,International Journal of Solids and Structures 24(8) (1988) 835–853.doi:10.1016/0020-7683(88)90051-0.
- [5]A.Sancho, M.Cox, T.Cartwright, G.Aldrich-Smith, P.Hooper, C.Davies,J.Dear, Experimental techniques for ductile damage characterisation,Procedia Structural Integrity 2 (2016) 966–973.doi:10.1016/j.prostr.2016.06.124.
- [6]J.Lemaitre, A Continuous Damage Mechanics Model for Ductile Fracture,Journal of Engineering Materials and Technology 107(1) (1985) 83–89.doi:10.1115/1.3225775.
- [7]A.L. Gurson, Continuum theory of ductile rupture by void nucleation andgrowth: Part i—yield criteria and flow rules for porous ductile media,Journal of Engineering Materials and Technology-transactions of The Asme 99(1977) 2–15.doi:10.1115/1.3443401.
- [8]J.Simo, J.-W. Ju, Strain- and stress-based continuum damage models. i.formulation, International Journal of Solids and Structures 23 (1987)821–840.doi:10.1016/0020-7683(87)90083-7.
- [9]M.Brünig, An anisotropic ductile damage model based on irreversiblethermodynamics, International Journal of Plasticity 19(10) (2003)1679–1713.doi:10.1016/S0749-6419(02)00114-6.
- [10]P.Steinmann, I.Carol, A framework for geometrically nonlinear continuumdamage mechanics, International Journal of Engineering Science 36(15) (1998)1793–1814.doi:10.1016/S0020-7225(97)00116-X.
- [11]A.Menzel, M.Ekh, P.Steinmann, K.Runesson, Anisotropic damage coupled toplasticity: Modelling based on the effective configuration concept,International Journal for Numerical Methods in Engineering 54(10) (2002)1409–1430.doi:10.1002/nme.470.
- [12]M.Ekh, A.Menzel, K.Runesson, P.Steinmann, Anisotropic damage with the mcreffect coupled to plasticity, International Journal of Engineering Science41(13) (2003) 1535–1551, damage and failure analysis of materials.doi:10.1016/S0020-7225(03)00032-6.
- [13]T.Wierzbicki, Y.Bao, Y.Bai, A New Experimental Technique forConstructing a Fracture Envelope of Metals under Multi-axialloading, Proceedings of the 2005 SEM Annual Conference and Exposition onExperimental and Applied Mechanics (2005).
- [14]Y.Bai, T.Wierzbicki, A new model of metal plasticity and fracture withpressure and lode dependence, International Journal of Plasticity 24(6)(2008) 1071–1096.doi:10.1016/j.ijplas.2007.09.004.
- [15]F.A. McClintock, A Criterion for Ductile Fracture by the Growth of Holes,Journal of Applied Mechanics 35(2) (1968) 363–371.doi:10.1115/1.3601204.
- [16]J.Rice, D.Tracey, On the ductile enlargement of voids in triaxial stressfields, Journal of the Mechanics and Physics of Solids 17(3) (1969)201–217.doi:10.1016/0022-5096(69)90033-7.
- [17]M.Brünig, S.Gerke, V.Hagenbrock, Micro-mechanical studies on the effect ofthe stress triaxiality and the Lode parameter on ductile damage,International Journal of Plasticity 50 (2013) 49–65.doi:10.1016/j.ijplas.2013.03.012.
- [18]A.Mattiello, R.Desmorat, Lode angle dependency due to anisotropic damage,International Journal of Damage Mechanics 30(2) (2021) 214–259.doi:10.1177/1056789520948563.
- [19]D.Anderson, C.Butcher, N.Pathak, M.Worswick, Failure parameteridentification and validation for a dual-phase 780 steel sheet, InternationalJournal of Solids and Structures 124 (2017) 89–107.doi:10.1016/j.ijsolstr.2017.06.018.
- [20]M.Basaran, Stress state dependent damage modeling with a focus on the lodeangle influence, Berichte aus dem Maschinenbau, Shaker, Aachen, 2011.
- [21]F.X.C. Andrade, M.Feucht, A.Haufe, F.Neukamm, An incremental stress statedependent damage model for ductile failure prediction, International Journalof Fracture 200(1-2) (2016) 127–150.doi:10.1007/s10704-016-0081-2.
- [22]F.Neukamm, Lokalisierung und Versagen von Blechstrukturen, no. Nr. 68 inBericht / Institut für Baustatik und Baudynamik der UniversitätStuttgart, Institut für Baustatik und Baudynamik der UniversitätStuttgart, Stuttgart, 2018.
- [23]M.Kuna, D.Z. Sun, Three-dimensional cell model analyses of void growth inductile materials, International Journal of Fracture 81(3) (1996) 235–258.doi:10.1007/BF00039573.
- [24]T.Pardoen, J.Hutchinson, An extended model for void growth and coalescence,Journal of the Mechanics and Physics of Solids 48(12) (2000) 2467–2512.doi:10.1016/S0022-5096(00)00019-3.
- [25]R.C. Lin, D.Steglich, W.Brocks, J.Betten, Performing RVE calculationsunder constant stress triaxiality for monotonous and cyclic loading,International Journal for Numerical Methods in Engineering 66(8) (2006)1331–1360, _eprint:https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1600.doi:10.1002/nme.1600.
- [26]K.S. Zhang, J.B. Bai, D.Fran, Numerical analysis of the in¯uence of theLode parameter on void growth, International Journal of Solids andStructures (2001).
- [27]X.Gao, J.Kim, Modeling of ductile fracture: Significance of voidcoalescence, International Journal of Solids and Structures 43(20) (2006)6277–6293.doi:10.1016/j.ijsolstr.2005.08.008.
- [28]R.Kiran, K.Khandelwal, A triaxiality and Lode parameter dependent ductilefracture criterion, Engineering Fracture Mechanics 128 (2014) 121–138.doi:10.1016/j.engfracmech.2014.07.010.
- [29]C.Cadet, J.Besson, S.Flouriot, S.Forest, P.Kerfriden, V.deRancourt,Ductile fracture of materials with randomly distributed voids, InternationalJournal of Fracture 230(1) (2021) 193–223.doi:10.1007/s10704-021-00562-7.
- [30]O.Hering, A.E. Tekkaya, Damage-induced performance variations of cold forgedparts, Journal of Materials Processing Technology 279 (2020) 116556.doi:10.1016/j.jmatprotec.2019.116556.
- [31]R.Gitschel, O.Hering, A.Schulze, A.ErmanTekkaya, Controlling DamageEvolution in Geometrically Identical Cold Forged Parts byCounterpressure, Journal of Manufacturing Science and Engineering 145(1)(2023) 011011.doi:10.1115/1.4056266.
- [32]J.Lemaitre, R.Desmorat, Engineering Damage Mechanics. Ductile, Creep, Fatigueand Brittle Failure, 2005.doi:10.1007/b138882.
- [33]A.E. Green, P.M. Naghdi, A general theory of an elastic-plasticcontinuum, Archive for Rational Mechanics and Analysis 18(4) (1965)251–281.doi:10.1007/BF00251666.
- [34]B.Halphen, Q.S. Nguyen, Sur les matériaux standardgénéralisés, 1975.
- [35]C.Frederick, P.Armstrong, A Mathematical Representation of theMultiaxial Bauscinger Effect, Materials at High Temperatures 24 (2007)1–26.doi:10.3184/096034007X207589.
- [36]K.Langenfeld, A.Schowtjak, R.Schulte, O.Hering, K.Möhring, T.Clausmeyer,R.Ostwald, F.Walther, A.Tekkaya, J.Mosler, Influence of anisotropicdamage evolution on cold forging, Production Engineering 14 (01 2020).doi:10.1007/s11740-019-00942-y.
- [37]K.Langenfeld, J.Mosler, A micromorphic approach for gradient-enhancedanisotropic ductile damage, Computer Methods in Applied Mechanics andEngineering 360 (2020) 112717.doi:10.1016/j.cma.2019.112717.
- [38]A.E. Tekkaya, N.BenKhalifa, O.Hering, R.Meya, S.Myslicki, F.Walther,Forming-induced damage and its effects on product properties, CIRP Annals66(1) (2017) 281–284.doi:10.1016/j.cirp.2017.04.113.
- [39]T.Zhu, H.Ding, C.Wang, Y.Liu, S.Xiao, G.Yang, B.Yang, ParametersCalibration of the GISSMO Failure Model for SUS301L-MT, ChineseJournal of Mechanical Engineering 36(1) (2023) 20.doi:10.1186/s10033-023-00844-2.
- [40]J.Kennedy, R.Eberhart, Particle swarm optimization, in: Proceedings ofICNN’95 - International Conference on Neural Networks, Vol.4, 1995, pp.1942–1948 vol.4.doi:10.1109/ICNN.1995.488968.
- [41]J.A. Nelder, R.Mead, A simplex method for function minimization, Comput. J. 7(1965) 308–313.doi:10.1093/comjnl/7.4.308.