Limits of isotropic damage models for complex load paths — beyond stress triaxiality and Lode angle parameter (2024)

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 investgated. 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: Journal of LaTeXTemplates

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 η𝜂\etaitalic_η and Lode angle parameter θ¯¯𝜃\bar{\theta}over¯ start_ARG italic_θ end_ARG

η𝜂\displaystyle\etaitalic_η=23σ1+σ2+σ3[σ1σ2]2+[σ2σ3]2+[σ3σ1]2absent23subscript𝜎1subscript𝜎2subscript𝜎3superscriptdelimited-[]subscript𝜎1subscript𝜎22superscriptdelimited-[]subscript𝜎2subscript𝜎32superscriptdelimited-[]subscript𝜎3subscript𝜎12\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}}}= divide start_ARG square-root start_ARG 2 end_ARG end_ARG start_ARG 3 end_ARG divide start_ARG italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT + italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT + italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_ARG start_ARG square-root start_ARG [ italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + [ italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + [ italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG(1)
θ¯¯𝜃\displaystyle\bar{\theta}over¯ start_ARG italic_θ end_ARG=12πarccos(L[L3][L+3][L2+3]3/2)withL=2σ2σ1σ3σ1σ3formulae-sequenceabsent12𝜋𝐿delimited-[]𝐿3delimited-[]𝐿3superscriptdelimited-[]superscript𝐿2332with𝐿2subscript𝜎2subscript𝜎1subscript𝜎3subscript𝜎1subscript𝜎3\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}}= 1 - divide start_ARG 2 end_ARG start_ARG italic_π end_ARG roman_arccos ( divide start_ARG italic_L [ italic_L - 3 ] [ italic_L + 3 ] end_ARG start_ARG [ italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + 3 ] start_POSTSUPERSCRIPT 3 / 2 end_POSTSUPERSCRIPT end_ARG ) with italic_L = divide start_ARG 2 italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_ARG start_ARG italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT - italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_ARG(2)

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

εp,eqsuperscript𝜀peq\displaystyle\varepsilon^{\mathrm{p,eq}}italic_ε start_POSTSUPERSCRIPT roman_p , roman_eq end_POSTSUPERSCRIPT=0t23𝜺˙p:𝜺˙pdt¯.absentsuperscriptsubscript0𝑡:23superscript˙𝜺psuperscript˙𝜺pdifferential-d¯𝑡.\displaystyle=\int\limits_{0}^{t}\sqrt{\dfrac{2}{3}\dot{\boldsymbol{%\varepsilon}}^{\text{p}}:\dot{\boldsymbol{\varepsilon}}^{\text{p}}}\,\mathrm{d%}\bar{t}\,\text{.}= ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_t end_POSTSUPERSCRIPT square-root start_ARG divide start_ARG 2 end_ARG start_ARG 3 end_ARG over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT : over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT end_ARG roman_d over¯ start_ARG italic_t end_ARG .(3)

L𝐿Litalic_L is the Lode parameter, σisubscript𝜎𝑖\sigma_{i}italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are the stress eigenvalues and 𝜺˙psuperscript˙𝜺p\dot{\boldsymbol{\varepsilon}}^{\text{p}}over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT 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 {η,θ¯,εp,eq}𝜂¯𝜃superscript𝜀peq\{\eta,\,\bar{\theta},\,\varepsilon^{\mathrm{p,eq}}\}{ italic_η , over¯ start_ARG italic_θ end_ARG , italic_ε start_POSTSUPERSCRIPT roman_p , roman_eq end_POSTSUPERSCRIPT }, 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 ψ𝜓\psiitalic_ψ is additively decomposed into an elastic part (ψesuperscript𝜓e\psi^{\text{e}}italic_ψ start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT) and a plastic part (ψpsuperscript𝜓p\psi^{\text{p}}italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT) according to

ψ𝜓\displaystyle\psiitalic_ψ=ψe(𝜺,𝜺p,𝒃)+ψp(𝒂,𝒌,𝒃,𝜺p),absentsuperscript𝜓e𝜺superscript𝜺p𝒃superscript𝜓p𝒂𝒌𝒃superscript𝜺p,\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{,}= italic_ψ start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ( bold_italic_ε , bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT , bold_italic_b ) + italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT ( bold_italic_a , bold_italic_k , bold_italic_b , bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT ) ,(4)
ψesuperscript𝜓e\displaystyle\psi^{\text{e}}italic_ψ start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT=λ2[𝒃:𝜺+e+𝑰:𝜺e]2+μ[𝒃:[𝜺+e𝒃𝜺+e]+𝜺e:𝜺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{,}= divide start_ARG italic_λ end_ARG start_ARG 2 end_ARG [ bold_italic_b : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + bold_italic_I : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + italic_μ [ bold_italic_b : [ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ⋅ bold_italic_b ⋅ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ] + bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] ,(5)
ψpsuperscript𝜓p\displaystyle\psi^{\text{p}}italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT=Hk2[𝒃:𝒌]2+Ha2𝒃:[𝒂𝒃𝒂].\displaystyle=\dfrac{H_{\mathrm{k}}}{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{.}= divide start_ARG italic_H start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG [ bold_italic_b : bold_italic_k ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG italic_H start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG bold_italic_b : [ bold_italic_a ⋅ bold_italic_b ⋅ bold_italic_a ] .(6)

The strain tensor 𝜺𝜺\boldsymbol{\varepsilon}bold_italic_ε is also additively decomposed into an elastic and a plastic part, i.e., 𝜺=𝜺e+𝜺p𝜺superscript𝜺esuperscript𝜺p\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{\text{e}}+\boldsymbol{%\varepsilon}^{\text{p}}bold_italic_ε = bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT + bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT, cf.[33].The model captures isotropic as well as kinematic hardening due to the state variables 𝒌𝒌\boldsymbol{k}bold_italic_k and 𝒂𝒂\boldsymbol{a}bold_italic_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}bold_italic_b. The integrity tensor can be interpreted as an inverse damage measure, i.e., 𝒃=𝑰𝒃𝑰\boldsymbol{b}=\boldsymbol{I}bold_italic_b = bold_italic_I represents the virgin material, while at least one of the eigenvalues of 𝒃𝒃\boldsymbol{b}bold_italic_b converge to zero for completely damaged material points. By choosing 𝒃=𝑰𝒃𝑰\boldsymbol{b}=\boldsymbol{I}bold_italic_b = bold_italic_I, Hooke’s law is recovered. Accordingly, λ𝜆\lambdaitalic_λ and μ𝜇\muitalic_μ are the Lame parameters. Furthermore, Hksubscript𝐻kH_{\mathrm{k}}italic_H start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT is the isotropic hardening modulus and Hasubscript𝐻aH_{\mathrm{a}}italic_H start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT 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

𝜺e=𝜺+e+𝜺e,𝜺+esuperscript𝜺esubscriptsuperscript𝜺esubscriptsuperscript𝜺e,subscriptsuperscript𝜺e\displaystyle\boldsymbol{\varepsilon}^{\text{e}}=\boldsymbol{\varepsilon}^{%\text{e}}_{+}+\boldsymbol{\varepsilon}^{\text{e}}_{-}\text{,}\quad\boldsymbol{%\varepsilon}^{\text{e}}_{+}bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT = bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT , bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT=i=13(εie)𝒎i.absentsuperscriptsubscript𝑖13subscriptsuperscript𝜀e𝑖subscript𝒎𝑖\displaystyle=\sum_{i=1}^{3}\mathcal{H}\!\left({\varepsilon^{\text{e}}_{i}}%\right)\boldsymbol{m}_{i}.= ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT caligraphic_H ( italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) bold_italic_m start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT .(7)

Here, εiesubscriptsuperscript𝜀e𝑖\varepsilon^{\text{e}}_{i}italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are the eigenvalues, 𝒎isubscript𝒎𝑖\boldsymbol{m}_{i}bold_italic_m start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are the eigenprojections of the elastic strain tensor and \mathcal{H}caligraphic_H is the Heaviside function.Within the numerical implementation, the discontinuous Heaviside function has been approximated in line with[12] by setting numerical parameters g0=0subscript𝑔00g_{0}=0italic_g start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, x0=0subscript𝑥00x_{0}=0italic_x start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0 and xR=106subscript𝑥𝑅superscript106x_{R}=10^{-6}italic_x start_POSTSUBSCRIPT italic_R end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT. A straightforward calculation of the thermodynamic driving forces leads to

𝝈𝝈\displaystyle\boldsymbol{\sigma}\phantom{{}^{\text{e}}}bold_italic_σ=ψ𝜺=λ[𝒃:𝜺+e+𝑰:𝜺e][𝒃:𝜺+e𝜺e+𝑰:𝜺e𝜺e]+2μ[[𝒃𝜺+e𝒃]:𝜺+e𝜺e+𝜺e:𝜺e𝜺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{,}= divide start_ARG ∂ italic_ψ end_ARG start_ARG ∂ bold_italic_ε end_ARG = italic_λ [ bold_italic_b : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + bold_italic_I : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] [ bold_italic_b : divide start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG + bold_italic_I : divide start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG ] + 2 italic_μ [ [ bold_italic_b ⋅ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ⋅ bold_italic_b ] : divide start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT end_ARG start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG + bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT : divide start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT end_ARG start_ARG ∂ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG ] ,(8)
𝜶𝜶\displaystyle\boldsymbol{\alpha}\phantom{{}^{\text{e}}}bold_italic_α=ψ𝒂=Ha𝒃𝒂𝒃,absent𝜓𝒂subscript𝐻a𝒃𝒂𝒃,\displaystyle=-\dfrac{\partial\psi\phantom{{}^{\text{e}}}}{\partial\boldsymbol%{a}}=-H_{\mathrm{a}}\,\boldsymbol{b}\cdot\boldsymbol{a}\cdot\boldsymbol{b}\,%\text{,}= - divide start_ARG ∂ italic_ψ end_ARG start_ARG ∂ bold_italic_a end_ARG = - italic_H start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT bold_italic_b ⋅ bold_italic_a ⋅ bold_italic_b ,(9)
𝜿𝜿\displaystyle\boldsymbol{\kappa}\phantom{{}^{\text{e}}}bold_italic_κ=ψ𝒌=Hk[𝒃:𝒌]𝒃,\displaystyle=-\dfrac{\partial\psi\phantom{{}^{\text{e}}}}{\partial\boldsymbol%{k}}=-H_{\mathrm{k}}\,\left[\boldsymbol{b}:\boldsymbol{k}\right]\,\boldsymbol{%b}\,\text{,}= - divide start_ARG ∂ italic_ψ end_ARG start_ARG ∂ bold_italic_k end_ARG = - italic_H start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT [ bold_italic_b : bold_italic_k ] bold_italic_b ,(10)
𝜷esuperscript𝜷e\displaystyle\boldsymbol{\beta}^{\text{e}}bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT=ψe𝒃=λ[𝒃:𝜺+e+𝑰:𝜺e]𝜺+e2μ𝜺+e𝒃𝜺+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{,}= - divide start_ARG ∂ italic_ψ start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG start_ARG ∂ bold_italic_b end_ARG = - italic_λ [ bold_italic_b : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT + bold_italic_I : bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - 2 italic_μ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ⋅ bold_italic_b ⋅ bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ,(11)
𝜷psuperscript𝜷p\displaystyle\boldsymbol{\beta}^{\text{p}}bold_italic_β start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT=ψp𝒃=Ha𝒂𝒃𝒂Hk[𝒃:𝒌]𝒌,\displaystyle=-\dfrac{\partial\psi^{\text{p}}}{\partial\boldsymbol{b}}=-H_{%\mathrm{a}}\,\boldsymbol{a}\cdot\boldsymbol{b}\cdot\boldsymbol{a}-H_{\mathrm{k%}}\,\left[\boldsymbol{b}:\boldsymbol{k}\right]\,\boldsymbol{k}\,\text{,}= - divide start_ARG ∂ italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT end_ARG start_ARG ∂ bold_italic_b end_ARG = - italic_H start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT bold_italic_a ⋅ bold_italic_b ⋅ bold_italic_a - italic_H start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT [ bold_italic_b : bold_italic_k ] bold_italic_k ,(12)
𝜷𝜷\displaystyle\boldsymbol{\beta}\phantom{{}^{\text{e}}}bold_italic_β=𝜷e+𝜷p,absentsuperscript𝜷esuperscript𝜷p,\displaystyle=\boldsymbol{\beta}^{\text{e}}+\boldsymbol{\beta}^{\text{p}}\,%\text{,}= bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT + bold_italic_β start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT ,(13)

where 𝝈𝝈\boldsymbol{\sigma}bold_italic_σ is the stress tensor, 𝜶𝜶\boldsymbol{\alpha}bold_italic_α is the back stress tensor, 𝜿𝜿\boldsymbol{\kappa}bold_italic_κ is the drag stress tensor and 𝜷𝜷\boldsymbol{\beta}bold_italic_β the energy release rate. Using the thermodynamic driving forces, the dissipation inequality reads in its reduced form

𝒟redsuperscript𝒟red\displaystyle\mathcal{D}^{\text{red}}caligraphic_D start_POSTSUPERSCRIPT red end_POSTSUPERSCRIPT=𝝈:𝜺˙ψ˙=𝝈:𝜺˙p+𝜶:𝒂˙+𝜿:𝒌˙+𝜷:𝒃˙0.:absent𝝈˙𝜺˙𝜓𝝈:superscript˙𝜺p𝜶:˙𝒂𝜿:˙𝒌𝜷:˙𝒃0\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.= bold_italic_σ : over˙ start_ARG bold_italic_ε end_ARG - over˙ start_ARG italic_ψ end_ARG = bold_italic_σ : over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT + bold_italic_α : over˙ start_ARG bold_italic_a end_ARG + bold_italic_κ : over˙ start_ARG bold_italic_k end_ARG + bold_italic_β : over˙ start_ARG bold_italic_b end_ARG ≥ 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

g𝑔\displaystyle gitalic_g=Φ(𝝈,𝜶,𝜿,𝒃)+Γα(𝜶,𝒃)+Γβ(𝜷e,𝒃)absentΦ𝝈𝜶𝜿𝒃subscriptΓ𝛼𝜶𝒃subscriptΓ𝛽superscript𝜷e𝒃\displaystyle=\Phi(\boldsymbol{\sigma},\,\boldsymbol{\alpha},\,\boldsymbol{%\kappa,\,\boldsymbol{b}})+\Gamma_{\alpha}(\boldsymbol{\alpha},\,\boldsymbol{b}%)+\Gamma_{\beta}(\boldsymbol{\beta}^{\text{e}},\,\boldsymbol{b})= roman_Φ ( bold_italic_σ , bold_italic_α , bold_italic_κ bold_, bold_italic_b ) + roman_Γ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ( bold_italic_α , bold_italic_b ) + roman_Γ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT ( bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT , bold_italic_b )(15)

and consists of the yield function ΦΦ\Phiroman_Φ and non-associated parts ΓαsubscriptΓ𝛼\Gamma_{\!\alpha}roman_Γ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT and ΓβsubscriptΓ𝛽\Gamma_{\beta}roman_Γ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT.The yield function is chosen as

ΦΦ\displaystyle\Phiroman_Φ=τ¯τy13𝒃1:𝜿Δτ[1exp(|𝒃1:𝜿|κu)],\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{,}= square-root start_ARG over¯ start_ARG italic_τ end_ARG end_ARG - italic_τ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT - divide start_ARG 1 end_ARG start_ARG 3 end_ARG bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT : bold_italic_κ - roman_Δ italic_τ [ 1 - roman_exp ( - divide start_ARG | bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT : bold_italic_κ | end_ARG start_ARG italic_κ start_POSTSUBSCRIPT roman_u end_POSTSUBSCRIPT end_ARG ) ] ,(16)
τ¯¯𝜏\displaystyle\bar{\tau}over¯ start_ARG italic_τ end_ARG=32𝒃1:[𝝉𝒃1𝝉]12[𝒃1:𝝉]2with𝝉=𝝈𝜶,\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{,}= divide start_ARG 3 end_ARG start_ARG 2 end_ARG bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT : [ bold_italic_τ ⋅ bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ⋅ bold_italic_τ ] - divide start_ARG 1 end_ARG start_ARG 2 end_ARG [ bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT : bold_italic_τ ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT with bold_italic_τ = bold_italic_σ - bold_italic_α ,(17)

where the equivalent stress τ¯¯𝜏\sqrt{\bar{\tau}}square-root start_ARG over¯ start_ARG italic_τ end_ARG end_ARG is of von Mises-type and depends on the relative stress tensor 𝝉𝝉\boldsymbol{\tau}bold_italic_τ.According to Eq.(16), the yield function captures linear and exponential isotropic hardening. The respective model parameters are denoted as τysubscript𝜏𝑦\tau_{y}italic_τ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT, ΔτΔ𝜏\Delta\tauroman_Δ italic_τ and κusubscript𝜅u\kappa_{\mathrm{u}}italic_κ start_POSTSUBSCRIPT roman_u end_POSTSUBSCRIPT.The non-associative parts of potential(15) are given as

ΓαsubscriptΓ𝛼\displaystyle\Gamma_{\alpha}roman_Γ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT=Ba2Ha𝒃1:[𝜶𝒃1𝜶],:absentsubscript𝐵𝑎2subscript𝐻𝑎superscript𝒃1delimited-[]𝜶superscript𝒃1𝜶,\displaystyle=\dfrac{B_{a}}{2\,H_{a}}\,\boldsymbol{b}^{-1}:\left[\boldsymbol{%\alpha}\cdot\boldsymbol{b}^{-1}\cdot\boldsymbol{\alpha}\right]\,\text{,}= divide start_ARG italic_B start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT end_ARG start_ARG 2 italic_H start_POSTSUBSCRIPT italic_a end_POSTSUBSCRIPT end_ARG bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT : [ bold_italic_α ⋅ bold_italic_b start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ⋅ bold_italic_α ] ,(18)
ΓβsubscriptΓ𝛽\displaystyle\Gamma_{\beta}roman_Γ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT=η12[𝒃m:𝜷e]2+η22𝒃m:[𝜷e𝒃m𝜷e]\displaystyle=\dfrac{\eta_{1}}{2}\,\left[\boldsymbol{b}^{m}:\boldsymbol{\beta}%^{\text{e}}\right]^{2}+\dfrac{\eta_{2}}{2}\,\boldsymbol{b}^{m}:\left[%\boldsymbol{\beta}^{\text{e}}\cdot\boldsymbol{b}^{m}\cdot\boldsymbol{\beta}^{%\text{e}}\right]= divide start_ARG italic_η start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG [ bold_italic_b start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT : bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + divide start_ARG italic_η start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_ARG start_ARG 2 end_ARG bold_italic_b start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT : [ bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ⋅ bold_italic_b start_POSTSUPERSCRIPT italic_m end_POSTSUPERSCRIPT ⋅ bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ](19)

and extend the model to non-linear kinematic hardening of Armstrong-Frederick type by means of model parameter Basubscript𝐵aB_{\mathrm{a}}italic_B start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT, cf.[35].The damage evolution is controlled by potential ΓβsubscriptΓ𝛽\Gamma_{\beta}roman_Γ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT and consists of an isotropic part (damage modulus η1subscript𝜂1\eta_{1}italic_η start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT) and an anisotropic part (damage modulus η2subscript𝜂2\eta_{2}italic_η start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT).The damage exponent m𝑚mitalic_m allows for greater flexibility when calibrating the model to experiments. Furthermore, potential ΓβsubscriptΓ𝛽\Gamma_{\beta}roman_Γ start_POSTSUBSCRIPT italic_β end_POSTSUBSCRIPT depends only on the elastic part of the energy release rate 𝜷𝜷\boldsymbol{\beta}bold_italic_β. 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

𝜺˙p=λ˙g𝝈,𝒂˙=λ˙g𝜶,𝒌˙=λ˙g𝜿,𝒃˙=λ˙g𝜷=λ˙g𝜷e,formulae-sequencesuperscript˙𝜺p˙𝜆𝑔𝝈formulae-sequence˙𝒂˙𝜆𝑔𝜶formulae-sequence˙𝒌˙𝜆𝑔𝜿˙𝒃˙𝜆𝑔𝜷˙𝜆𝑔superscript𝜷e,\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{,}over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT = over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_σ end_ARG , over˙ start_ARG bold_italic_a end_ARG = over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_α end_ARG , over˙ start_ARG bold_italic_k end_ARG = over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_κ end_ARG , over˙ start_ARG bold_italic_b end_ARG = over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_β end_ARG = over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_β start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG ,(20)

which concludes the constitutive model.

Isotropic model

The isotropic version of the anisotropic ECC-Model follows straightforward by setting η2=0subscript𝜂20\eta_{2}=0italic_η start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT = 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}caligraphic_G is specified as

𝒢𝒢\displaystyle\mathcal{G}caligraphic_G=𝒢e+𝝈:𝜺pψp,:absentsuperscript𝒢e𝝈superscript𝜺psuperscript𝜓p,\displaystyle=\mathcal{G}^{\text{e}}+\boldsymbol{\sigma}:\boldsymbol{%\varepsilon}^{\text{p}}-\psi^{\text{p}}\,\text{,}= caligraphic_G start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT + bold_italic_σ : bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT - italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT ,(21)
𝒢esuperscript𝒢e\displaystyle\mathcal{G}^{\text{e}}caligraphic_G start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT=1+2ν2E[𝑯:[𝝈+dev𝑯𝝈+dev]+𝝈dev:𝝈dev]+12ν2E[σhyd21ωDhyd+σhyd2],\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-\omega\,D^{\mathrm{hyd}}}+\left%\langle-\sigma^{\mathrm{hyd}}\right\rangle^{2}\right]\,\text{,}= divide start_ARG 1 + 2 italic_ν end_ARG start_ARG 2 italic_E end_ARG [ bold_italic_H : [ bold_italic_σ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ⋅ bold_italic_H ⋅ bold_italic_σ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ] + bold_italic_σ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT : bold_italic_σ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ] + divide start_ARG 1 - 2 italic_ν end_ARG start_ARG 2 italic_E end_ARG [ divide start_ARG ⟨ italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT end_ARG + ⟨ - italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ] ,(22)
ψpsuperscript𝜓p\displaystyle\psi^{\text{p}}italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT=12K𝒂:𝒂+Rinf[p+exp(cp)c]:absent12𝐾𝒂𝒂subscript𝑅infdelimited-[]𝑝exp𝑐𝑝𝑐\displaystyle=\frac{1}{2}\,K\,\boldsymbol{a}:\boldsymbol{a}+R_{\mathrm{inf}}\,%\left[p+\frac{\mathrm{exp}\left(-c\,p\right)}{c}\right]\,= divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_K bold_italic_a : bold_italic_a + italic_R start_POSTSUBSCRIPT roman_inf end_POSTSUBSCRIPT [ italic_p + divide start_ARG roman_exp ( - italic_c italic_p ) end_ARG start_ARG italic_c end_ARG ](23)

and depends on model parameters Young’s modulus E𝐸Eitalic_E, Poisson’s ratio ν𝜈\nuitalic_ν, damage parameter ω𝜔\omegaitalic_ω, the isotropic hardening parameters Rinfsubscript𝑅infR_{\mathrm{inf}}italic_R start_POSTSUBSCRIPT roman_inf end_POSTSUBSCRIPT and c𝑐citalic_c, and the kinematic hardening modulus K𝐾Kitalic_K. The additive decomposition of the total strains into elastic and plastic parts is again considered as 𝜺=𝜺e+𝜺p𝜺superscript𝜺esuperscript𝜺p\boldsymbol{\varepsilon}=\boldsymbol{\varepsilon}^{\text{e}}+\boldsymbol{%\varepsilon}^{\text{p}}bold_italic_ε = bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT + bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT. Furthermore, the stress tensor is decomposed into a deviatoric (𝝈devsuperscript𝝈dev\boldsymbol{\sigma}^{\mathrm{dev}}bold_italic_σ start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT) and a hydrostatic part (σhyd𝑰superscript𝜎hyd𝑰\sigma^{\mathrm{hyd}}\,\boldsymbol{I}italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT bold_italic_I) and subsequently into their respective positive and negative parts in order to capture the MCR-effect. Strain-like variables 𝒂𝒂\boldsymbol{a}bold_italic_a and p𝑝pitalic_p are suitable for kinematic and isotropic hardening, respectively. The damage evolution is captured by means of second-order tensor 𝑫𝑫\boldsymbol{D}bold_italic_D, which enters energy(22) in terms of variables

Dhydsuperscript𝐷hyd\displaystyle D^{\mathrm{hyd}}italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT13tr(𝑫),𝑯[𝑰𝑫]12.formulae-sequenceabsent13tr𝑫,𝑯superscriptdelimited-[]𝑰𝑫12.\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{.}≔ divide start_ARG 1 end_ARG start_ARG 3 end_ARG roman_tr ( bold_italic_D ) , bold_italic_H ≔ [ bold_italic_I - bold_italic_D ] start_POSTSUPERSCRIPT - divide start_ARG 1 end_ARG start_ARG 2 end_ARG end_POSTSUPERSCRIPT .

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

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

𝜺esuperscript𝜺e\displaystyle\phantom{-}\boldsymbol{\varepsilon}^{\text{e}}bold_italic_ε start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT=𝒢𝝈=1+νE[dev(𝑯𝝈+dev𝑯)+𝝈dev]+12νE[σhyd1ωDhydσhyd]𝑰,absent𝒢𝝈1𝜈𝐸delimited-[]dev𝑯superscriptsubscript𝝈dev𝑯superscriptsubscript𝝈dev12𝜈𝐸delimited-[]delimited-⟨⟩superscript𝜎hyd1𝜔superscript𝐷hyddelimited-⟨⟩superscript𝜎hyd𝑰,\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-\omega\,D^{\mathrm{hyd}}}-\left\langle-%\sigma^{\mathrm{hyd}}\right\rangle\right]\,\boldsymbol{I}\,\text{,}= divide start_ARG ∂ caligraphic_G end_ARG start_ARG ∂ bold_italic_σ end_ARG = divide start_ARG 1 + italic_ν end_ARG start_ARG italic_E end_ARG [ roman_dev ( bold_italic_H ⋅ bold_italic_σ start_POSTSUBSCRIPT + end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ⋅ bold_italic_H ) + bold_italic_σ start_POSTSUBSCRIPT - end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ] + divide start_ARG 1 - 2 italic_ν end_ARG start_ARG italic_E end_ARG [ divide start_ARG ⟨ italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ end_ARG start_ARG 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT end_ARG - ⟨ - italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ ] bold_italic_I ,(24)
𝜶𝜶\displaystyle\phantom{-}\boldsymbol{\alpha}bold_italic_α=𝒢𝒂=23K𝒂,absent𝒢𝒂23𝐾𝒂,\displaystyle=-\dfrac{\partial\mathcal{G}}{\partial\boldsymbol{a}}=\frac{2}{3}%\,K\,\boldsymbol{a}\,\text{,}= - divide start_ARG ∂ caligraphic_G end_ARG start_ARG ∂ bold_italic_a end_ARG = divide start_ARG 2 end_ARG start_ARG 3 end_ARG italic_K bold_italic_a ,(25)
q𝑞\displaystyle\phantom{-}qitalic_q=𝒢p=Rinf[1exp(cp)],absent𝒢𝑝subscript𝑅infdelimited-[]1exp𝑐𝑝,\displaystyle=-\dfrac{\partial\mathcal{G}}{\partial p}=R_{\mathrm{inf}}\,\left%[1-\mathrm{exp}\left(-c\,p\right)\right]\,\text{,}= - divide start_ARG ∂ caligraphic_G end_ARG start_ARG ∂ italic_p end_ARG = italic_R start_POSTSUBSCRIPT roman_inf end_POSTSUBSCRIPT [ 1 - roman_exp ( - italic_c italic_p ) ] ,(26)
𝒀𝒀\displaystyle\phantom{-}\boldsymbol{Y}bold_italic_Y=𝒢𝑫=1+2νE[𝝈+dev𝑯𝝈+dev]:𝑯𝑫ω12ν6Eσhyd2[1ωDhyd]2𝑰.:absent𝒢𝑫12𝜈𝐸delimited-[]subscriptsuperscript𝝈dev𝑯subscriptsuperscript𝝈dev𝑯𝑫𝜔12𝜈6𝐸superscriptdelimited-⟨⟩superscript𝜎hyd2superscriptdelimited-[]1𝜔superscript𝐷hyd2𝑰.\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}}-\omega\,\frac{1-2\,\nu}{6\,E}\,\frac{\left\langle%\sigma^{\mathrm{hyd}}\right\rangle^{2}}{\left[1-\omega\,D^{\mathrm{hyd}}\right%]^{2}}\,\boldsymbol{I}\,\text{.}= - divide start_ARG ∂ caligraphic_G end_ARG start_ARG ∂ bold_italic_D end_ARG = - divide start_ARG 1 + 2 italic_ν end_ARG start_ARG italic_E end_ARG [ bold_italic_σ start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ⋅ bold_italic_H ⋅ bold_italic_σ start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ] : divide start_ARG ∂ bold_italic_H end_ARG start_ARG ∂ bold_italic_D end_ARG - italic_ω divide start_ARG 1 - 2 italic_ν end_ARG start_ARG 6 italic_E end_ARG divide start_ARG ⟨ italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG [ 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG bold_italic_I .(27)

The transformation of Gibbs energy 𝒢𝒢\mathcal{G}caligraphic_G into Helmholtz energy ψ𝜓\psiitalic_ψ reads

ψ𝜓\displaystyle\psiitalic_ψ=𝝈:𝜺𝒢=𝝈:𝜺𝒢e𝝈:𝜺p+ψp:absent𝝈𝜺𝒢𝝈:𝜺superscript𝒢e𝝈:superscript𝜺psuperscript𝜓p\displaystyle=\boldsymbol{\sigma}:\boldsymbol{\varepsilon}-\mathcal{G}=%\boldsymbol{\sigma}:\boldsymbol{\varepsilon}-\mathcal{G}^{\text{e}}-%\boldsymbol{\sigma}:\boldsymbol{\varepsilon}^{\text{p}}+\psi^{\text{p}}\,= bold_italic_σ : bold_italic_ε - caligraphic_G = bold_italic_σ : bold_italic_ε - caligraphic_G start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT - bold_italic_σ : bold_italic_ε start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT + italic_ψ start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT(28)

and allows to derive the reduced dissipation inequality as

𝒟redsuperscript𝒟red\displaystyle\mathcal{D}^{\mathrm{red}}caligraphic_D start_POSTSUPERSCRIPT roman_red end_POSTSUPERSCRIPT=𝝈:𝜺˙ψ˙=𝝈:𝜺˙p𝒀:𝑫˙𝜶:𝒂˙qp˙0.:absent𝝈˙𝜺˙𝜓𝝈:superscript˙𝜺p𝒀:˙𝑫𝜶:˙𝒂𝑞˙𝑝0.\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}}-q\,\dot{p}\geq 0%\,\text{.}= bold_italic_σ : over˙ start_ARG bold_italic_ε end_ARG - over˙ start_ARG italic_ψ end_ARG = bold_italic_σ : over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT - bold_italic_Y : over˙ start_ARG bold_italic_D end_ARG - bold_italic_α : over˙ start_ARG bold_italic_a end_ARG - italic_q over˙ start_ARG italic_p end_ARG ≥ 0 .(29)

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

