An anisotropic,brittle damage model for finite strains with a generic damage tensor regularization

Abstract

This paper establishes a generic framework for the nonlocal modeling of anisotropic damage at finite strains. By the combination of two recent works, the new framework allows for the flexible incorporation of different established hyperelastic finite strain material formulations into anisotropic damage whilst ensuring mesh-independent results by employing a generic set of micromorphic gradient-extensions. First, the anisotropic damage model, generally satisfying the damage growth criterion, is investigated for the specific choice of a neo-Hookean material on a single element. Next, the model is applied with different gradient-extensions in structural simulations of an asymmetrically notched specimen to identify an efficient choice in the form of a volumetric–deviatoric regularization. Thereafter, the generic framework, which is without loss of generality here specified for a neo-Hookean material with a volumetric–deviatoric gradient-extension, successfully serves for the complex simulation of a pressure-loaded rotor blade. The codes of the material subroutines are accessible to the public at https://doi.org/10.5281/zenodo.11171630.

Introduction

Motivation. The prediction of structural failure requires the precise local modeling of damage and the accurate nonlocal regularization of the softening phenomenon in structural simulations. Both aspects, the modeling and regularization of softening due to damage, constitute fundamental fields of research in continuum mechanics that continuously generate enhanced local and nonlocal solutions. This motivates the proposition of a general modeling framework that combines a flexible local anisotropic damage model with different nonlocal gradient-extensions in this work. The following paragraphs only provide an overview of recent contributions to both fields. For a general overview of damage modeling, we refer to the works (e.g. Besson, 2010; Lemaitre and Desmorat, 2005; Murakami, 2012; Voyiadjis and Kattan, 2009).

Modeling of softening. Mattiello and Desmorat (2021) investigate different evolution laws for symmetric second-order damage tensors with respect to the Lode angle dependency. Dorn and Wulfinghoff (2021) assume a multiplicative decomposition of the deformation gradient into normal and shear crack contributions that directly yields a tension–compression asymmetry. For low cycle fatigue, Ferreira et al. (2022) propose a damage evolution law depending on the stress triaxiality and test the model in material point studies for multiaxial and nonproportional loading paths. In Reese et al. (2021), a local formulation for anisotropic damage is introduced that can incorporate arbitrary hyperelastic energies and generally fulfills the damage growth criterion (Wulfinghoff et al., 2017). Modeling of anisotropic damage based on the compliance tensor is investigated in Görthofer et al. (2022). Petrini et al. (2023) use a dynamic phase-field model with a fourth-order degradation tensor to account for damage anisotropy at infinitesimal strains. Based on the work by Basaran (2023); Voyiadjis and Kattan (2017, 2024) present an unsymmetrical decomposition of the tensorial damage variable to separately account for crack and void induced damage. To capture stiffness recovery after damage evolution, Shojaei and Voyiadjis (2023) lay the theoretical groundwork for statistical continuum damage healing mechanics by considering the sealing and healing effects. Further models concerned with the modeling of damage healing are presented (e.g. Han et al., 2023, 2024). Loiseau et al. (2023) employ a data-driven approach to identify an accurate tensorial representation of anisotropic damage to account for micro-cracking based on a virtual set of beam–particle simulations. Moreover, novel acoustic testing methodologies for damage evolution are proposed in Yu et al. (2024). The multiphysical extension of damage modeling is presented in Zhou et al. (2022) in the context of chemo-mechanical damage.

Regularization of softening. Jirásek and Desmorat (2019) investigate the regularization performance of a new nonlocal integral approach with a damage-dependent interaction distance for pure damage models, damage-plasticity models, and damage with inelastic strain models. Ahmed et al. (2021) present a damage model for concrete that captures tensile, compressive, and shear damage with separate scalar damage variables, where each quantity is regularized by a gradient-extension with three individual length scales. In Zhang et al. (2022), a convenient approach for the implementation of gradient-extended damage into the finite element software Abaqus is presented and validated in complex structural simulations, where in-built Abaqus features, such as contact or element deletion, are additionally exploited. In Holthusen et al. (2020, 2022a), a general notation for micromorphic gradient-extensions is formulated based on a micromorphic tuple, which is investigated for anisotropic damage with different gradient-extensions in van der Velden et al. (2024). Sprave and Menzel (2023) formulate a ductile anisotropic damage model at finite strains and analyze the isolated and combined regularization of damage and plasticity. The regularization of ductile damage in the logarithmic strain space is studied in Friedlein et al. (2023) with a gradient-extension of the plastic hardening variable or the local damage variable. For damage in semicrystalline polymers, Satouri et al. (2022) examine the regularization effects of the damage variable and the hardening state variable. A regularization of damage through peridynamic modeling is investigated in Wu et al. (2022).

Current and future works. In this work, we combine the elastic energy of the anisotropic damage model of Reese et al. (2021) with the micromorphic gradient-extensions investigated in van der Velden et al. (2024) to obtain a generic framework for regularized anisotropic damage at finite strains. The formulation enables the incorporation of arbitrary hyperelastic finite strain energies and, thereby, yields a flexible local damage model. The generic gradient-extensions provide additional nonlocal flexibility in terms of the regularized quantities and the number of nonlocal degrees of freedom. In future works, the iCANN framework of Holthusen et al. (2024) can provide a functional basis for the elastic energy as well as the gradient-extension and can identify the optimal set of parameters.

Outline of the work. The constitutive framework is presented in Section “Constitutive modeling.” The numerical examples are provided in Section “Numerical examples” with single-element studies in Section “Single-element studies,” the study of an asymmetrically notched specimen in Section “Asymmetrically notched specimen,” and the study of a rotor blade specimen in Section “Rotor blade.” The conclusions are drawn in Section “Conclusion.”

Notational conventions. In this work, italic characters a, A denote scalars (i.e. zeroth-order tensors) and bold-face italic characters $\mathit{b}$ , $\mathit{B}$ refer to first- and second-order tensors. The operators $\mathrm{Div}(\bullet )$ and $\mathrm{Grad}(\bullet )$ denote the divergence and gradient operation of a quantity with respect to the reference configuration. The operator $(\bullet {)}^{\text{T}}$ denotes the transpose of a second-order tensor. A $\cdot$ defines the single contraction and a : the double contraction of two tensors. The time derivative of a quantity is given by $\dot{(\bullet )}$ and a fixed variable value by $(\bullet {)}^{\mathrm{\prime }}$ .

Constitutive modeling

Balance equations

The constitutive framework is based on the micromorphic approach of Forest (2009, 2016) and is formulated in the reference configuration. It comprises the balance of linear momentum

$$\mathrm{Div}\left(\mathit{F}\mathit{S}\right)+{\mathit{f}}_0=\mathbf{0}\phantom{{\mathit{t}}_0{\mathit{u}}^{\mathrm{\prime }}}\text{in }{\mathrm{\Omega }}_0$$

(1)

$$\mathit{F}\mathit{S}\cdot {\mathit{n}}_0={\mathit{t}}_0\phantom{\mathbf{0}{\mathit{u}}^{\mathrm{\prime }}}\text{on }{\mathrm{\Gamma }}_{t0}$$

(2)

$$\mathit{u}={\mathit{u}}^{\mathrm{\prime }}\phantom{\mathbf{0}{\mathit{t}}_0}\text{on }{\mathrm{\Gamma }}_{u0}$$

(3)

and the micromorphic balance equation

$$\mathrm{Div}({\mathbf{\Xi }}_{0_i}-{\mathbf{\Xi }}_{0_e})-{\mathit{\xi }}_{0_i}+{\mathit{\xi }}_{0_e}=\mathbf{0}\phantom{{\mathit{\xi }}_{0_c}{\overline{\mathbf{d}}}^{\mathrm{\prime }}}\text{in }{\mathrm{\Omega }}_0$$

(4)