g𝑔\displaystyle gitalic_g=Φ(𝝈,𝜶,𝑫,q)+Γα(𝜶)absentΦ𝝈𝜶𝑫𝑞subscriptΓ𝛼𝜶\displaystyle=\Phi\!\left({\boldsymbol{\sigma},\boldsymbol{\alpha},\boldsymbol%{D},q}\right)+\Gamma_{\!\alpha}\!\left({\boldsymbol{\alpha}}\right)= roman_Φ ( bold_italic_σ , bold_italic_α , bold_italic_D , italic_q ) + roman_Γ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ( bold_italic_α )(30)

with von Mises-type yield function

ΦΦ\displaystyle\Phiroman_Φ=τ¯τyqwithτ¯=32𝝉:𝝉,𝝉=𝑯𝝈dev𝑯𝜶:formulae-sequenceabsent¯𝜏subscript𝜏y𝑞with¯𝜏32𝝉𝝉,𝝉𝑯superscript𝝈dev𝑯𝜶\displaystyle=\sqrt{\bar{\tau}}-\tau_{\mathrm{y}}-q\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}= square-root start_ARG over¯ start_ARG italic_τ end_ARG end_ARG - italic_τ start_POSTSUBSCRIPT roman_y end_POSTSUBSCRIPT - italic_q with over¯ start_ARG italic_τ end_ARG = divide start_ARG 3 end_ARG start_ARG 2 end_ARG bold_italic_τ : bold_italic_τ , bold_italic_τ = bold_italic_H ⋅ bold_italic_σ start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT ⋅ bold_italic_H - bold_italic_α(31)

and initial yield limit τysubscript𝜏𝑦\tau_{y}italic_τ start_POSTSUBSCRIPT italic_y end_POSTSUBSCRIPT. The non-associative part of potential(30) reads

ΓαsubscriptΓ𝛼\displaystyle\Gamma_{\!\alpha}roman_Γ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT=3γ4K𝜶:𝜶:absent3𝛾4𝐾𝜶𝜶\displaystyle=\frac{3\,\gamma}{4\,K}\,\boldsymbol{\alpha}:\boldsymbol{\alpha}= divide start_ARG 3 italic_γ end_ARG start_ARG 4 italic_K end_ARG bold_italic_α : bold_italic_α(32)

and extends the model to non-linear kinematic hardening due to model parameter γ𝛾\gammaitalic_γ. The evolution equations then result in

𝜺˙psuperscript˙𝜺p\displaystyle\dot{\boldsymbol{\varepsilon}}^{\text{p}}over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT=λ˙g𝝈,𝒂˙=λ˙g𝜶,p˙=λ˙gq=λ˙.formulae-sequenceabsent˙𝜆𝑔𝝈,formulae-sequence˙𝒂˙𝜆𝑔𝜶,˙𝑝˙𝜆𝑔𝑞˙𝜆.\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{p}=-\dot{\lambda}\,\dfrac{\partial g}%{\partial q}=\dot{\lambda}\,\text{.}= over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_σ end_ARG , over˙ start_ARG bold_italic_a end_ARG = - over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ bold_italic_α end_ARG , over˙ start_ARG italic_p end_ARG = - over˙ start_ARG italic_λ end_ARG divide start_ARG ∂ italic_g end_ARG start_ARG ∂ italic_q end_ARG = over˙ start_ARG italic_λ end_ARG .(33)

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