$$({\mathbf{\Xi }}_{0_i}-{\mathbf{\Xi }}_{0_e})\cdot {\mathit{n}}_0={\mathit{\xi }}_{0_c}\phantom{\mathbf{0}{\overline{\mathbf{d}}}^{\mathrm{\prime }}}\text{on }{\mathrm{\Gamma }}_{c0}$$

(5)

$$\overline{\mathbf{d}}={\overline{\mathbf{d}}}^{\mathrm{\prime }}\phantom{\mathbf{0}{\mathit{\xi }}_{0_c}}\text{o}n{\mathrm{\Gamma }}_{\overline{d}0}$$

(6)

with their corresponding Neumann and Dirichlet boundary conditions. In equations (1) to (6), $\mathit{F}$ denotes the deformation gradient, $\mathit{S}$ the second Piola–Kirchhoff stress tensor, ${\mathit{f}}_0$ the volume forces, ${\mathit{n}}_0$ the outward normal vector, ${\mathit{t}}_0$ the prescribed tractions on the mechanical Neumann boundary ${\mathrm{\Gamma }}_{t0}$ , $\mathit{u}$ the mechanical displacements, and ${\mathit{u}}^{\mathrm{\prime }}$ their fixed values on the mechanical Dirichlet boundary ${\mathrm{\Gamma }}_{u0}$ . Furthermore, ${\mathbf{\Xi }}_{0_i}$ and ${\mathit{\xi }}_{0_i}$ denote the micromorphic internal forces, ${\mathbf{\Xi }}_{0_e}$ and ${\mathit{\xi }}_{0_e}$ the micromorphic external forces, ${\mathit{\xi }}_{0_c}$ the micromorphic contact forces on the micromorphic Neumann boundary ${\mathrm{\Gamma }}_{c0}$ , $\overline{\mathbf{d}}$ the tuple of the nonlocal field variables, and ${\overline{\mathbf{d}}}^{\mathrm{\prime }}$ its fixed value on the micromorphic Dirichlet boundary ${\mathrm{\Gamma }}_{\overline{d}0}$ . Analogously to Brepols et al. (2020); Holthusen et al. (2020, 2022a), the micromorphic external and contact forces are neglected in this work, that is, ${\mathbf{\Xi }}_{0_e}=\mathbf{0}$ , ${\mathit{\xi }}_{0_e}=\mathbf{0}$ and ${\mathit{\xi }}_{0_c}=\mathbf{0}$ , and micromorphic Dirichlet boundary conditions are not considered, that is, ${\mathrm{\Gamma }}_{\overline{d}0}=\mathrm{\varnothing }$ . With these simplifications and using the test functions $\delta \mathit{u}$ and $\delta \overline{\mathbf{d}}$ , the following weak forms with the virtual Green–Lagrange strain $\delta \mathit{E}$ are obtained (cf. Holthusen et al., 2022a)

$$g_u(\mathit{u},\overline{\mathbf{d}},\delta \mathit{u}):={\int }_{{\mathrm{\Omega }}_0}\mathit{S}:\delta \mathit{E}\mathrm{d}V-{\int }_{{\mathrm{\Omega }}_0}{\mathit{f}}_0\cdot \delta \mathit{u}\mathrm{d}V-{\int }_{{\mathrm{\Gamma }}_{t0}}{\mathit{t}}_0\cdot \delta \mathit{u}\mathrm{d}A=0,$$

(7)

$$g_{\overline{d}}(\mathit{u},\overline{\mathbf{d}},\delta \overline{\mathbf{d}}):={\int }_{{\mathrm{\Omega }}_0}{\mathit{\xi }}_{0_i}\cdot \delta \overline{\mathbf{d}}\mathrm{d}V+{\int }_{{\mathrm{\Omega }}_0}{\mathbf{\Xi }}_{0_i}:\mathrm{Grad}\left(\delta \overline{\mathbf{d}}\right)\mathrm{d}V=0$$

(8)

Helmholtz free energy

The Helmholtz free energy consists of four parts that additively compose the total energy

$$\psi (\mathit{C},\mathit{D},{\xi }_d,\mathbf{d},\overline{\mathbf{d}},\mathrm{Grad}\left(\overline{\mathbf{d}}\right))={\psi }_e(\mathit{C},\mathit{D})+{\psi }_d\left({\xi }_d\right)+{\psi }_h\left(\mathit{D}\right)+{\psi }_{\overline{d}}(\mathbf{d},\overline{\mathbf{d}},\mathrm{Grad}\left(\overline{\mathbf{d}}\right))$$

(9)

where ${\psi }_e$ denotes the elastic energy part depending on the right Cauchy–Green tensor $\mathit{C}$ and the second-order damage tensor $\mathit{D}$ , ${\psi }_d$ denotes the isotropic damage hardening energy part depending on the accumulated damage hardening variable ${\xi }_d$ , ${\psi }_h$ denotes the additional kinematic damage hardening energy part depending on $\mathit{D}$ , and ${\psi }_{\overline{d}}$ denotes the micromorphic energy contribution part depending on the local micromorphic tuple d, the nonlocal counterpart $\overline{\mathbf{d}}$ and its gradient $\mathrm{Grad}\left(\overline{\mathbf{d}}\right)$ .

Remark—Novelty

In van der Velden et al. (2024), a framework for the generic regularization of anisotropic damage at finite strains was presented, where multiple invariant-based micromorphic tuples were investigated. However, the former work utilized one specific elastic energy formulated in logarithmic strains that may not preserve the convexity of the model (cf. Neff et al., 2015). Moreover, the damage tensor directly affects only the deviatoric part of the logarithmic strains in a mixed invariant in the aforementioned elastic energy (Holthusen et al., 2022a). Therefore, the research question arises whether the generic regularization of van der Velden et al. (2024) can also be applied to a general class of functions for the elastic energy. Here, we strive to combine the generic regularization of van der Velden et al. (2024) with a generic class of elastic energies accounting for anisotropic damage based on the right Cauchy–Green tensor $\mathit{C}$ . Now, we formulate our elastic energy in line with Reese et al. (2021) according to

$${\psi }_e(\mathit{C},\mathit{D})=f_{\text{d}am}(\mathit{C},\mathit{D}){\psi }_{\star }\left(\mathit{C}\right)$$

(10)

where $f_{\text{d}am}(\mathit{C},\mathit{D})$ denotes a combined isotropic and anisotropic degradation function and ${\psi }_{\star }(\mathit{C})$ denotes the general class of functions depending on the right Cauchy–Green tensor $\mathit{C}$ . This generic formulation of the elastic energy allows the incorporation of arbitrarily established finite strains formulations such as neo-Hookean, Mooney–Rivlin, or Arruda–Boyce models (cf. Steinmann et al., 2012) into the framework of gradient-extended anisotropic damage. So far, the generic formulation of the elastic energy of Reese et al. (2021) was only investigated in material point studies. Hence, we propose a framework with a generic regularization and a generic elastic energy formulation for anisotropic damage at finite strains.

Isothermal Clausius–Duhem inequality

The isothermal Clausius–Duhem inequality including the micromorphic parts (cf. Forest, 2009, 2016) reads

$$-\dot{\psi }+\frac{1}{2}\mathit{S}:\dot{\mathit{C}}+{\mathit{\xi }}_{0_i}\cdot \dot{\overline{\mathbf{d}}}+{\mathbf{\Xi }}_{0_i}:\mathrm{Grad}\left(\dot{\overline{\mathbf{d}}}\right)\ge 0$$

(11)

with the rate of the Helmholtz free energy being