𝑫˙˙𝑫\displaystyle\dot{\boldsymbol{D}}over˙ start_ARG bold_italic_D end_ARG=[PY¯]s[𝜺˙+p𝜺˙p],absentsuperscriptdelimited-[]𝑃¯𝑌𝑠delimited-[]subscriptsuperscript˙𝜺psubscriptsuperscript˙𝜺p,\displaystyle=\left[P\,\bar{Y}\right]^{s}\,\left[\dot{\boldsymbol{\varepsilon}%}^{\text{p}}_{+}-\dot{\boldsymbol{\varepsilon}}^{\text{p}}_{-}\right]\,\text{,}= [ italic_P over¯ start_ARG italic_Y end_ARG ] start_POSTSUPERSCRIPT italic_s end_POSTSUPERSCRIPT [ over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT - over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - end_POSTSUBSCRIPT ] ,(34)
withY¯with¯𝑌\displaystyle\text{with }\bar{Y}with over¯ start_ARG italic_Y end_ARG=1+ν2Etr([𝑯𝝈+dev𝑯]2)+36ν2Eωσhyd2[1ωDhyd]2absent1𝜈2𝐸trsuperscriptdelimited-[]𝑯subscriptsuperscript𝝈dev𝑯236𝜈2𝐸𝜔superscriptdelimited-⟨⟩superscript𝜎hyd2superscriptdelimited-[]1𝜔superscript𝐷hyd2\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}\,\omega\,\frac{\left\langle{\sigma^{\mathrm{hyd}}}\right%\rangle^{2}}{\left[1-\omega\,D^{\mathrm{hyd}}\right]^{2}}= divide start_ARG 1 + italic_ν end_ARG start_ARG 2 italic_E end_ARG roman_tr ( [ bold_italic_H ⋅ bold_italic_σ start_POSTSUPERSCRIPT roman_dev end_POSTSUPERSCRIPT start_POSTSUBSCRIPT + end_POSTSUBSCRIPT ⋅ bold_italic_H ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ) + divide start_ARG 3 - 6 italic_ν end_ARG start_ARG 2 italic_E end_ARG italic_ω divide start_ARG ⟨ italic_σ start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ⟩ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG [ 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG(35)

and depends on model parameters P𝑃Pitalic_P and s𝑠sitalic_s. This evolution is formulated in terms of a scalar-valued energy release rate Y¯¯𝑌\bar{Y}over¯ start_ARG italic_Y end_ARG and the additive decomposition of the plastic strain rate tensor 𝜺˙psuperscript˙𝜺p\dot{\boldsymbol{\varepsilon}}^{\text{p}}over˙ start_ARG bold_italic_ε end_ARG start_POSTSUPERSCRIPT p end_POSTSUPERSCRIPT 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}}over˙ start_ARG bold_italic_D end_ARG=[P^Y¯]s^ε˙p,eq𝑰absentsuperscriptdelimited-[]^𝑃¯𝑌^𝑠superscript˙𝜀peq𝑰\displaystyle=\left[\hat{P}\,\bar{Y}\right]^{\hat{s}}\,\dot{\varepsilon}^{%\mathrm{p,eq}}\,\boldsymbol{I}= [ over^ start_ARG italic_P end_ARG over¯ start_ARG italic_Y end_ARG ] start_POSTSUPERSCRIPT over^ start_ARG italic_s end_ARG end_POSTSUPERSCRIPT over˙ start_ARG italic_ε end_ARG start_POSTSUPERSCRIPT roman_p , roman_eq end_POSTSUPERSCRIPT bold_italic_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 ε110.13subscript𝜀110.13\varepsilon_{11}\!\approx 0.13italic_ε start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ≈ 0.13, which represents a typical material used in bulk forming such as forward extrusion. The strains remain spatially hom*ogeneous 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 λ𝜆\lambdaitalic_λ and μ𝜇\muitalic_μ and the same initial yield stress τysubscript𝜏y\tau_{\mathrm{y}}italic_τ start_POSTSUBSCRIPT roman_y end_POSTSUBSCRIPT. 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 nameSymbolValueUnit
First Lame parameterλ𝜆\lambdaitalic_λ118870118870118870118870MPamegapascal\mathrm{MPa}roman_MPa
Second Lame parameterμ𝜇\muitalic_μ79249792497924979249MPamegapascal\mathrm{MPa}roman_MPa
Yield stressτysubscript𝜏y\tau_{\mathrm{y}}italic_τ start_POSTSUBSCRIPT roman_y end_POSTSUBSCRIPT308.260308.260308.260308.260MPamegapascal\mathrm{MPa}roman_MPa
Kinematic hardening modulusHasubscript𝐻aH_{\mathrm{a}}italic_H start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT7728.8637728.8637728.8637728.863MPamegapascal\mathrm{MPa}roman_MPa
Armstrong-Frederick extensionBasubscript𝐵aB_{\mathrm{a}}italic_B start_POSTSUBSCRIPT roman_a end_POSTSUBSCRIPT38.21838.21838.21838.218-
Isotropic hardening modulusHksubscript𝐻kH_{\mathrm{k}}italic_H start_POSTSUBSCRIPT roman_k end_POSTSUBSCRIPT1.8291041.829superscript1041.829\cdot 10^{-4}1.829 ⋅ 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPTMPamegapascal\mathrm{MPa}roman_MPa
Isotropic hardening incrementΔτyΔsubscript𝜏y\Delta\tau_{\mathrm{y}}roman_Δ italic_τ start_POSTSUBSCRIPT roman_y end_POSTSUBSCRIPT2.2611022.261superscript1022.261\cdot 10^{-2}2.261 ⋅ 10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPTMPamegapascal\mathrm{MPa}roman_MPa
Isotropic hardening parameterκusubscript𝜅u\kappa_{\mathrm{u}}italic_κ start_POSTSUBSCRIPT roman_u end_POSTSUBSCRIPT0.1590.1590.1590.159MPamegapascal\mathrm{MPa}roman_MPa
Anisotropic (ECC-a model)
Isotropic damage modulusη1subscript𝜂1\eta_{1}italic_η start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT0.00.00.00.0MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Anisotropic damage modulusη2subscript𝜂2\eta_{2}italic_η start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT14.50314.50314.50314.503MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Damage exponentm𝑚mitalic_m11.21711.21711.21711.217-
Isotropic (ECC-i model)
Isotropic damage modulusη^1subscript^𝜂1\hat{\eta}_{1}over^ start_ARG italic_η end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT14.40814.40814.40814.408MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Anisotropic damage modulusη^2subscript^𝜂2\hat{\eta}_{2}over^ start_ARG italic_η end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT0.00.00.00.0MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Damage exponentialm^^𝑚\hat{m}over^ start_ARG italic_m end_ARG11.37311.37311.37311.373-
Parameter nameSymbolValueUnit
First Lame parameterλ𝜆\lambdaitalic_λ118875.0118875.0118875.0118875.0MPamegapascal\mathrm{MPa}roman_MPa
Second Lame parameterμ𝜇\muitalic_μ79250.079250.079250.079250.0MPamegapascal\mathrm{MPa}roman_MPa
Yield stressτysubscript𝜏y\tau_{\mathrm{y}}italic_τ start_POSTSUBSCRIPT roman_y end_POSTSUBSCRIPT308.26308.26308.26308.26MPamegapascal\mathrm{MPa}roman_MPa
Isotropic hardening parameterc𝑐citalic_c0.18540.18540.18540.1854-
Kinematic hardening modulusK𝐾Kitalic_K42155.742155.742155.742155.7MPamegapascal\mathrm{MPa}roman_MPa
Nonlinear kinematic parameterγ𝛾\gammaitalic_γ7133.57133.57133.57133.5MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Saturation hardening stressRinfsubscript𝑅infR_{\mathrm{inf}}italic_R start_POSTSUBSCRIPT roman_inf end_POSTSUBSCRIPT26452.026452.026452.026452.0MPamegapascal\mathrm{MPa}roman_MPa
Anisotropic (LEM-a model)
Damage parameterω𝜔\omegaitalic_ω1.01.01.01.0-
Damage scaleP𝑃Pitalic_P1256.71256.71256.71256.7MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Damage exponents𝑠sitalic_s0.20.20.20.2-
Isotropic (LEM-i model)
Damage parameterω^^𝜔\hat{\omega}over^ start_ARG italic_ω end_ARG1.01.01.01.0-
Damage scaleP^^𝑃\hat{P}over^ start_ARG italic_P end_ARG623.1623.1623.1623.1MPa1superscriptMPa1\mathrm{M}\mathrm{Pa}^{-1}roman_MPa start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT
Damage exponents^^𝑠\hat{s}over^ start_ARG italic_s end_ARG0.20.20.20.2-

\psfrag{tensx}[c][c]{\scalebox{0.7}{Tensile strain $\varepsilon_{11}$ [-]}}\psfrag{tensy}[c][c]{\scalebox{0.7}{Tensile stress $\sigma_{11}$ [$\mathrm{MPa}$]}}\psfrag{exp}[l][l]{\scalebox{0.7}{Experiment}}\psfrag{modelekh}[l][l]{\scalebox{0.7}{\text{ECC-a} model}}\psfrag{modeleki}[l][l]{\scalebox{0.7}{\text{ECC-i} model}}\includegraphics[width=346.89731pt]{Chapter/02_Fundamentals/Subsections/Figures/TensFit1.eps}

\psfrag{tensx}[c][c]{\scalebox{0.7}{Tensile strain $\varepsilon_{11}$ [-]}}\psfrag{tensy}[c][c]{\scalebox{0.7}{Tensile stress $\sigma_{11}$ [$\mathrm{MPa}$]}}\psfrag{exp}[l][l]{\scalebox{0.7}{Experiment}}\psfrag{modellad}[l][l]{\scalebox{0.7}{\text{LEM-a} model}}\psfrag{modellid}[l][l]{\scalebox{0.7}{\text{LEM-i} model}}\includegraphics[width=346.89731pt]{Chapter/02_Fundamentals/Subsections/Figures/TensFit2.eps}

\psfrag{shearx}[c][c]{\scalebox{0.7}{Shear strain $\varepsilon_{12}$ [-]}}\psfrag{sheary}[c][c]{\scalebox{0.7}{Shear stress $\sigma_{12}$ [$\mathrm{MPa}$]}}\psfrag{modelekh}[l][l]{\scalebox{0.7}{\text{ECC-a} model}}\psfrag{modeleki}[l][l]{\scalebox{0.7}{\text{ECC-i} model}}\psfrag{modellad}[l][l]{\scalebox{0.7}{\text{LEM-a} model}}\psfrag{modellid}[l][l]{\scalebox{0.7}{\text{LEM-i} model}}\includegraphics[width=346.89731pt]{Chapter/02_Fundamentals/Subsections/Figures/TensFitShear.eps}

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 ε11=0.07subscript𝜀110.07\varepsilon_{11}=0.07italic_ε start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT = 0.07. However, the respective error is below 10MPatimes10megapascal10\text{\,}\mathrm{MPa}start_ARG 10 end_ARG start_ARG times end_ARG start_ARG roman_MPa end_ARG and hence, less than 2.5%percent2.52.5\%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%percent88\%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:

𝜺=[ε11ε12ε13symε22ε23symsymε33][𝒆1,𝒆2,𝒆3]𝝈=[σ11σ12σ13symσ22σ23symsymσ33][𝒆1,𝒆2,𝒆3]formulae-sequence𝜺subscriptmatrixsubscript𝜀11missing-subexpressionsubscript𝜀12missing-subexpressionsubscript𝜀13symmissing-subexpressionsubscript𝜀22missing-subexpressionsubscript𝜀23symmissing-subexpressionsymmissing-subexpressionsubscript𝜀33subscript𝒆1subscript𝒆2subscript𝒆3𝝈subscriptmatrixsubscript𝜎11missing-subexpressionsubscript𝜎12missing-subexpressionsubscript𝜎13symmissing-subexpressionsubscript𝜎22missing-subexpressionsubscript𝜎23symmissing-subexpressionsymmissing-subexpressionsubscript𝜎33subscript𝒆1subscript𝒆2subscript𝒆3\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]}bold_italic_ε = [ start_ARG start_ROW start_CELL italic_ε start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT 13 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT 23 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT 33 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT bold_italic_σ = [ start_ARG start_ROW start_CELL italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 12 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 13 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 23 end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 33 end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT(37)

In Eq.(37), the components are defined with respect to a chosen basis 𝒆isubscript𝒆𝑖\boldsymbol{e}_{i}bold_italic_e start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT (here, a cartesian setting is adopted). Loading is implemented by enforcing a strain in 11111111-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 σij=0subscript𝜎𝑖𝑗0\sigma_{ij}=0italic_σ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT = 0 except for σ11subscript𝜎11\sigma_{11}italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT.

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

Lpwl(t;𝑷)subscript𝐿pwl𝑡𝑷\displaystyle L_{\mathrm{pwl}}\!\left({t;\boldsymbol{P}}\right)italic_L start_POSTSUBSCRIPT roman_pwl end_POSTSUBSCRIPT ( italic_t ; bold_italic_P )=Pi+PiPjtitjtift[ti,tj]formulae-sequenceabsentsubscript𝑃𝑖subscript𝑃𝑖subscript𝑃𝑗subscript𝑡𝑖subscript𝑡𝑗𝑡if𝑡subscript𝑡𝑖subscript𝑡𝑗\displaystyle=P_{i}+\frac{P_{i}-P_{j}}{t_{i}-t_{j}}\,t\quad\text{if }t\in\left%[t_{i},\,t_{j}\right]= italic_P start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT + divide start_ARG italic_P start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_P start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG start_ARG italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT - italic_t start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT end_ARG italic_t if italic_t ∈ [ italic_t start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT , italic_t start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ](38)

or by Bézier curves defined as

LBez(t;𝑸)subscript𝐿Bez𝑡𝑸\displaystyle L_{\mathrm{Bez}}\!\left({t;\boldsymbol{Q}}\right)italic_L start_POSTSUBSCRIPT roman_Bez end_POSTSUBSCRIPT ( italic_t ; bold_italic_Q )=i=0N(Ni)[1t][Ni]tiQiabsentsuperscriptsubscript𝑖0𝑁binomial𝑁𝑖superscriptdelimited-[]1𝑡delimited-[]𝑁𝑖superscript𝑡𝑖subscript𝑄𝑖\displaystyle=\sum_{i=0}^{N}{N\choose i}\,\left[1-t\right]^{\left[N-i\right]}%\,t^{i}\,Q_{i}= ∑ start_POSTSUBSCRIPT italic_i = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_N end_POSTSUPERSCRIPT ( binomial start_ARG italic_N end_ARG start_ARG italic_i end_ARG ) [ 1 - italic_t ] start_POSTSUPERSCRIPT [ italic_N - italic_i ] end_POSTSUPERSCRIPT italic_t start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT italic_Q start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT(39)

with control points Pisubscript𝑃𝑖P_{i}italic_P start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and Qisubscript𝑄𝑖Q_{i}italic_Q start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT, respectively. The time is normalized within the interval t[0;1]𝑡01t\in[0;1]italic_t ∈ [ 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., 𝔼e:e=𝕀:superscript𝔼esuperscripte𝕀\mathbb{E}^{\text{e}}:\mathbb{C}^{\text{e}}=\mathbb{I}blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : blackboard_C start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT = blackboard_I with 𝕀𝕀\mathbb{I}blackboard_I being the symmetric fourth-order identity. Starting from the Helmholtz energy ψ𝜓\psiitalic_ψ, the elastic moduli of the ECC-model are obtained as

[𝔼ECCe]ijkl:=2ψεijεkl=λ[bijbkl]+μ[bikbjl+bilbjk]assignsubscriptdelimited-[]superscriptsubscript𝔼ECCe𝑖𝑗𝑘𝑙superscript2𝜓subscript𝜀𝑖𝑗𝜀𝑘𝑙𝜆delimited-[]subscript𝑏𝑖𝑗subscript𝑏𝑘𝑙𝜇delimited-[]subscript𝑏𝑖𝑘subscript𝑏𝑗𝑙subscript𝑏𝑖𝑙subscript𝑏𝑗𝑘\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][ blackboard_E start_POSTSUBSCRIPT ECC end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ] start_POSTSUBSCRIPT italic_i italic_j italic_k italic_l end_POSTSUBSCRIPT := divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_ψ end_ARG start_ARG ∂ italic_ε start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ∂ italic_ε italic_k italic_l end_ARG = italic_λ [ italic_b start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_b start_POSTSUBSCRIPT italic_k italic_l end_POSTSUBSCRIPT ] + italic_μ [ italic_b start_POSTSUBSCRIPT italic_i italic_k end_POSTSUBSCRIPT italic_b start_POSTSUBSCRIPT italic_j italic_l end_POSTSUBSCRIPT + italic_b start_POSTSUBSCRIPT italic_i italic_l end_POSTSUBSCRIPT italic_b start_POSTSUBSCRIPT italic_j italic_k end_POSTSUBSCRIPT ](40)

with 𝔼ECCe:ECCe=𝕀:superscriptsubscript𝔼ECCesuperscriptsubscriptECCe𝕀\mathbb{E}_{\mathrm{\text{ECC}}}^{\text{e}}:\mathbb{C}_{\mathrm{\text{ECC}}}^{%\text{e}}=\mathbb{I}blackboard_E start_POSTSUBSCRIPT ECC end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : blackboard_C start_POSTSUBSCRIPT ECC end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT = blackboard_I. Likewise, the LEM-model yields the elastic compliance

[LEMe]ijkl:=2𝒢σijσkl=1+νE[12[HikHjl+HilHjk]13[δijHkmHml+HinHnjδkl]+19[HopHop+31ωDhyd]δijδkl]νE11ωDhydδijδklassignsubscriptdelimited-[]superscriptsubscriptLEMe𝑖𝑗𝑘𝑙superscript2𝒢subscript𝜎𝑖𝑗subscript𝜎𝑘𝑙1𝜈𝐸delimited-[]12delimited-[]subscript𝐻𝑖𝑘subscript𝐻𝑗𝑙subscript𝐻𝑖𝑙subscript𝐻𝑗𝑘13delimited-[]subscript𝛿𝑖𝑗subscript𝐻𝑘𝑚subscript𝐻𝑚𝑙subscript𝐻𝑖𝑛subscript𝐻𝑛𝑗subscript𝛿𝑘𝑙19delimited-[]subscript𝐻𝑜𝑝subscript𝐻𝑜𝑝31𝜔superscript𝐷hydsubscript𝛿𝑖𝑗subscript𝛿𝑘𝑙𝜈𝐸11𝜔superscript𝐷hydsubscript𝛿𝑖𝑗subscript𝛿𝑘𝑙\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-\,\omega\,D^{\mathrm{hyd}}}\right]\,\delta_{ij}\,\delta_{kl}%\right]-\frac{\nu}{E}\,\frac{1}{1-\omega\,D^{\mathrm{hyd}}}\,\delta_{ij}\,%\delta_{kl}\end{split}start_ROW start_CELL [ blackboard_C start_POSTSUBSCRIPT LEM end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ] start_POSTSUBSCRIPT italic_i italic_j italic_k italic_l end_POSTSUBSCRIPT := divide start_ARG ∂ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT caligraphic_G end_ARG start_ARG ∂ italic_σ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ∂ italic_σ start_POSTSUBSCRIPT italic_k italic_l end_POSTSUBSCRIPT end_ARG = end_CELL start_CELL divide start_ARG 1 + italic_ν end_ARG start_ARG italic_E end_ARG [ divide start_ARG 1 end_ARG start_ARG 2 end_ARG [ italic_H start_POSTSUBSCRIPT italic_i italic_k end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_j italic_l end_POSTSUBSCRIPT + italic_H start_POSTSUBSCRIPT italic_i italic_l end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_j italic_k end_POSTSUBSCRIPT ] - divide start_ARG 1 end_ARG start_ARG 3 end_ARG [ italic_δ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_k italic_m end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_m italic_l end_POSTSUBSCRIPT + italic_H start_POSTSUBSCRIPT italic_i italic_n end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_n italic_j end_POSTSUBSCRIPT italic_δ start_POSTSUBSCRIPT italic_k italic_l end_POSTSUBSCRIPT ] end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL + divide start_ARG 1 end_ARG start_ARG 9 end_ARG [ italic_H start_POSTSUBSCRIPT italic_o italic_p end_POSTSUBSCRIPT italic_H start_POSTSUBSCRIPT italic_o italic_p end_POSTSUBSCRIPT + divide start_ARG 3 end_ARG start_ARG 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT end_ARG ] italic_δ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_δ start_POSTSUBSCRIPT italic_k italic_l end_POSTSUBSCRIPT ] - divide start_ARG italic_ν end_ARG start_ARG italic_E end_ARG divide start_ARG 1 end_ARG start_ARG 1 - italic_ω italic_D start_POSTSUPERSCRIPT roman_hyd end_POSTSUPERSCRIPT end_ARG italic_δ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT italic_δ start_POSTSUBSCRIPT italic_k italic_l end_POSTSUBSCRIPT end_CELL end_ROW(41)

where δijsubscript𝛿𝑖𝑗\delta_{ij}italic_δ start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT is the Kronecker-Delta and 𝔼LEMe:LEMe=𝕀:superscriptsubscript𝔼LEMesuperscriptsubscriptLEMe𝕀\mathbb{E}_{\mathrm{\text{LEM}}}^{\text{e}}:\mathbb{C}_{\mathrm{\text{LEM}}}^{%\text{e}}=\mathbb{I}blackboard_E start_POSTSUBSCRIPT LEM end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : blackboard_C start_POSTSUBSCRIPT LEM end_POSTSUBSCRIPT start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT = blackboard_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

f𝔼𝔼e::𝔼e\displaystyle f_{\mathbb{E}}\coloneqq\sqrt{\mathbb{E}^{\text{e}}::\mathbb{E}^{%\text{e}}}italic_f start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT ≔ square-root start_ARG blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : : blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG(42)

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