$$\dot{\psi }={\displaystyle \frac{\mathrm{\partial }{\psi }_e}{\mathrm{\partial }\mathit{C}}:\dot{\mathit{C}}+(\underset{=:-{\mathit{Y}}_e}{\underbrace{{\displaystyle \frac{\mathrm{\partial }{\psi }_e}{\mathrm{\partial }\mathit{D}}}}}+\underset{=:{\mathit{Y}}_h}{\underbrace{{\displaystyle \frac{\mathrm{\partial }{\psi }_h}{\mathrm{\partial }\mathit{D}}}}}+\underset{=:{\mathit{Y}}_{\overline{d}}}{\underbrace{{\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\mathit{D}}}}}):\dot{\mathit{D}}+\underset{=:-R_d}{\underbrace{{\displaystyle \frac{\mathrm{\partial }{\psi }_d}{\mathrm{\partial }{\xi }_d}}}}{\dot{\xi }}_d+{\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\overline{\mathbf{d}}}\cdot \dot{\overline{\mathbf{d}}}+{\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\mathrm{Grad}\left(\overline{\mathbf{d}}\right)}:\mathrm{Grad}\left(\dot{\overline{\mathbf{d}}}\right).}}}$$

(12)

In equation (12), the rate of the micromorphic tuple $\dot{\mathbf{d}}$ does not explicitly enter the rate of the Helmholtz free energy $\dot{\psi }$ . Instead, it implicitly contributes through the rate of the damage tensor $\dot{\mathit{D}}$ due to the dependency $\mathbf{d}=\mathbf{d}(\mathit{D})$ in the nonlocal driving force ${\mathit{Y}}_{\overline{d}}$ . The insertion of equation (12) in equation (11) yields

$$\begin{array}{rl}& (\frac{1}{2}\mathit{S}-{\displaystyle \frac{\mathrm{\partial }{\psi }_e}{\mathrm{\partial }\mathit{C}}}):\dot{\mathit{C}}+(\underset{=:\mathit{Y}}{\underbrace{{\mathit{Y}}_e-{\mathit{Y}}_h-{\mathit{Y}}_{\overline{d}}}}):\dot{\mathit{D}}+R_d{\dot{\xi }}_d\\ & +({\mathit{\xi }}_{0_i}-{\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\overline{\mathbf{d}}}})\cdot \dot{\overline{\mathbf{d}}}+({\mathbf{\Xi }}_{0_i}-{\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\mathrm{Grad}\left(\overline{\mathbf{d}}\right)}}):\mathrm{Grad}\left(\dot{\overline{\mathbf{d}}}\right)\ge 0\end{array}$$

(13)

from which the state laws (cf. Brepols et al., 2020; Holthusen et al., 2022a) for the mechanical and generalized micromorphic stresses are obtained as

$$\mathit{S}=2{\displaystyle \frac{\mathrm{\partial }{\psi }_e}{\mathrm{\partial }\mathit{C}},}$$

(14)

$${\mathit{\xi }}_{0_i}={\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\overline{\mathbf{d}}},}$$

(15)

$${\mathbf{\Xi }}_{0_i}={\displaystyle \frac{\mathrm{\partial }{\psi }_{\overline{d}}}{\mathrm{\partial }\mathrm{Grad}\left(\overline{\mathbf{d}}\right)}.}$$

(16)

After the definition of the state laws, the reduced dissipation inequality follows from equation (13) with the thermodynamic driving forces $\mathit{Y}:=-\mathrm{\partial }\psi /\mathrm{\partial }\mathit{D}$ and $R_d$ as

$$\mathit{Y}:\dot{\mathit{D}}+R_d{\dot{\xi }}_d\ge 0.$$

(17)

Damage onset criterion and evolution equations

The damage onset criterion incorporates distortional damage hardening and is formulated analogously to Holthusen et al. (2022a, 2022b) reading

$${\mathrm{\Phi }}_d:=\sqrt{3}\sqrt{{\mathit{Y}}_{+}:\mathbb{A}:{\mathit{Y}}_{+}}-(Y_0-R_d)\le 0$$

(18)

where $\mathbb{A}$ is a fourth-order interaction tensor

$$\mathbb{A}={({(\mathbf{I}-\mathit{D})}^{c_d}\otimes {(\mathbf{I}-\mathit{D})}^{c_d})}^{{}_{\mathrm{T}}^{23}}$$

(19)

and ${\mathit{Y}}_{+}$ the positive semi-definite part of the damage driving force with $⟨\bullet ⟩=\max (\bullet ,0)$

$${\mathit{Y}}_{+}=\sum _{i=1}^3⟨Y_i⟩{\mathit{n}}_i^{\mathit{Y}}\otimes {\mathit{n}}_i^{\mathit{Y}}$$

(20)

where $Y_i$ and ${\mathit{n}}_i^{\mathit{Y}}$ denote the eigenvalues and eigenvectors of $\mathit{Y}$ . The associative evolution equations for the internal variables read

$$\dot{\mathit{D}}={\dot{\gamma }}_d{\displaystyle \frac{\mathrm{\partial }{\mathrm{\Phi }}_d}{\mathrm{\partial }\mathit{Y}},}$$

(21)

$${\dot{\xi }}_d={\dot{\gamma }}_d{\displaystyle \frac{\mathrm{\partial }{\mathrm{\Phi }}_d}{\mathrm{\partial }{R}_d}}$$

(22)

with the Karush-Kuhn–Tucker conditions

$${\dot{\gamma }}_d\ge 0,{\mathrm{\Phi }}_d\le 0,{\dot{\gamma }}_d{\mathrm{\Phi }}_d=0.$$

(23)

Specific Helmholtz free energies and micromorphic tuples

The elastic Helmholtz free energy is formulated in line with Reese et al. (2021), where a hyperelastic energy formulation is multiplied with a linear combination of two degradation functions, and reads

$${\psi }_e=((1-k_{\mathrm{ani}})f_{\mathrm{iso}}+k_{\mathrm{ani}}{f}_{\mathrm{ani}}){\psi }_{\star }$$

(24)

where the material parameter $k_{\mathrm{ani}}\in [0,1]$ defines the degree of damage anisotropy and where the degradation functions $f_{\mathrm{iso}}(\mathit{D})$ and $f_{\mathrm{ani}}(\mathit{C},\mathit{D})$ account for isotropic and anisotropic damage, respectively. This versatile formulation allows for the straightforward incorporation of different established hyperelastic energies ${\psi }_{\star }$ into the framework of anisotropic damage. Here, a neo-Hookean energy is considered with

$${\psi }_{\star }={\psi }_{\mathrm{NH}}=\frac{\mu }{2}(\mathrm{tr}\left(\mathit{C}\right)-3-2\ln \left(\sqrt{det\mathit{C}}\right))+\frac{\mathrm{\Lambda }}{4}(det\mathit{C}-1-2\ln \left(\sqrt{det\mathit{C}}\right))$$

(25)

where $\mathrm{\Lambda }$ and $\mu$ denote the first and second Lamé constant, respectively. The degradation functions read

$$f_{\mathrm{iso}}:={(1-\frac{\mathrm{tr}\left(\mathit{D}\right)}{3})}^{e_d}$$

(26)

and

$$f_{\mathrm{ani}}:={(1-\frac{\mathrm{tr}({\mathit{C}}^2\mathit{D})}{\mathrm{tr}({\mathit{C}}^2)})}^{f_d}$$

(27)

with the exponents $e_d$ and $f_d$ being additional material parameters introduced for further flexibility. In contrast to Reese et al. (2021), the anisotropic degradation function $f_{\mathrm{ani}}(\mathit{C},\mathit{D})$ is in this work formulated with respect to the right Cauchy–Green tensor $\mathit{C}$ instead of the Green–Lagrange strain tensor $\mathit{E}$ to avoid a division by zero after an unloading of the material.

Remark—fulfillment of damage growth criterion

This ansatz for formulating the elastic energy yields a general fulfillment of the damage growth criterion as derived in “Appendix A.1 General fulfillment of the damage growth criterion,” which may not be achieved by an isochoric–volumetric energy split (see “Appendix A.2 Isochoric violation of the damage growth criterion”) without considering a logarithmic strain energy formulation (cf. Holthusen et al., 2022a).

The isotropic damage hardening energy contains a nonlinear (cf. Reese et al., 2021) and a linear contribution