f𝔼diffsuperscriptsubscript𝑓𝔼diff\displaystyle f_{\mathbb{E}}^{\mathrm{diff}}italic_f start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_diff end_POSTSUPERSCRIPT=[𝔼e::𝔼e𝔼¯e::𝔼¯e]2.\displaystyle=\left[{\sqrt{\mathbb{E}^{\text{e}}::\mathbb{E}^{\text{e}}}-\sqrt%{\bar{\mathbb{E}}^{\text{e}}::\bar{\mathbb{E}}^{\text{e}}}}\right]^{2}.= [ square-root start_ARG blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : : blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG - square-root start_ARG over¯ start_ARG blackboard_E end_ARG start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : : over¯ start_ARG blackboard_E end_ARG start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT end_ARG ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT .(43)

The undamaged configuration represents a natural reference solution, i.e., 𝔼¯e=𝔼e(t=0)superscript¯𝔼esuperscript𝔼e𝑡0\bar{\mathbb{E}}^{\text{e}}=\mathbb{E}^{\text{e}}(t=0)over¯ start_ARG blackboard_E end_ARG start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT = blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT ( italic_t = 0 ).As an alternative to the Frobenius-norm(42), the directional elastic stiffness and compliance

E𝒓=[𝒓𝒓]:𝔼e:[𝒓𝒓]:subscript𝐸𝒓delimited-[]tensor-product𝒓𝒓superscript𝔼e:delimited-[]tensor-product𝒓𝒓\displaystyle E_{\boldsymbol{r}}=\left[\boldsymbol{r}\otimes\boldsymbol{r}%\right]:\mathbb{E}^{\text{e}}:\left[\boldsymbol{r}\otimes\boldsymbol{r}\right]italic_E start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT = [ bold_italic_r ⊗ bold_italic_r ] : blackboard_E start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : [ bold_italic_r ⊗ bold_italic_r ](44)
C𝒓=[𝒓𝒓]:e:[𝒓𝒓]:subscript𝐶𝒓delimited-[]tensor-product𝒓𝒓superscripte:delimited-[]tensor-product𝒓𝒓\displaystyle C_{\boldsymbol{r}}=\left[\boldsymbol{r}\otimes\boldsymbol{r}%\right]:\mathbb{C}^{\text{e}}:\left[\boldsymbol{r}\otimes\boldsymbol{r}\right]italic_C start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT = [ bold_italic_r ⊗ bold_italic_r ] : blackboard_C start_POSTSUPERSCRIPT e end_POSTSUPERSCRIPT : [ bold_italic_r ⊗ bold_italic_r ](45)

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

ξ𝔼subscript𝜉𝔼\displaystyle\xi_{\mathbb{E}}italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT=min𝒓E𝒓(𝒓)E0.absentsubscript𝒓subscript𝐸𝒓𝒓subscript𝐸0\displaystyle=\frac{\min_{\boldsymbol{r}}E_{\boldsymbol{r}}\!\left({%\boldsymbol{r}}\right)}{E_{0}}.= divide start_ARG roman_min start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT italic_E start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT ( bold_italic_r ) end_ARG start_ARG italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG .(46)

Here, E0subscript𝐸0E_{0}italic_E start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT corresponds to the initial state at t=0𝑡0t=0italic_t = 0, which is isotropic and independent of 𝒓𝒓\boldsymbol{r}bold_italic_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𝒓subscript𝐶𝒓C_{\boldsymbol{r}}italic_C start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT can be compared to the reference resulting in

ζ(𝒓)=C0C𝒓.𝜁𝒓subscript𝐶0subscript𝐶𝒓\displaystyle\zeta(\boldsymbol{r})=\dfrac{C_{0}}{C_{\boldsymbol{r}}}.italic_ζ ( bold_italic_r ) = divide start_ARG italic_C start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_C start_POSTSUBSCRIPT bold_italic_r end_POSTSUBSCRIPT end_ARG .(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

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 [𝒆r,𝒆θ,𝒆z]subscript𝒆rsubscript𝒆𝜃subscript𝒆𝑧\left[\boldsymbol{e}_{\mathrm{r}},\boldsymbol{e}_{\theta},\boldsymbol{e}_{z}\right][ bold_italic_e start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ], cf.Fig.2.

\psfrag{uz}[c][c]{$u$}\psfrag{gz}[c][c]{$\gamma$}\psfrag{g}[c][c]{$\theta$}\psfrag{r}[c][c]{$R$}\psfrag{l}[c][c]{$L$}\psfrag{e1}[c][c]{$\boldsymbol{e}_{r}$}\psfrag{e2}[c][c]{$\boldsymbol{e}_{\theta}$}\psfrag{e3}[c][c]{$\boldsymbol{e}_{z}$}\psfrag{e12}[c][c]{$\theta,\!z$}\includegraphics[width=325.215pt]{Chapter/0X_NumericalResults/TensionShear/Figures/TensionTorsion.eps}

𝜺𝜺\displaystyle\boldsymbol{\varepsilon}bold_italic_ε=[εrrεrθεrzsymεθθεθzsymsymεzz][𝒆r,𝒆θ,𝒆z]absentsubscriptmatrixsubscript𝜀𝑟𝑟missing-subexpressionsubscript𝜀𝑟𝜃missing-subexpressionsubscript𝜀𝑟𝑧symmissing-subexpressionsubscript𝜀𝜃𝜃missing-subexpressionsubscript𝜀𝜃𝑧symmissing-subexpressionsymmissing-subexpressionsubscript𝜀𝑧𝑧subscript𝒆rsubscript𝒆𝜃subscript𝒆𝑧\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]}= [ start_ARG start_ROW start_CELL italic_ε start_POSTSUBSCRIPT italic_r italic_r end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT italic_r italic_θ end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT italic_r italic_z end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT italic_θ italic_θ end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT
𝝈𝝈\displaystyle\boldsymbol{\sigma}bold_italic_σ=[σrrσrθσrzsymσθθσθzsymsymσzz][𝒆r,𝒆θ,𝒆z]absentsubscriptmatrixsubscript𝜎𝑟𝑟missing-subexpressionsubscript𝜎𝑟𝜃missing-subexpressionsubscript𝜎𝑟𝑧symmissing-subexpressionsubscript𝜎𝜃𝜃missing-subexpressionsubscript𝜎𝜃𝑧symmissing-subexpressionsymmissing-subexpressionsubscript𝜎𝑧𝑧subscript𝒆rsubscript𝒆𝜃subscript𝒆𝑧\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]}= [ start_ARG start_ROW start_CELL italic_σ start_POSTSUBSCRIPT italic_r italic_r end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT italic_r italic_θ end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT italic_r italic_z end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT italic_θ italic_θ end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_θ end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT italic_z end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT

A virtual tension-torsion sequence serves as the reference load path. Within this test, torsional strain εθzsubscript𝜀𝜃𝑧\varepsilon_{\theta z}italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT is increased first up to εθz=0.03subscript𝜀𝜃𝑧0.03\varepsilon_{\theta z}=0.03italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT = 0.03 (while εzz=0=const.subscript𝜀𝑧𝑧0const.\varepsilon_{zz}=0=\text{const.}italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT = 0 = const.), and subsequently tensile strain εzzsubscript𝜀𝑧𝑧\varepsilon_{zz}italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT is increased up to εzz=0.06subscript𝜀𝑧𝑧0.06\varepsilon_{zz}=0.06italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT = 0.06 (while εθz=0.03=const.subscript𝜀𝜃𝑧0.03const.\varepsilon_{\theta z}=0.03=\text{const.}italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT = 0.03 = const.). The final strain components at normalized time t=1𝑡1t=1italic_t = 1 are εzz=0.06subscript𝜀𝑧𝑧0.06\varepsilon_{zz}=0.06italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT = 0.06 and εθz=0.03subscript𝜀𝜃𝑧0.03\varepsilon_{\theta z}=0.03italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT = 0.03. The evolution of the stress-triaxiality corresponding to this reference load path is denoted by ηref(t)superscript𝜂ref𝑡\eta^{\text{ref}}(t)italic_η start_POSTSUPERSCRIPT ref end_POSTSUPERSCRIPT ( italic_t ) and the Frobenius-norm of the final elastic stiffness tensor is f𝔼ref(t=1)superscriptsubscript𝑓𝔼ref𝑡1f_{\mathbb{E}}^{\text{ref}}(t=1)italic_f start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ref end_POSTSUPERSCRIPT ( italic_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 εθzsubscript𝜀𝜃𝑧\varepsilon_{\theta z}italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT and εzzsubscript𝜀𝑧𝑧\varepsilon_{zz}italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT are discretized by means of fifth-order Bézier curves, see Eq.(39). While the initial conditions are εθz(t=0)=εzz(t=0)=0subscript𝜀𝜃𝑧𝑡0subscript𝜀𝑧𝑧𝑡00\varepsilon_{\theta z}(t=0)=\varepsilon_{zz}(t=0)=0italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT ( italic_t = 0 ) = italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT ( italic_t = 0 ) = 0, the final state is defined by εθz(t=1)=0.03subscript𝜀𝜃𝑧𝑡10.03\varepsilon_{\theta z}(t=1)=0.03italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT ( italic_t = 1 ) = 0.03 and εzz(t=1)=0.06subscript𝜀𝑧𝑧𝑡10.06\varepsilon_{zz}(t=1)=0.06italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT ( italic_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:

maxQif𝔼(t=1)s.t.ηopt(t)ηref(t)t[0,1]formulae-sequencesubscript𝑄𝑖maxsubscript𝑓𝔼𝑡1𝑠𝑡superscript𝜂opt𝑡superscript𝜂ref𝑡for-all𝑡01\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]start_UNDERACCENT italic_Q start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_UNDERACCENT start_ARG max end_ARG italic_f start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT ( italic_t = 1 ) italic_s . italic_t . italic_η start_POSTSUPERSCRIPT roman_opt end_POSTSUPERSCRIPT ( italic_t ) ≥ italic_η start_POSTSUPERSCRIPT roman_ref end_POSTSUPERSCRIPT ( italic_t ) ∀ italic_t ∈ [ 0 , 1 ](48)

Accordingly, we search for the load path in terms of control points Qisubscript𝑄𝑖Q_{i}italic_Q start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT that maximizes the average elastic stiffness. This path has to fulfill the constraint ηopt(t)ηref(t)tsuperscript𝜂opt𝑡superscript𝜂ref𝑡for-all𝑡\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)\,\forall\,titalic_η start_POSTSUPERSCRIPT roman_opt end_POSTSUPERSCRIPT ( italic_t ) ≥ italic_η start_POSTSUPERSCRIPT roman_ref end_POSTSUPERSCRIPT ( italic_t ) ∀ italic_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 ηopt(t)ηref(t)tsuperscript𝜂opt𝑡superscript𝜂ref𝑡for-all𝑡\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)\,\forall\,titalic_η start_POSTSUPERSCRIPT roman_opt end_POSTSUPERSCRIPT ( italic_t ) ≥ italic_η start_POSTSUPERSCRIPT roman_ref end_POSTSUPERSCRIPT ( italic_t ) ∀ italic_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).

\begin{overpic}[width=424.94574pt]{Chapter/0X_NumericalResults/TensionShear/%Figures/finxi.eps}\put(25.0,23.0){\rotatebox{90.0}{\scriptsize Reference}}\put(31.0,23.0){\rotatebox{90.0}{\scriptsize Optimized}}\end{overpic}

\psfrag{p}[c][c]{\raisebox{5.69046pt}{\scalebox{0.7}{Equiv. plastic strain $\varepsilon^{\text{p,eq}}$ $\left[-\right]$}}}\psfrag{1Ia}[c][c]{\scalebox{0.6}{\text{ECC-a}}}\psfrag{2Ib}[c][c]{\scalebox{0.6}{\text{ECC-i}}}\psfrag{3IIa}[c][c]{\scalebox{0.6}{\text{LEM-a}}}\psfrag{4IIb}[c][c]{\scalebox{0.6}{\text{LEM-i}}}\includegraphics[width=420.61192pt]{Chapter/0X_NumericalResults/TensionShear/Figures/finp.eps}

\psfrag{stiffness I-a}[c][c]{\scalebox{0.7}{Elastic stiffness ratio $\xi_{\mathbb{E}}\,\left[-\right]$}}\psfrag{p}[c][c]{\raisebox{-17.07182pt}{\scalebox{0.7}{\shortstack[c]{Equivalent plastic strain\\$\varepsilon^{\text{p,eq}}$ $\left[-\right]$}}}}\includegraphics[width=429.28616pt]{Chapter/0X_NumericalResults/TensionShear/Figures/xi.eps}

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 εp,eqsuperscript𝜀p,eq\varepsilon^{\text{p,eq}}italic_ε start_POSTSUPERSCRIPT p,eq end_POSTSUPERSCRIPT 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 ξ𝔼subscript𝜉𝔼\xi_{\mathbb{E}}italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT-axis and the εp,eqsuperscript𝜀p,eq\varepsilon^{\text{p,eq}}italic_ε start_POSTSUPERSCRIPT p,eq end_POSTSUPERSCRIPT-axis, respectively. Although the curves associated with optimal paths are below their respective reference (i.e. less remaining stiffness at distinct values of εp,eqsuperscript𝜀p,eq\varepsilon^{\text{p,eq}}italic_ε start_POSTSUPERSCRIPT p,eq end_POSTSUPERSCRIPT), the final remaining stiffness can be increased by reducing the path length in εp,eqsuperscript𝜀p,eq\varepsilon^{\text{p,eq}}italic_ε start_POSTSUPERSCRIPT p,eq end_POSTSUPERSCRIPT. Therefore, the equivalent plastic strain εp,eqsuperscript𝜀p,eq\varepsilon^{\text{p,eq}}italic_ε start_POSTSUPERSCRIPT p,eq end_POSTSUPERSCRIPT 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 ηopt(t)ηref(t)superscript𝜂opt𝑡superscript𝜂ref𝑡\eta^{\mathrm{opt}}(t)\geq\eta^{\mathrm{ref}}(t)italic_η start_POSTSUPERSCRIPT roman_opt end_POSTSUPERSCRIPT ( italic_t ) ≥ italic_η start_POSTSUPERSCRIPT roman_ref end_POSTSUPERSCRIPT ( italic_t ) is indeed fulfilled at any time.

\psfrag{Strain11}[c][c]{\scalebox{0.8}{Tensile strain $\varepsilon_{zz}$}}\psfrag{Strain12}[c][c]{\scalebox{0.8}{Shear strain $\varepsilon_{\theta z}$}}\includegraphics[width=433.62pt]{Chapter/0X_NumericalResults/TensionShear/Figures/StrainPaths.eps}

\psfrag{Time}[c][c]{\scalebox{0.8}{Normalized time $t\,\left[-\right]$}}\psfrag{triax opt}[c][c]{\scalebox{0.8}{Stress triaxiality $\eta\,\left[-\right]$}}\includegraphics[width=433.62pt]{Chapter/0X_NumericalResults/TensionShear/Figures/Triaxiality.eps}