$${\psi }_d=r_d({\xi }_d+\frac{\exp -s_d{\xi }_d-1}{s_d})+\frac{1}{2}{H}_d{\xi }_d^2$$

(28)

with the damage hardening material parameters $r_d$ , $s_d$ , and $H_d$ . The additional damage hardening energy ${\psi }_h$ results in kinematic damage hardening. Analogously to Fassin et al. (2019a, 2019b); Holthusen et al. (2022a), the energy is formulated in terms of the eigenvalues $D_i$ of the damage tensor and ensures that these do not exceed a value of one. It reads

$${\psi }_h=K_h\sum _{i=1}^3(-\frac{{(1-D_i)}^{1-(1/n_h)}}{1-(1/n_h)}-D_i+\frac{1}{1-(1/n_h)})$$

(29)

with the material parameters $K_h$ and $n_h$ . The micromorphic energy contribution contains a summation over the number of nonlocal degrees of freedom $n_{\overline{d}}$ according to the size of the micromorphic tuple

$${\psi }_{\overline{d}}=\frac{1}{2}\sum _{i=1}^{n_{\overline{d}}}{H}_i{(d_i-{\overline{d}}_i)}^2+\frac{1}{2}\sum _{i=1}^{n_{\overline{d}}}{A}_i\mathrm{Grad}\left({\overline{d}}_i\right)\cdot \mathrm{Grad}\left({\overline{d}}_i\right)$$

(30)

and includes the micromorphic penalty parameters $H_i$ and the micromorphic gradient parameters $A_i$ .

The micromorphic tuples used in this work are investigated in van der Velden et al. (2024) and describe a full and two reduced regularizations of the damage tensor. The micromorphic tuple of model A reads

$$\begin{array}{rl}{\mathbf{d}}^{\mathit{A}}=& ((\mathrm{tr}(\mathit{D}{\mathit{M}}_1),\mathrm{tr}(\mathit{D}{\mathit{M}}_2),\mathrm{tr}(\mathit{D}{\mathit{M}}_3),\mathrm{tr}(\mathit{D}{\mathit{M}}_4),\mathrm{tr}(\mathit{D}{\mathit{M}}_5),\mathrm{tr}(\mathit{D}{\mathit{M}}_6))\end{array}$$

(31)

where the structural tensors are defined using the Cartesian basis vectors ${\mathit{e}}_1$ , ${\mathit{e}}_2$ , and ${\mathit{e}}_3$ as

$$\begin{array}{rl}{\mathit{M}}_1& ={\mathit{e}}_1\otimes {\mathit{e}}_1,{\mathit{M}}_2={\mathit{e}}_2\otimes {\mathit{e}}_2,{\mathit{M}}_3={\mathit{e}}_3\otimes {\mathit{e}}_3,\\ {\mathit{M}}_4& ={\mathit{e}}_1\otimes {\mathit{e}}_2,{\mathit{M}}_5={\mathit{e}}_1\otimes {\mathit{e}}_3,{\mathit{M}}_6={\mathit{e}}_2\otimes {\mathit{e}}_3.\end{array}$$

(32)

Due to the symmetry of $\mathit{D}$ and the properties of the trace operator and double contraction with $\text{tr}({\mathit{A}}^{\text{T}}\mathit{B})=\mathit{A}:\mathit{B}$ and $\text{sym}(\mathit{A}):\mathit{B}=\text{sym}(\mathit{A}):\text{sym}(\mathit{B})$ , the structural tensors for the micromorphic tuple of model A could alternatively be defined as

$$\begin{array}{rl}{\mathit{M}}_1& =\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_1),{\mathit{M}}_2=\mathrm{sym}({\mathit{e}}_2\otimes {\mathit{e}}_2),{\mathit{M}}_3=\mathrm{sym}({\mathit{e}}_3\otimes {\mathit{e}}_3),\\ {\mathit{M}}_4& =\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_2),{\mathit{M}}_5=\mathrm{sym}({\mathit{e}}_1\otimes {\mathit{e}}_3),{\mathit{M}}_6=\mathrm{sym}({\mathit{e}}_2\otimes {\mathit{e}}_3).\end{array}$$

(33)

Model A controls all six independent components of the second-order damage tensor by a nonlocal regularization field and, thus, serves with the full regularization as a reference solution for the reduced micromorphic tuples. Model B employs a reduced principal traces regularization with three nonlocal degrees of freedom and stems from Holthusen et al. (2022a) with

$${\mathbf{d}}^{\mathit{B}}=(\mathrm{tr}(\mathit{D}),\mathrm{tr}({\mathit{D}}^2),\mathrm{tr}({\mathit{D}}^3)).$$

(34)

Model C utilizes a reduced volumetric–deviatoric regularization with two nonlocal degrees of freedom as proposed by Holthusen et al. (2022b) and investigated in van der Velden et al. (2024) with

$${\mathbf{d}}^{\mathit{C}}=(\frac{\mathrm{tr}(\mathit{D})}{3},\mathrm{tr}(\mathrm{dev}{\left(\mathit{D}\right)}^2)).$$

(35)

Explicit thermodynamic conjugate driving forces

Considering the specific forms of the Helmholtz free energy and the micromorphic tuples given in Section “Specific Helmholtz free energies and micromorphic tuples,” the explicit forms of the thermodynamic conjugate driving forces are presented here according to the state laws and definitions of Section “Isothermal Clausius–Duhem inequality.” The second Piola–Kirchhoff stress follows from equations (14) and (24)

$$\mathit{S}=((1-k_{\mathrm{ani}})f_{\mathrm{iso}}+k_{\mathrm{ani}}{f}_{\mathrm{ani}}){\mathit{S}}_{\mathrm{NH}}+2k_{\mathrm{ani}}{\displaystyle \frac{\mathrm{\partial }{f}_{\mathrm{ani}}}{\mathrm{\partial }\mathit{C}}{\psi }_{\mathrm{NH}}}$$

(36)

with ${\mathit{S}}_{\mathrm{NH}}:=2\mathrm{\partial }{\psi }_{\mathrm{NH}}/\mathrm{\partial }\mathit{C}$ . As stated in Reese et al. (2021), the partial derivative $\mathrm{\partial }{f}_{\mathrm{ani}}/\mathrm{\partial }\mathit{C}$ should vanish in the undamaged state, that is, $\mathit{D}=\mathbf{0}$ , and completely damaged state, that is, $\mathit{D}=\mathit{I}$ . Due to the choice of a modified anisotropic degradation function in equation (27), the analytical solution of $\mathrm{\partial }{f}_{\mathrm{ani}}/\mathrm{\partial }\mathit{C}{|}_{\mathit{D}=\mathbf{0}}$ and $\mathrm{\partial }{f}_{\mathrm{ani}}/\mathrm{\partial }\mathit{C}{|}_{\mathit{D}=\mathit{I}}$ is presented in “Appendix A.3 Derivative of anisotropic degradation function.”

The elastic damage driving force follows from equations (12) and (24)

$${\mathit{Y}}_e=-\left((1-k_{\mathrm{ani}}){\displaystyle \frac{\mathrm{\partial }{f}_{\mathrm{iso}}}{\mathrm{\partial }\mathit{D}}+k_{\mathrm{ani}}{\displaystyle \frac{\mathrm{\partial }{f}_{\mathrm{ani}}}{\mathrm{\partial }\mathit{D}}}}\right){\psi }_{\mathrm{NH}}$$

(37)

where the partial derivatives of the degradation functions with respect to the damage tensor are given in equations (59) and (60). The kinematic damage hardening driving force is formulated in the eigensystem of the damage tensor, implemented analogously to Holthusen et al. (2022a), and follows from equations (12) and (29)

$${\mathit{Y}}_h=K_h\sum _{i=1}^3(\frac{1}{{(1-D_i)}^{1/n_h}}-1){\mathit{n}}_i^{\mathit{D}}\otimes {\mathit{n}}_i^{\mathit{D}}$$

(38)