\psfrag{Time}[c][c]{\scalebox{0.8}{Normalized time $t\,\left[-\right]$}}\psfrag{p}[t][c]{\scalebox{0.8}{Cumul. eq. pl. strain $p\,\left[-\right]$}}\psfrag{lode I-b}[c][c]{\scalebox{0.8}{Lode angle par. $\bar{\theta}\,\left[-\right]$}}\includegraphics[width=433.62pt]{Chapter/0X_NumericalResults/TensionShear/Figures/LoAnPar.eps}

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 εθzsubscript𝜀𝜃𝑧\varepsilon_{\theta z}italic_ε start_POSTSUBSCRIPT italic_θ italic_z end_POSTSUBSCRIPT and axial component εzzsubscript𝜀𝑧𝑧\varepsilon_{zz}italic_ε start_POSTSUBSCRIPT italic_z italic_z end_POSTSUBSCRIPT 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𝑡1t=1italic_t = 1 is given in Fig.6. It shows the directional elastic stiffness ratio 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=52.03227pt]{Chapter/0X_NumericalResults/TensionShear/Figures/CoordinateSystemPSFRAG.eps}\psfrag{E1}[c][c]{$\zeta$}\includegraphics[width=26.01613pt]{Chapter/0X_NumericalResults/TensionShear/Figures/StiffnessEkh_Leg.eps}

\psfrag{A}[c][c]{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}Ref.}}\psfrag{B}[c][c]{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}Opt.}}\includegraphics[width=385.92152pt]{Chapter/0X_NumericalResults/TensionShear/Figures/StiffnessEkh_named.eps}

\psfrag{A}[c][c]{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}Ref.}}\psfrag{B}[c][c]{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}Opt.}}\includegraphics[width=385.92152pt]{Chapter/0X_NumericalResults/TensionShear/Figures/StiffnessLad_named.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 relative directional stiffness ratio of ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 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 ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 0.8 is reached. For the whole load path, the stress triaxiality and the Lode angle parameter are kept constant (proportional loading). 17×13171317\times 1317 × 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

𝝈=[σ1100symσ22(σ11,η,θ¯)0symsymσ33(σ11,η,θ¯)][𝒆1,𝒆2,𝒆3]𝝈subscriptmatrixsubscript𝜎11missing-subexpression0missing-subexpression0symmissing-subexpressionsubscript𝜎22subscript𝜎11𝜂¯𝜃missing-subexpression0symmissing-subexpressionsymmissing-subexpressionsubscript𝜎33subscript𝜎11𝜂¯𝜃subscript𝒆1subscript𝒆2subscript𝒆3\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]}bold_italic_σ = [ start_ARG start_ROW start_CELL italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT end_CELL start_CELL end_CELL start_CELL 0 end_CELL start_CELL end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT ( italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT , italic_η , over¯ start_ARG italic_θ end_ARG ) end_CELL start_CELL end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 33 end_POSTSUBSCRIPT ( italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT , italic_η , over¯ start_ARG italic_θ end_ARG ) end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT(49)

where η𝜂\etaitalic_η is the (prescribed) stress triaxiality and θ¯¯𝜃\bar{\theta}over¯ start_ARG italic_θ end_ARG the (prescribed) Lode angle parameter. Closed-form solutions for σ22(σ11,η,L(θ¯))subscript𝜎22subscript𝜎11𝜂𝐿¯𝜃\sigma_{22}(\sigma_{11},\eta,L\!\left({\bar{\theta}}\right))italic_σ start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT ( italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT , italic_η , italic_L ( over¯ start_ARG italic_θ end_ARG ) ) and σ33(σ11,η,L(θ¯))subscript𝜎33subscript𝜎11𝜂𝐿¯𝜃\sigma_{33}(\sigma_{11},\eta,L\!\left({\bar{\theta}}\right))italic_σ start_POSTSUBSCRIPT 33 end_POSTSUBSCRIPT ( italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT , italic_η , italic_L ( over¯ start_ARG italic_θ end_ARG ) ) are given inA. Thus, the only free loading parameter is σ11subscript𝜎11\sigma_{11}italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT. 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.8. 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=390.25534pt]{Chapter/0X_NumericalResults/FractureSurfaces/Figures/FS_EKH.eps}

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

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

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

\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.8 (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 Y¯¯𝑌\bar{Y}over¯ start_ARG italic_Y end_ARG 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., 𝝈=σ11𝒆1𝒆1𝝈tensor-productsubscript𝜎11subscript𝒆1subscript𝒆1\boldsymbol{\sigma}=\sigma_{11}\,\boldsymbol{e}_{1}\otimes\boldsymbol{e}_{1}bold_italic_σ = italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ⊗ bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. It corresponds to a stress triaxiality of η=1/3𝜂13\eta=1/3italic_η = 1 / 3 and Lode angle parameter of θ¯=1¯𝜃1\bar{\theta}=1over¯ start_ARG italic_θ end_ARG = 1. Inspired by uniaxial tension, the following more general stress path is considered:

𝝈(t)=[σ11(t)00symσ22(σ11(t),η(t))0symsym0][𝒆1,𝒆2,𝒆3]𝝈𝑡subscriptmatrixsubscript𝜎11𝑡missing-subexpression0missing-subexpression0symmissing-subexpressionsubscript𝜎22subscript𝜎11𝑡𝜂𝑡missing-subexpression0symmissing-subexpressionsymmissing-subexpression0subscript𝒆1subscript𝒆2subscript𝒆3\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]}bold_italic_σ ( italic_t ) = [ start_ARG start_ROW start_CELL italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_t ) end_CELL start_CELL end_CELL start_CELL 0 end_CELL start_CELL end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL italic_σ start_POSTSUBSCRIPT 22 end_POSTSUBSCRIPT ( italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_t ) , italic_η ( italic_t ) ) end_CELL start_CELL end_CELL start_CELL 0 end_CELL end_ROW start_ROW start_CELL sym end_CELL start_CELL end_CELL start_CELL sym end_CELL start_CELL end_CELL start_CELL 0 end_CELL end_ROW end_ARG ] start_POSTSUBSCRIPT [ bold_italic_e start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT , bold_italic_e start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] end_POSTSUBSCRIPT(50)

Clearly, by setting η=1/3𝜂13\eta=1/3italic_η = 1 / 3 and θ¯=1¯𝜃1\bar{\theta}=1over¯ start_ARG italic_θ end_ARG = 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., η(t=0)=η(t=1)=1/3𝜂𝑡0𝜂𝑡113\eta(t=0)=\eta(t=1)=1/3italic_η ( italic_t = 0 ) = italic_η ( italic_t = 1 ) = 1 / 3, it reaches its maximum deviation from this state at t=0.5𝑡0.5t=0.5italic_t = 0.5 with a stress triaxiality of η(t=0.5)=1/37/30=1/10𝜂𝑡0.513730110\eta(t=0.5)=1/3-7/30=1/10italic_η ( italic_t = 0.5 ) = 1 / 3 - 7 / 30 = 1 / 10 for the first load path and η(t=0.5)=1/3+7/30=17/30𝜂𝑡0.5137301730\eta(t=0.5)=1/3+7/30=17/30italic_η ( italic_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 σ11subscript𝜎11\sigma_{11}italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT and η𝜂\etaitalic_η. Finally, it is noted that the time-scaling of σ11(t)subscript𝜎11𝑡\sigma_{11}(t)italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_t ) is iteratively adapted until the state σ11(t=1)subscript𝜎11𝑡1\sigma_{11}(t=1)italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_t = 1 ) and η(t=1)𝜂𝑡1\eta(t=1)italic_η ( italic_t = 1 ) correspond to uniaxial tension and a damage threshold of ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 0.8.

The aforementioned two different load paths are shown in Fig.9 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.9, these paths are non-proportional and even more important, they lead to different equivalent plastic strains at the final time t=1𝑡1t=1italic_t = 1 corresponding to ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 0.8, see the zoom-in in Fig.9 (right). The deviation of the equivalent plastic strain for these paths is above 30%percent3030\%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.10.

Limits of isotropic damage models for complex load paths — beyond stress triaxiality and Lode angle parameter (1)
\begin{overpic}[width=216.81pt]{Chapter/0X_NumericalResults/FractureSurfaces/%Figures/barP2FS.eps}\put(23.5,24.0){\rotatebox{90.0}{\scriptsize Proportional}}\put(28.0,28.0){\rotatebox{90.0}{\scriptsize Non-prop. 1}}\put(32.5,22.0){\rotatebox{90.0}{\scriptsize Non-prop. 2}}\end{overpic}

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 σ11(t)subscript𝜎11𝑡\sigma_{11}(t)italic_σ start_POSTSUBSCRIPT 11 end_POSTSUBSCRIPT ( italic_t ) is iteratively adapted to synchronize η(t=1)𝜂𝑡1\eta(t=1)italic_η ( italic_t = 1 ) with the predefined equivalent plastic strain ε¯p,eq=0.0462superscript¯𝜀peq0.0462\bar{\varepsilon}^{\mathrm{p,eq}}=0.0462over¯ start_ARG italic_ε end_ARG start_POSTSUPERSCRIPT roman_p , roman_eq end_POSTSUPERSCRIPT = 0.0462. This specific value corresponds to the respective location on the damage iso-surface (uniaxial tension, proportional load path, ECC-a, η=1/3𝜂13\eta=1/3italic_η = 1 / 3, θ¯=1¯𝜃1\bar{\theta}=1over¯ start_ARG italic_θ end_ARG = 1).

The results in Fig.11 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 (ξ𝔼=1subscript𝜉𝔼1\xi_{\mathbb{E}}=1italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 1) as a reference point, lower stress triaxialities increase the remaining stiffness ξ𝔼subscript𝜉𝔼\xi_{\mathbb{E}}italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT by more than 9%percent99\,\%9 %, while the higher stress triaxialities decrease it by 4%percent44\,\%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.

Limits of isotropic damage models for complex load paths — beyond stress triaxiality and Lode angle parameter (2)

\psfrag{D}[c][c]{\scalebox{0.8}{Elastic stiffness ratio $\xi_{\mathbb{E}}$ [-]}}\psfrag{EP}[c][c]{\scalebox{0.8}{Equivalent plastic strain $\varepsilon^{\text{p,eq}}$ [-]}}\includegraphics[width=433.62pt]{Chapter/0X_NumericalResults/FractureSurfaces/Figures/P2FS3_EKH_D_max.eps}

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 η𝜂\etaitalic_η and θ¯¯𝜃\bar{\theta}over¯ start_ARG italic_θ end_ARG). This is precisely the idea analyzed in this paragraph. For that purpose, two different stress tensors are considered. While 𝝈𝝈\boldsymbol{\sigma}bold_italic_σ corresponds to proportional loading, 𝝈~~𝝈\tilde{\boldsymbol{\sigma}}over~ start_ARG bold_italic_σ end_ARG is associated with a more complex load path. Both share the same eigenvalues σi(t)=σ~i(t)subscript𝜎𝑖𝑡subscript~𝜎𝑖𝑡\sigma_{i}\!\left({t}\right)=\tilde{\sigma}_{i}\!\left({t}\right)italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_t ) = over~ start_ARG italic_σ end_ARG start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( italic_t ) (and thus the same triaxiality and Lode parameter). But they have different directions in terms of the associated eigenvectors. 𝝈𝝈\boldsymbol{\sigma}bold_italic_σ and 𝝈~~𝝈\tilde{\boldsymbol{\sigma}}over~ start_ARG bold_italic_σ end_ARG can be implemented by means of the spectral decomposition, i.e.,

𝝈𝝈\displaystyle\boldsymbol{\sigma}bold_italic_σ=i=13σi𝑵i𝑵iabsentsuperscriptsubscript𝑖13tensor-productsubscript𝜎𝑖superscript𝑵𝑖superscript𝑵𝑖\displaystyle=\sum_{i=1}^{3}\sigma_{i}\,\boldsymbol{N}^{i}\otimes\boldsymbol{N%}^{i}= ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT bold_italic_N start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ⊗ bold_italic_N start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT(51)
𝝈~~𝝈\displaystyle\tilde{\boldsymbol{\sigma}}over~ start_ARG bold_italic_σ end_ARG=i=13σi[𝑹𝑵i][𝑹𝑵i]with𝑹(α,β,γ)=𝑹𝒆¯¯1(α)𝑹𝒆¯2(β)𝑹𝒆3(γ).formulae-sequenceabsentsuperscriptsubscript𝑖13tensor-productsubscript𝜎𝑖delimited-[]𝑹superscript𝑵𝑖delimited-[]𝑹superscript𝑵𝑖with𝑹𝛼𝛽𝛾subscript𝑹subscript¯¯𝒆1𝛼subscript𝑹subscript¯𝒆2𝛽subscript𝑹subscript𝒆3𝛾.\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{.}= ∑ start_POSTSUBSCRIPT italic_i = 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT [ bold_italic_R ⋅ bold_italic_N start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ] ⊗ [ bold_italic_R ⋅ bold_italic_N start_POSTSUPERSCRIPT italic_i end_POSTSUPERSCRIPT ] with bold_italic_R ( italic_α , italic_β , italic_γ ) = bold_italic_R start_POSTSUBSCRIPT over¯ start_ARG over¯ start_ARG bold_italic_e end_ARG end_ARG start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_α ) ⋅ bold_italic_R start_POSTSUBSCRIPT over¯ start_ARG bold_italic_e end_ARG start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_β ) ⋅ bold_italic_R start_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT end_POSTSUBSCRIPT ( italic_γ ) .(52)