where $D_i$ and ${\mathit{n}}_i^{\mathit{D}}$ denote the eigenvalues and eigenvectors of $\mathit{D}$ . The nonlocal damage driving force reads generally (see van der Velden et al., 2024) with the definition in equation (12)

$${\mathit{Y}}_{\overline{d}}=\sum _{i=1}^{n_{\overline{d}}}{H}_i(d_i-{\overline{d}}_i){\displaystyle \frac{\mathrm{\partial }{d}_i}{\mathrm{\partial }\mathit{D}}.}$$

(39)

The explicit forms follow with the definitions for the micromorphic tuples in equations (31), (34), and (35) and read for model A

$$\begin{array}{rl}{\mathit{Y}}_{\overline{d}}^{\mathit{A}}& =H_1(\mathrm{tr}(\mathit{D}{\mathit{M}}_1)-{\overline{d}}_1)\mathrm{sym}\left({\mathit{M}}_1\right)+H_2(\mathrm{tr}(\mathit{D}{\mathit{M}}_2)-{\overline{d}}_2)\mathrm{sym}\left({\mathit{M}}_2\right)\\ & +H_3(\mathrm{tr}(\mathit{D}{\mathit{M}}_3)-{\overline{d}}_3)\mathrm{sym}\left({\mathit{M}}_3\right)+H_4(\mathrm{tr}(\mathit{D}{\mathit{M}}_4)-{\overline{d}}_4)\mathrm{sym}\left({\mathit{M}}_4\right)\\ & +H_5(\mathrm{tr}(\mathit{D}{\mathit{M}}_5)-{\overline{d}}_5)\mathrm{sym}\left({\mathit{M}}_5\right)+H_6(\mathrm{tr}(\mathit{D}{\mathit{M}}_6)-{\overline{d}}_6)\mathrm{sym}\left({\mathit{M}}_6\right),\end{array}$$

(40)

for model B

$${\mathit{Y}}_{\overline{d}}^{\mathit{B}}=H_1(\mathrm{tr}(\mathit{D})-{\overline{d}}_1)\mathit{I}+H_2(\mathrm{tr}({\mathit{D}}^2)-{\overline{d}}_2)2\mathit{D}+H_3(\mathrm{tr}({\mathit{D}}^3)-{\overline{d}}_3)3{\mathit{D}}^2,$$

(41)

and for model C

$${\mathit{Y}}_{\overline{d}}^{\mathit{C}}=\frac{H_1}{3}(\frac{\mathrm{tr}(\mathit{D})}{3}-{\overline{d}}_1)\mathit{I}+H_2(\mathrm{tr}(\mathrm{dev}{\left(\mathit{D}\right)}^2)-{\overline{d}}_2)(2\mathit{D}-\frac{2}{3}\mathrm{tr}\left(\mathit{D}\right)\mathit{I}).$$

(42)

The isotropic damage hardening driving force follows from equations (12) and (28)

$$R_d=-[r_d(1-\exp (-s_d{\xi }_d))+H_d{\xi }_d]$$

(43)

and the generalized micromorphic stresses from equations (15), (16), and (30)

$${\left({\mathit{\xi }}_{0_i}\right)}_k=-H_k(d_k-{\overline{d}}_k),k\in \{1,\dots ,n_{\overline{d}}\},$$

(44)

$${\left({\mathbf{\Xi }}_{0_i}\right)}_k=A_k\mathrm{Grad}\left({\overline{d}}_k\right),k\in \{1,\dots ,n_{\overline{d}}\}.$$

(45)

Numerical aspects

The evolution equations, equations (21) to (23), are discretized in time by the backward Euler method. Moreover, the local residual to be solved in an inelastic damage step reads

$$\left(\begin{array}{c}{r}_1\\ {\mathit{R}}_2\end{array}\right):=\left(\begin{array}{c}3{\mathit{Y}}_{+}:\mathbb{A}:{\mathit{Y}}_{+}-(Y_0-R_d{)}^2\\ \mathit{D}-{\mathit{D}}_n-\mathrm{\Delta }{\gamma }_{d}3({\mathbb{P}}_{+}:\mathbb{A}:{\mathit{Y}}_{+}+{\mathit{Y}}_{+}:\mathbb{A}:{\mathbb{P}}_{+})/(2(Y_0-R_d))\end{array}\right)=\left(\begin{array}{c}0\\ \mathbf{0}\end{array}\right)$$

(46)