where 𝑵isubscript𝑵𝑖\boldsymbol{N}_{i}bold_italic_N start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT are the eigenvectors of 𝝈𝝈\boldsymbol{\sigma}bold_italic_σ and 𝑹𝒆isubscript𝑹subscript𝒆𝑖\boldsymbol{R}_{\boldsymbol{e}_{i}}bold_italic_R start_POSTSUBSCRIPT bold_italic_e start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUBSCRIPT are rotation matrices depending on the Euler angles α𝛼\alphaitalic_α, β𝛽\betaitalic_β and γ𝛾\gammaitalic_γ. Four different mechanical tests are considered:

  1. 1.

    Uniaxial tension: η=1/3𝜂13\eta=1/3italic_η = 1 / 3 and θ¯=1¯𝜃1\bar{\theta}=1over¯ start_ARG italic_θ end_ARG = 1

  2. 2.

    Combined tension and shear: η=1/6𝜂16\eta=1/6italic_η = 1 / 6 and θ¯=0.506¯𝜃0.506\bar{\theta}=0.506over¯ start_ARG italic_θ end_ARG = 0.506

  3. 3.

    Simple shear: η=0𝜂0\eta=0italic_η = 0 and θ¯=0¯𝜃0\bar{\theta}=0over¯ start_ARG italic_θ end_ARG = 0

  4. 4.

    Combined compression and shear: η=1/6𝜂16\eta=-1/6italic_η = - 1 / 6 and θ¯=0.506¯𝜃0.506\bar{\theta}=-0.506over¯ start_ARG italic_θ end_ARG = - 0.506

While the first eigenvalue σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT is increased, the second as well as the third follow from the prescribed and constant parameters η𝜂\etaitalic_η and θ¯¯𝜃\bar{\theta}over¯ start_ARG italic_θ end_ARG, seeA for details. For the proportional load paths, stress tensor(51) is considered with constant eigenvectors 𝑵isubscript𝑵𝑖\boldsymbol{N}_{i}bold_italic_N start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT is increased up to a damage threshold of ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 0.8. In the case of the non-proportional load paths, the Euler angles α𝛼\alphaitalic_α, β𝛽\betaitalic_β and γ𝛾\gammaitalic_γ evolve as illustrated in Fig.12.In line with the example studied in the previous paragraph, the temporal discretization of σ1(t)subscript𝜎1𝑡\sigma_{1}(t)italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ( italic_t ) is chosen such that all Euler angles vanish at final time t=1𝑡1t=1italic_t = 1 when damage threshold ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 0.8 is reached, i.e., the time scaling of σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT and the Euler angles is synchronized.

The bar charts in Fig.13 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 ξ𝔼=0.8subscript𝜉𝔼0.8\xi_{\mathbb{E}}=0.8italic_ξ start_POSTSUBSCRIPT blackboard_E end_POSTSUBSCRIPT = 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{time}[c][c]{\scalebox{0.8}{Time $t\,[-]$}}\psfrag{angle}[c][c]{\scalebox{0.8}{Angle $[\text{Rad}]$}}\psfrag{-0.4pi}[c][c]{\scalebox{0.8}{$-0.4\,\pi$}}\psfrag{-0.2pi}[c][c]{\scalebox{0.8}{$-0.2\,\pi$}}\psfrag{0pi}[c][c]{\scalebox{0.8}{$0$}}\psfrag{0.2pi}[c][c]{\scalebox{0.8}{$0.2\,\pi$}}\psfrag{0.4pi}[c][c]{\scalebox{0.8}{$0.4\,\pi$}}\psfrag{alp}[l][l]{\scalebox{0.8}{$\alpha$}}\psfrag{bet}[l][l]{\scalebox{0.8}{$\beta$}}\psfrag{gam}[l][l]{\scalebox{0.8}{$\gamma$}}\psfrag{0}[c][c]{\scalebox{0.8}{$0$}}\psfrag{0.5}[c][c]{\scalebox{0.8}{$0.5$}}\psfrag{1}[c][c]{\scalebox{0.8}{$1$}}\includegraphics[width=303.53267pt]{Chapter/0X_NumericalResults/IdenticalParameters/Figures/UniaxialTension/Angles.eps}

Limits of isotropic damage models for complex load paths — beyond stress triaxiality and Lode angle parameter (3)\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}

\begin{overpic}[width=216.81pt]{Chapter/0X_NumericalResults/FractureSurfaces/%Figures/barEAC.eps}\put(58.5,10.0){\rotatebox{90.0}{\scriptsize Tension}}\put(61.8,11.0){\rotatebox{90.0}{\scriptsize Tension + Shear}}\put(65.2,11.0){\rotatebox{90.0}{\scriptsize Shear}}\put(69.0,13.0){\rotatebox{90.0}{\scriptsize Compression + Shear}}\end{overpic}

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 Parameter realization

Stress triaxiality(1)and Lode parameter(2) can be reformulated in terms of eigenvalues

[σ2]1,2subscriptdelimited-[]subscript𝜎212\displaystyle\left[{\sigma_{2}}\right]_{1,2}[ italic_σ start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT=χσ1±χ2ω|σ1|absentplus-or-minus𝜒subscript𝜎1superscript𝜒2𝜔subscript𝜎1\displaystyle=-\chi\,\sigma_{1}\pm\sqrt{\chi^{2}-\omega}\left|\sigma_{1}\right|= - italic_χ italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ± square-root start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ω end_ARG | italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT |(53)
[σ3]1,2subscriptdelimited-[]subscript𝜎312\displaystyle\left[{\sigma_{3}}\right]_{1,2}[ italic_σ start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT ] start_POSTSUBSCRIPT 1 , 2 end_POSTSUBSCRIPT=[a+bχ]σ1bχ2ω|σ1|absentminus-or-plusdelimited-[]𝑎𝑏𝜒subscript𝜎1𝑏superscript𝜒2𝜔subscript𝜎1\displaystyle=\left[{a+b\,\chi}\right]\,\sigma_{1}\mp b\,\sqrt{\chi^{2}-\omega%}\,\left|\sigma_{1}\right|= [ italic_a + italic_b italic_χ ] italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT ∓ italic_b square-root start_ARG italic_χ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_ω end_ARG | italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT |(54)

with

a𝑎\displaystyle aitalic_a=L+1L1,b=2L1formulae-sequenceabsent𝐿1𝐿1𝑏2𝐿1\displaystyle=\frac{L+1}{L-1},\quad b=\frac{2}{L-1}= divide start_ARG italic_L + 1 end_ARG start_ARG italic_L - 1 end_ARG , italic_b = divide start_ARG 2 end_ARG start_ARG italic_L - 1 end_ARG

and

χ𝜒\displaystyle\chiitalic_χ=18η2[LL22]4L[L3]18η2[L23]2[L3]2,ω=18η2[L2L+2]8L218η2[L23]2[L3]2.formulae-sequenceabsent18superscript𝜂2delimited-[]𝐿superscript𝐿224𝐿delimited-[]𝐿318superscript𝜂2delimited-[]superscript𝐿232superscriptdelimited-[]𝐿32𝜔18superscript𝜂2delimited-[]superscript𝐿2𝐿28superscript𝐿218superscript𝜂2delimited-[]superscript𝐿232superscriptdelimited-[]𝐿32.\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{.}= divide start_ARG 18 italic_η start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ italic_L - italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 2 ] - 4 italic_L [ italic_L - 3 ] end_ARG start_ARG 18 italic_η start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 3 ] - 2 [ italic_L - 3 ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG , italic_ω = divide start_ARG 18 italic_η start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - italic_L + 2 ] - 8 italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG 18 italic_η start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT [ italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 3 ] - 2 [ italic_L - 3 ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG .

They allow to enforce both the stress triaxiality as well as the Lode parameter (and hence the Lode angle parameter) and to apply loading by increasing σ1subscript𝜎1\sigma_{1}italic_σ start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT. Clearly, special care is required, if L1𝐿1L\rightarrow 1italic_L → 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.
Limits of isotropic damage models for complex load paths — beyond stress triaxiality and Lode angle parameter (2024)

FAQs

What is lode angle and triaxiality? ›

The stress triaxiality and the Lode angle parameter are two well established stress invariants for the charac- terization of damage evolution. This work assesses the limits of this tuple by using it for damage predictions via ductile damage models from continuum damage mechanics.

What is stress triaxiality in Abaqus? ›

Stress Triaxiality is defined as the ratio of (HYDROSTATIC PRESSURE OR MEAN STRESS) to the EQUIVALENT STRESS. When I check the abaqus field output variables this option is available but only in the recent versions (6.8 & 6.9). The input file usage according to the documentation is TRIAX.

What is the lode angle parameter? ›

The Lode parameter is a function of the third invariant of the stress deviator and is used to distinguish between the different shear stress states in three dimensions (3-D), ranging from axisymmetric tension to biaxial tension with axisymmetric compression and passing through in-plane shear.

What is the normalized lode angle? ›

Normalized Lode angle θ is called as Lode parameter. From its parameter, it is possible to know whether the current stress state is tension or pure shear or compression.

How to measure stress triaxiality? ›

The formula used to define the stress triaxiality is the ratio of hydrostatic (mean) stress (σm) to von Mises (equivalent) stress (σe), i.e., σm/ σe 1.

Why is stress triaxiality important? ›

Stress triaxiality has important applications in fracture mechanics and can often be used to predict the type of fracture (i.e. ductile or brittle) within the region defined by that stress state.

What is the damage model in ABAQUS? ›

ABAQUS offers a general capability for modeling progressive damage and failure in engineering structures. – Material failure refers to the complete loss of load carrying capacity that results from progressive degradation of the material stiffness. – Stiffness degradation is modeled using damage mechanics.

What is the lode angle in soil mechanics? ›

Lode angle parameter helps understanding behaviour of material in tension, shear and compression when evaluating stress triaxiallity at different strain rates.

What is lode in geology? ›

In geology, a lode is a deposit of metalliferous ore that fills or is embedded in a fissure (or crack) in a rock formation. Lode deposits are distinguished primarily from placer deposits, where the ore has been eroded out from its original depositional environment and redeposited by sedimentary forces.

What is the difference between torque angle and load angle? ›

But the fact is that Load angle or power angle is the angle between terminal voltage and excitation voltage. Torque angle is the angle between excitation field mmf and resultant air gap mmf…. in a synchronous generator. The load angle and torque angles differ due to the presence of armature leakage emf.

What is the stress triaxiality for uniaxial tension? ›

The theoretical values of stress triaxiality and Lode parameter for uniaxial tension are 1/3 and 1, respectively. In simple and pure shear loading, both values are zero.

Top Articles
Latest Posts
Article information

Author: Tuan Roob DDS

Last Updated:

Views: 5960

Rating: 4.1 / 5 (62 voted)

Reviews: 85% of readers found this page helpful

Author information

Name: Tuan Roob DDS

Birthday: 1999-11-20

Address: Suite 592 642 Pfannerstill Island, South Keila, LA 74970-3076

Phone: +9617721773649

Job: Marketing Producer

Hobby: Skydiving, Flag Football, Knitting, Running, Lego building, Hunting, Juggling

Introduction: My name is Tuan Roob DDS, I am a friendly, good, energetic, faithful, fantastic, gentle, enchanting person who loves writing and wants to share my knowledge and understanding with you.