where the squared damage onset function ${\mathrm{\Phi }}_d^2$ is considered in $r_1$ . The time discretized evolution equation of the damage tensor, equation (21), is modified by replacing $\sqrt{3}\sqrt{{\mathit{Y}}_{+}:\mathbb{A}:{\mathit{Y}}_{+}}=Y_0-R_d$ in ${\mathit{R}}_2$ (see Challamel et al., 2005; Holthusen et al., 2022a). Furthermore, $\mathrm{\Delta }{\gamma }_d$ denotes the time discretized damage multiplier, and ${\mathbb{P}}_{+}$ denotes a positive projection tensor (see Holthusen et al., 2022a). The local residual is solved in combination with a line search algorithm that the interested reader may find in the provided material subroutines (https://doi.org/10.5281/zenodo.11171630).

Moreover, an artificial viscosity ${\eta }_v$ that serves to prevent snap-back in structural simulations is considered in the material formulation. It is implemented via a modification of the mechanical tangent $\hat{\mathbb{C}}=\mathrm{\partial }\hat{\mathit{S}}/\mathrm{\partial }\hat{\mathit{E}}$ and the mechanical stress $\hat{\mathit{S}}$ at the integration point level after the local iteration is converged. Here, $\widehat{(\bullet )}$ denotes the Nye notation of the corresponding tensorial quantity. The modified quantities read

$${\hat{\mathbb{C}}}_{\text{t}ot}={\hat{\mathbb{C}}}_{\text{m}ech}+{\hat{\mathbb{C}}}_{{\eta }_v}$$

(47)

$${\hat{\mathit{S}}}_{\text{t}ot}={\hat{\mathit{S}}}_{\text{m}ech}+{\hat{\mathit{S}}}_{{\eta }_v}$$

(48)

with

$${\hat{\mathbb{C}}}_{{\eta }_v}:=\frac{{\eta }_v}{\mathrm{\Delta }t}\left[\begin{array}{cccccc}2& & & & & \\ & 2& & & & \\ & & 2& & & \\ & & & 1& & \\ & & & & 1& \\ & & & & & 1\end{array}\right]$$

(49)

and

$${\hat{\mathit{S}}}_{{\eta }_v}={\hat{\mathbb{C}}}_{{\eta }_v}(\hat{\mathit{E}}-{\hat{\mathit{E}}}_n)$$

(50)

where $\mathrm{\Delta }t$ denotes the current time increment and the vector of the strain-like quantities reads in Nye notation $\hat{\mathit{E}}=[E_{xx},E_{yy},E_{zz},2E_{xy},2E_{xz},2E_{yz}{]}^{\text{T}}$ . For further details, we kindly refer to the provided material subroutines (https://doi.org/10.5281/zenodo.11171630).

Numerical examples

First, the numerical examples show the local study without using any gradient-extension of a single finite element that is loaded by uniaxial tension, uniaxial strain, simple shear, and torsion in Section “Single-element studies” to analyze the behavior at the material point level. Thereafter in Section “Asymmetrically notched specimen,” the gradient-extended finite element formulation is applied for the structural simulation of an asymmetrically notched specimen using models A, B, and C to confirm the accuracy and efficiency of the reduced volumetric–deviatoric regularization of model C. Additionally, this example is investigated using a local formulation without gradient-extension and also using further reduced regularizations based on a single component of the damage tensor. Finally in Section “Rotor blade,” model C is employed for the complex three-dimensional structural simulation of a pressure-loaded rotor blade specimen.

The material parameters are provided in Table 1, where the parameters of Set 1 are utilized for the single-element studies in Section “Single-element studies” and the parameters of Set 2 for the structural simulations in Sections “Asymmetrically notched specimen” and “Rotor blade.” The Taylor series sampling point $a_h$ is not introduced during the presentation of the constitutive modeling in Section “Constitutive modeling,” but is required for the implementation of the kinematic damage driving force ${\mathit{Y}}_h$ as elaborated in Holthusen et al. (2022a).

Table 1.

Material and numerical parameters.

Symbol	Material parameter	Set 1	Set 2	Unit
$\mathrm{\Lambda }$	First Lamé constant	5000	25,000	MPa
$\mu$	Second Lamé constant	7500	55,000	MPa
$k_{\mathrm{ani}}$	Degree damage anisotropy	0 – 1	1	–
$e_d$	Exponent $f_{\mathrm{iso}}$	2	2	–
$f_d$	Exponent $f_{\mathrm{ani}}$	1	1	–
$Y_0$	Initial damage threshold	10	2.5	MPa
$c_d$	Distortional hardening exponent	1	1	–
$H_d$	Linear isotropic hardening prefactor	1	1	MPa
$r_d$	Nonlinear isotropic hardening prefactor	10	5	MPa
$s_d$	Nonlinear isotropic hardening scaling factor	100	100	–
$K_h$	Kinematic hardening prefactor	0.1	0.1	MPa
$n_h$	Kinematic hardening exponent	2	2	–
$A_i$	Internal length scales	0	300 – 3000	MPa mm2

Symbol	Numerical parameter	Value	Value	Unit
$a_h$	Taylor series sampling point	0.999999	0.999999	–
$H_i$	Micromorphic penalty parameters	0	104	MPa
${\eta }_v$	Artificial viscosity	1	1	MPa s

The finite elements are eight-node hexahedrons with full integration that could be substituted in future investigations with the reduced integration-based elements of Barfusz et al. (2021a, 2021b, 2022). For the finite element simulations, we utilize the software FEAP (Taylor and Govindjee, 2020), for the finite element meshes of the rotor blade in Section “Rotor blade” the software HyperMesh (HyperWorks, 2022), and for processing the contour plots of the simulations the software ParaView (Ahrens et al., 2005).

Single-element studies

The studies of a single finite element serve for the investigation of the material behavior under a homogeneous loading state and, thus, employ its local formulation by neglecting the micromorphic regularization by setting $A_i=0[\text{MPa}{\text{mm}}^2]$ and $H_i=0[\text{MPa}]$ .

Uniaxial tension is applied in Figure 1(a) for different values of the damage anisotropy parameter with $k_{\mathrm{ani}}\in (0.00[-],0.33[-],0.67[-],1.00[-])$ , where $k_{\mathrm{ani}}=0.00[-]$ models purely isotropic and $k_{\mathrm{ani}}=1.00[-]$ purely anisotropic material degradation. The normalized force-displacement curves show an increasing maximum retention force followed by a steeper force reduction for an increase of the damage anisotropy parameter $k_{\mathrm{ani}}$ . A uniaxial strain state is applied in Figure 1(b), where the curves show a qualitatively similar behavior to Figure 1(a). However, quantitatively, the maximum retention force for $k_{\mathrm{ani}}=1.00[-]$ is $5.93[\mathrm{\%}]$ higher compared to uniaxial tension due to the constrained lateral contraction. Simple shear is applied in Figure 1(c) and yields a response qualitatively similar to uniaxial tension and uniaxial strain, but the maximum retention force for $k_{\mathrm{ani}}=1.00[-]$ is $34.27[\mathrm{\%}]$ smaller compared to uniaxial tension. Finally, a torsional load is applied in Figure 1(d) also for different values of the damage anisotropy parameter and confirms the previous observations of Figure 1(a) to (c) that for this model a higher material strength corresponds to higher values of damage anisotropy. However, it should be noted that torsion in Figure 1(d), unlike the loadings in Figure 1(a) to (c), does not represent a homogeneous strain state.

Figure 1.

Single-element studies for uniaxial tension, uniaxial strain, simple shear, and torsion. The force-displacement curves are normalized with respect to the maximum force for uniaxial tension for $k_{\mathrm{ani}}=1.00$ with $F_{\text{max}}=8.1251\times {10}^2[\text{N}]$ . The moment-twist curves are normalized with respect to the maximum moment for $k_{\mathrm{ani}}=1.00$ with $M_{\text{max}}=-2.1427\times {10}^2[\text{N}\text{mm}]$ : (a) uniaxial tension; (b) uniaxial strain; (c) simple shear; and (d) torsion.

Asymmetrically notched specimen

The first structural example considers an asymmetrically notched specimen in Figure 2 under plane strain conditions with the dimensions $l=100[\text{mm}]$ , $l_1=40[\text{mm}]$ , $l_2=20[\text{mm}]$ , $w=36[\text{mm}]$ , $r=5[\text{mm}]$ and a thickness of $1[\text{mm}]$ that was previously studied in e.g. Brepols et al., 2017. Analogously to van der Velden et al. (2024), this example is investigated using the gradient-extensions of models A, B, and C. The gradient parameter of model A is arbitrarily chosen as $A_i^{\mathit{A}}=1000[\text{MPa}{\text{mm}}^2]$ and the parameters of models B and C are identified as $A_i^{\mathit{B}}=300[\text{MPa}{\text{mm}}^2]$ and $A_i^{\mathit{C}}=3000[\text{MPa}{\text{mm}}^2]$ to obtain the same structural load bearing capacity.

Figure 2.

Geometry and boundary value problem for the asymmetrically notched specimen.

The force-displacement curves in Figure 3(a) to (c) shows the mesh convergence studies for models A, B, and C with the meshes stemming from Holthusen et al. (2022a). All models yield excellent convergence behavior with negligible differences between the results of the coarsest and finest mesh. A model comparison of the structural response with respect to the force-displacement curves is provided in Figure 3(d) and shows a fine agreement between model A with the full regularization and model C with the reduced volumetric–deviatoric regularization. For model B, the vertical drop in the force-displacement curve is shifted significantly to the right compared to models A and C and, thus, entails a larger amount of dissipated energy in the failure process, which is in line with the results of van der Velden et al. (2024).

Figure 3.

Force-displacement curves for the mesh convergence study and the model comparison of the asymmetrically notched specimen. The forces are normalized with respect to the maximum force of model C (13,955 elements) with $F_{\text{max}}=4.4031\times {10}^4[\text{N}]$ . In Figure 3(a) to (c), the curves obtained with the coarse and fine meshes are essentially congruent. In Figure 3(d), the squares indicate the snapshots of Figure 4 and the circles indicate the snapshots of Figure 5: (a) model A; (b) model B; (c) model C; and (d) model comparison (13,955 elements).

A comparison of the damage contour plots of models A, B, and C for the asymmetrically notched specimen is given in Figure 4. The crack width and damage-affected zone of model B (Figure 4(b)) are thicker than those of models A and C (Figure 4(a) and (c)) and, thereby, yield a higher energy dissipation. Nevertheless, close agreement between the damage patterns of models A and C is observed.

Figure 4.

Damage contour plots for the asymmetrically notched specimen with models A, B, and C (13,955 elements) at the end of the simulation: (a) model A; (b) model B; and (c) model C.

The evolution of the normal ( $D_{xx}$ , $D_{yy}$ ) and shear ( $D_{xy}$ ) components of the damage tensor for model C is presented in Figure 5, where the position of the snapshots is indicated by the circles in Figure 3(d). Damage initiates at both notches (Figure 5(a)), then it forms a damage shear band zone (Figure 5(b)), in which finally the crack forms (Figure 5(c)).

Figure 5.

Damage evolution contour plots for the asymmetrically notched specimen with model C (13955 elements): (a) initial; (b) intermediate; and (c) final.

Thereafter, the behavior of the local anisotropic damage without regularization is investigated for the asymmetrically notched specimen. Analogously to Fassin et al. (2019a), the micromorphic parameters are set to $A_i=0[\text{MPa}{\text{mm}}^2]$ and $H_i=0[\text{MPa}]$ and, moreover, zero Dirichlet conditions are applied to all nonlocal degrees of freedom, that is, $\overline{\mathbf{d}}=\mathbf{0}$ . Figure 6 shows the corresponding force-displacement curves for different finite element mesh discretizations. Using the local formulation, a converged solution is neither achieved with respect to the maximum load-bearing capacity of the specimen nor with respect to the amount of dissipated energy.

Figure 6.

Force-displacement curves for the local damage model for the asymmetrically notched specimen. The forces are normalized with respect to the maximum force of the finest mesh (13,955 elements) with $F_{\text{max}}=1.5477\times {10}^4[\text{N}]$ .

Furthermore, the damage contour plots for the study of the local anisotropic damage model are presented in Figure 7. For all meshes, two cracks form at the notches and propagate horizontally through the specimen, where all cracks localize into a single row of elements. They do not exhibit a tendency to form a shear crack and to coalesce, which contradicts the results of the gradient-extended solution in Figures 3 to 5 and the experimental investigations of, for example, Ambati et al. (2016). Hence, the utilized artificial viscosity ${\eta }_v=1[\text{MPa}\text{s}]$ does not yield regularizing effects, and the excellent mesh convergence in Figure 3(a) to (c) is related to the gradient-extensions of models A, B, and C.

Figure 7.

Damage contour plots with the local damage model for the asymmetrically notched specimen at $u/l\times {\text{10}}^2=0.4[-]$ . The full specimen is shown with the coarsest and finest mesh in Figure 7(a) and (b) and the center of the specimen is shown for all meshes in Figure 7(c) to (h). For this study, the damage variables are averaged over all integration points per element. (a) 1624 elements; (b) 13,955 elements; (c) 1624 elements; (d) 3592 elements; (e) 6651 elements; (f) 9667 elements; (g) 12,704 elements; and (h) 13,955 elements.

Next, we investigate the possibility of a further reduction of the micromorphic tuple by using only a single degree of freedom. Therefore, model A is again utilized for the simulation of the asymmetrically notched specimen, but each time only a single component of the micromorphic tuple is activated. For example, “ ${\overline{d}}_1$ active” refers to the micromorphic gradient parameters $A_1^{\mathit{A}}=1000[\text{MPa}{\text{mm}}^2]$ , $A_2^{\mathit{A}}=A_3^{\mathit{A}}=A_4^{\mathit{A}}=A_5^{\mathit{A}}=A_6^{\mathit{A}}=0[\text{MPa}{\text{mm}}^2]$ , the micromorphic penalty parameters $H_1^{\mathit{A}}={10}^4[\text{MPa}]$ , $H_2^{\mathit{A}}=H_3^{\mathit{A}}=H_4^{\mathit{A}}=H_5^{\mathit{A}}=H_6^{\mathit{A}}=0[\text{MPa}]$ , and the Dirichlet boundary conditions ${\overline{d}}_2={\overline{d}}_3={\overline{d}}_4={\overline{d}}_5={\overline{d}}_6=0[-]$ . For model A, ${\overline{d}}_1$ is associated with the regularization of $D_{xx}$ , ${\overline{d}}_2$ with $D_{yy}$ , ${\overline{d}}_3$ with $D_{zz}$ , ${\overline{d}}_4$ with $D_{xy}$ , ${\overline{d}}_5$ with $D_{xz}$ , and ${\overline{d}}_6$ with $D_{yz}$ .

Figure 8 shows the corresponding force-displacement curves for the regularization of a single component of the damage tensor as well as the reference solution where all six independent components are regularized (from Figure 3(a)). The regularization of a single component does not suffice for the regularization of the shear crack, since all solutions obtained by a single-component regularization underestimate the maximum force carried by the specimen. The regularization of the normal component $D_{xx}$ by ${\overline{d}}_1$ , that is, the damage tensor component in the loading direction, yields an underestimation of $16.86[\mathrm{\%}]$ . The regularization of the normal components $D_{yy}$ and $D_{zz}$ by ${\overline{d}}_2$ and ${\overline{d}}_3$ yield an underestimation of $19.77[\mathrm{\%}]$ and $19.55[\mathrm{\%}]$ and the regularization of the shear components $D_{xy}$ , $D_{xz}$ , and $D_{yz}$ by ${\overline{d}}_4$ , ${\overline{d}}_5$ , and ${\overline{d}}_6$ each yields an underestimation of $62.80[\mathrm{\%}]$ .

Figure 8.

Force-displacement curves for the damage model with regularization of single components of the damage tensor using model A for the asymmetrically notched specimen. The forces are normalized with respect to the maximum force of the reference regularization using the mesh with 1624 elements with $F_{\text{max}}=4.4114\times {10}^4[\text{N}]$ .

The damage contour plots at $u/l\times {\text{10}}^2=0.9[-]$ and $u/l\times {\text{10}}^2=2.0[-]$ are provided in Figures 9 and 10. According to the force-displacement curve of the reference solution in Figure 3(a), an undamaged state exists at $u/l\times {\text{10}}^2=0.9[-]$ and the completely damage state at $u/l\times {\text{10}}^2=2.0[-]$ . The regularization of single damage tensor components influences the specific component’s evolution significantly, but overall fails to obtain the reference solution. If, for example, ${\overline{d}}_1$ is active, it diminishes the damage initiation of $D_{xx}$ , whereas $D_{yy}$ already shows signs of localization at the beginning of the failure process (Figure 9(a)). Upon further loading, $D_{xx}$ still displays a diffuse damage zone at the edges of the crack, but eventually the crack localizes and propagates horizontally instead of slantingly (Figure 10(a)). The evolution of $D_{yy}$ at further loading yields the same crack pattern, but exhibits a thicker fully damaged zone with little diffusive character at its edges. The activation of ${\overline{d}}_2$ instead of ${\overline{d}}_1$ yields a converse behavior of $D_{xx}$ and $D_{yy}$ (Figures 9(b) and 10(b)). The activation of ${\overline{d}}_3$ with a regularization of $D_{zz}$ yields a pronounced concentration behavior for $D_{xx}$ and $D_{yy}$ (Figures 9(c) and 10(c)) analogously to the previously unregularized quantities ( $D_{yy}$ in Figures 9(a) and 10(a) and $D_{xx}$ in Figures 9(b) and 10(b)). The activation of ${\overline{d}}_4$ , ${\overline{d}}_5$ , and ${\overline{d}}_6$ with a regularization of $D_{xy}$ , $D_{xz}$ , and $D_{yz}$ yields a localization of the normal components $D_{xx}$ and $D_{yy}$ into a single row of elements (Figures 9(d) and 10(d)).

Figure 9.

Damage contour plots for the damage model with regularization of single components of the damage tensor using model A for the asymmetrically notched specimen (1624 elements) at $u/l\times {\text{10}}^2=0.9[-]$ in Figure 8. For this study, the damage variables are averaged over all integration points per element: (a) ${\overline{d}}_1$ active; (b) ${\overline{d}}_2$ active; (c) ${\overline{d}}_3$ active; (d) ${\overline{d}}_4$ active; and (e) reference.

Figure 10.

Damage contour plots for the damage model with regularization of single components of the damage tensor using model A for the asymmetrically notched specimen (1624 elements) at $u/l\times {\text{10}}^2=2.0[-]$ in Figure 8. For this study, the damage variables are averaged over all integration points per element: (a) ${\overline{d}}_1$ active; ${\overline{d}}_2$ active; (c) ${\overline{d}}_3$ active; (d) ${\overline{d}}_4$ active; and (e) reference.

Rotor blade

After confirming the accuracy of the reduced volumetric–deviatoric regularization of model C in Section “Asymmetrically notched specimen,” it is applied for the three-dimensional simulation of a rotor blade specimen. This aims at studying the performance of the generic anisotropic damage framework, here specified for a neo-Hookean material regularized by model C, to predict damage evolution on a complex structural example. The number of nodes increases significantly in three-dimensional simulations and, thus, only model C is employed for the simulation with two nonlocal degrees of freedom.

The geometry is inspired by previous works in van der Velden et al. (2023), where the electrochemical machining process is simulated for the manufacturing of a rotor blade (Figure 11(a)) that can be assembled to an entire rotor (Figure 11(b)).

Figure 11.

Motivation for the rotor blade specimen. In Figure 11(a), a single blade is manufactured by electrochemical machining from a solid metal workpiece (dark gray). The tool (light gray) defines the shape and the electrolyte (blue) serves as an electrical conductor for a direct current. (a) Blade production (van der Velden et al., 2023); (b) rotor.

Figure 12 provides the side and the bottom view of the geometry and boundary value problem where the geometrical dimensions, angles, and positions read $l=20[\text{mm}]$ , $r_1=2[\text{mm}]$ , $r_2=20[\text{mm}]$ , $r_3=25[\text{mm}]$ , $r_4=0.5[\text{mm}]$ , $w_1=13.5[\text{mm}]$ , $w_2=4.882[\text{mm}]$ , ${\phi }_1=60{[}^{\circ }]$ , ${\phi }_2=90{[}^{\circ }]$ , ${\phi }_3=100.137{[}^{\circ }]$ , ${\phi }_4=117.444{[}^{\circ }]$ , $P_1=(-9.0|-15.5884|0.0)[\text{mm}]$ , $P_2=(-13.5|-23.3827|0.0)[\text{mm}]$ . The rear edge of the rotor blade is clamped at $z=0[\text{mm}]$ and the loading case of a pressure load applied from the bottom side is considered. This simplified loading neglects the lift generation on the upper (suction) side of the rotor blade. Furthermore, centrifugal effects due to the blade’s rotation are not included in the simulation.

Figure 12.

Geometry and boundary value problem for the rotor blade specimen. Side view (x–y) and bottom view (x–z).

In Figure 13, the finite element meshes are shown, which are refined toward the clamped edge and toward the side with the smaller radius due to results of preliminary studies that revealed damage initiation in these regions.

Figure 13.

Finite element meshes for the rotor blade specimen (bottom view): (a) coarsest mesh (21,560 elements) and (b) finest mesh (116,565 elements).

The normalized force-displacement curves in Figure 14 depict the sum of the forces in the y-direction at the clamped edge over the deflection in the y-direction of point A (see Figure 12) that is located at position $P_A=(-13.5|1.6173|20.0)[\text{mm}]$ . The mesh convergence study yields a close agreement in the force-displacement curves between the results of all meshes. Moreover, the coarsest mesh with 21,560 elements underestimates the structural load-bearing capacity of the rotor blade that is obtained with the finest mesh with 116,565 elements by just $0.30[\mathrm{\%}]$ .

Figure 14.

Force-displacement curves for the mesh convergence study of the rotor blade specimen using model C. The forces are normalized with respect to the maximum force of the finest mesh (116,565 elements) with $F_{\text{max}}=1.5738\times {10}^4[\text{N}]$ . The orange boxes indicate the points of comparison for mesh convergence in Figures 15 and 16. The black circles indicate the points of evaluation for the damage evolution in Figures 17 and 18.

The excellent coarse mesh accuracy, which has been observed with respect to the load-bearing capacity, is also reflected in the comparison of the damage contour plots in Figures 15 (bottom view) and 16 (top view), where high agreement between the normal components of the damage tensor obtained with the coarsest and finest mesh is observed. Figure 16 reveals the concentrated evolution of the damage tensor components $D_{xx}$ and $D_{yy}$ in the middle of the upper side of the blade at the clamped end. In Figure 15, the component $D_{zz}$ predominantly evolves at the edges of the lower side at the clamped end. The coarse mesh, in Figures 15(a) and 16(a), predicts these damage patterns accurately and the fine mesh, in Figures 15(b) and 16(b), corroborates these results.

Figure 15.

Contour plots for the mesh convergence of the normal components of the damage tensor for the rotor blade specimen (bottom view): (a) 21,560 elements and (b) 116,565 elements.

Figure 16.

Contour plots for the mesh convergence of the normal components of the damage tensor for the rotor blade specimen (top view): (a) 21,560 elements and (b) 116,565 elements.

The damage evolution process is presented in Figures 17 and 18 from damage initiation (Figures 17(a) and 18(a)) to crack opening (Figures 17(d) and 18(d)) on the deformed configuration. For all normal components, damage initiates at the clamped edge (Figure 17(a)) and, upon further deformation, progresses with a diffuse damage zone in the z-direction (Figure 17(b) and (c)). Finally, crack opening occurs by a significant concentration of $D_{zz}$ at the intersection of the free outer edge with radius $r_4$ and the clamped edge (Figure 17(d)). The corresponding evolution of the shear components is presented in Figure 18 and shows that their peak values occur at the clamped edge.

Figure 17.

Contour plots for the evolution of the normal components of the damage tensor for the rotor blade specimen (side view). The contours are plotted on the deformed configuration and the opaque solid shapes indicate the reference configuration: (a) damage initiation; (b) intermediate 1; (c) intermediate 2; and (d) crack opening.

Figure 18.

Contour plots for the evolution of the shear components of the damage tensor for the rotor blade specimen (side view). The contours are plotted on the deformed configuration and the opaque solid shapes indicate the reference configuration: (a) damage initiation; (b) intermediate 1; (c) intermediate 2; and (d) crack opening.

Summary of the numerical results

The most important findings of the numerical examples with the local formulation, model A (full regularization), model B (reduced principal traces regularization), and model C (reduced volumetric–deviatoric regularization) include:

The material strength of the local model increases with the degree of damage anisotropy for this specific choice of the neo-Hookean elastic energy formulation (Figure 1).

Models A, B, and C prove excellent coarse mesh accuracy (Figures 3(a) to (c), and 14).

Models A and C yield high consistency in the structural force-displacements curves (Figure 3(d)) and damage contour plots (Figure 4(a) and (c)), while Model B yields a delayed force reduction (Figure 3(d)) with higher energy dissipation and a thicker damage zone (Figure 4(b)).

Utilizing the local anisotropic damage model without gradient-extension results in a false crack path prediction (Figure 7).

A single-component regularization of the damage tensor proves to be inadequate and yields different degrees of crack localization (Figures 8 to 10).

With two nonlocal degrees of freedom, model C is an efficient and accurate formulation with extraordinary coarse mesh accuracy (Figures 14 to 18).

Conclusion

In this work, we successfully introduced a generic framework for nonlocal anisotropic damage at finite strains. Due to the design of the degradation functions, arbitrarily established hyperelastic finite strain energy formulations can be incorporated into the anisotropic damage model. Furthermore, the generic micromorphic gradient-extension provides a versatile regularization that can prevent any desired number of local quantities from localization.

Initially, the behavior of a specific local anisotropic damage model utilizing a neo-Hookean energy was investigated in single-element studies and yielded an increased material strength for an increase in the degree of damage anisotropy. Thereafter, the model was applied with different gradient-extensions for the simulation of an asymmetrically notched specimen. It confirmed the accurate regularization capabilities of the volumetric–deviatoric regularization stemming from previous works for a new local anisotropic damage model. Moreover, the regularization of single components of the damage tensor was studied for the asymmetrically notched specimen and yielded localization regardless of the selection of the gradient-extended quantity. Eventually, the anisotropic damage model with the volumetric–deviatoric regularization was used for the simulation of a rotor blade specimen that was subjected to a pressure load.

Subsequent works could investigate the anisotropic damage behavior of further hyperelastic finite strain energy formulations utilizing the presented framework by straightforwardly replacing the neo-Hookean energy in the elastic term. Additionally, future studies should include unloading as well as reversion of the loadings in the numerical examples along with experimental comparisons.