# On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states

F. Vicentini, Camilla Zolesi, Pietro Carrara, Corrado Maurini, Laura De Lorenzis

► To cite this version:

F. Vicentini, Camilla Zolesi, Pietro Carrara, Corrado Maurini, Laura De Lorenzis. On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states. International Journal of Fracture, 2024, 10.1007/s10704-024-00763-w. hal-04231075v2

HAL Id: hal-04231075

https://hal.sorbonne-universite.fr/hal-04231075v2

Submitted on 23 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states

F. Vicentini, Camilla Zolesi, Pietro Carrara, Corrado Maurini, Laura De Lorenzis

► To cite this version:

F. Vicentini, Camilla Zolesi, Pietro Carrara, Corrado Maurini, Laura De Lorenzis. On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states. International Journal of Fracture, 2024, 10.1007/s10704-024-00763-w. hal-04231075v2

HAL Id: hal-04231075

https://hal.sorbonne-universite.fr/hal-04231075v2

Submitted on 23 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states

F. Vicentini · C. Zolesi · P. Carrara · C. Maurini · L. De Lorenzis

Received: 6 October 2023 / Accepted: 4 January 2024
© The Author(s) 2024

**Abstract** Phase-field models of brittle fracture are typically endowed with a decomposition of the elastic strain energy density in order to realistically describe fracture under multi-axial stress states. In this contribution, we identify the essential requirements for this decomposition to correctly describe both nucleation and propagation of cracks. Discussing the evolution of the elastic domains in the strain and stress spaces as damage evolves, we highlight the links between the nucleation and propagation conditions and the modulation of the elastic energy with the phase-field variable. In light of the identified requirements, we review some of the existing energy decompositions, showcasing their merits and limitations, and conclude that none of them is able to fulfil all requirements. As a partial remedy to this outcome, we propose a new energy decomposition, denoted as star-convex model, which involves a minimal modification of the volumetric-deviatoric

decomposition. Predictions of the star-convex model are compared with those of the existing models with different numerical tests encompassing both nucleation and propagation.

**Keywords** Nucleation · Propagation · Fracture · Multi-axial · Phase-field

## 1 Introduction

Francfort and Marigo revisited Griffith's criterion from a variational perspective by casting it into a global energy minimization framework; thus, they formulated a general principle circumventing the need for ad-hoc criteria to handle arbitrarily complex crack topologies in 2D or 3D (Francfort and Marigo 1998). The direct application of the ensuing (so-called *sharp-crack*) model is limited by the difficulty to handle displacement jumps in the numerical setting. Phase-field approaches come into play as a regularization of the sharp-crack model (Bourdin et al. 2000; Ambrosio 1992) and offer a smeared description of the crack as the localization of an auxiliary variable, i.e., the phase-field variable. The regularized model is prone to a simple numerical treatment using standard finite element discretizations with smooth basis functions; it introduces a regularization length $\ell$, which defines the typical width of the localization bands representing the approximation of the cracks.

F. Vicentini and C. Zolesi have contributed equally.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10704-024-00763-w.

F. Vicentini · P. Carrara · L. De Lorenzis
Computational Mechanics Group, Eidgenössische Technische Hochschule Zürich, Tannenstrasse 3, 8092 Zürich, Switzerland

C. Zolesi · C. Maurini (✉)
CNRS, Institut Jean Le Rond d'Alembert, Sorbonne Université, UMR 7190, 75005 Paris, France
e-mail: corrado.maurini@sorbonne-universite.fr
# On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states

F. Vicentini · C. Zolesi · P. Carrara · C. Maurini · L. De Lorenzis

Received: 6 October 2023 / Accepted: 4 January 2024
© The Author(s) 2024

**Abstract** Phase-field models of brittle fracture are typically endowed with a decomposition of the elastic strain energy density in order to realistically describe fracture under multi-axial stress states. In this contribution, we identify the essential requirements for this decomposition to correctly describe both nucleation and propagation of cracks. Discussing the evolution of the elastic domains in the strain and stress spaces as damage evolves, we highlight the links between the nucleation and propagation conditions and the modulation of the elastic energy with the phase-field variable. In light of the identified requirements, we review some of the existing energy decompositions, showcasing their merits and limitations, and conclude that none of them is able to fulfil all requirements. As a partial remedy to this outcome, we propose a new energy decomposition, denoted as star-convex model, which involves a minimal modification of the volumetric-deviatoric

decomposition. Predictions of the star-convex model are compared with those of the existing models with different numerical tests encompassing both nucleation and propagation.

**Keywords** Nucleation · Propagation · Fracture · Multi-axial · Phase-field

## 1 Introduction

Francfort and Marigo revisited Griffith's criterion from a variational perspective by casting it into a global energy minimization framework; thus, they formulated a general principle circumventing the need for ad-hoc criteria to handle arbitrarily complex crack topologies in 2D or 3D (Francfort and Marigo 1998). The direct application of the ensuing (so-called *sharp-crack*) model is limited by the difficulty to handle displacement jumps in the numerical setting. Phase-field approaches come into play as a regularization of the sharp-crack model (Bourdin et al. 2000; Ambrosio 1992) and offer a smeared description of the crack as the localization of an auxiliary variable, i.e., the phase-field variable. The regularized model is prone to a simple numerical treatment using standard finite element discretizations with smooth basis functions; it introduces a regularization length $\ell$, which defines the typical width of the localization bands representing the approximation of the cracks.

F. Vicentini and C. Zolesi have contributed equally.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10704-024-00763-w.

F. Vicentini · P. Carrara · L. De Lorenzis
Computational Mechanics Group, Eidgenössische Technische Hochschule Zürich, Tannenstrasse 3, 8092 Zürich, Switzerland

C. Zolesi · C. Maurini (✉)
CNRS, Institut Jean Le Rond d'Alembert, Sorbonne Université, UMR 7190, 75005 Paris, France
e-mail: corrado.maurini@sorbonne-universite.fr
Phase-field models are naturally able to predict the evolution of damage from a pristine material, which can be regarded as the capability to predict *nucleation* of cracks. From a theoretical perspective, this ability can be justified by treating the auxiliary regularization variable as a *damage* variable and interpreting the model as a gradient damage model endowed with a finite internal length *ℓ* (Pham et al. 2011b). Gradient damage modeling in the quasi-static rate-independent setting defines the system evolution through the three principles of *stability*, *irreversibility* and *energy balance* (Marigo et al. 2016). In particular, the stability requirement is equivalent to *local* energy minimization, in contrast to the starting sharp-crack problem which instead calls for *global* energy minimization. Evolution following local energy minima is not only numerically convenient when dealing with large-scale problems, but also more physically reasonable as it does not entail crossing energy barriers (Bourdin et al. 2008; Negri 2010; Baldelli and Maurini 2021).

Following the three principles of gradient dam-
age model evolution, localization of damage natu-
rally occurs upon the loss of stability of the homoge-
neous solution whose threshold depends on the length
ℓ (Pham et al. 2011b). Thus, the internal length can
be tuned to calibrate the uniaxial tensile strength of
the material (Tanné et al. 2018; Wu et al. 2017). On
the other hand, once the elastic parameters, the frac-
ture toughness of the material and ℓ are fixed, the stan-
dard phase-field model does not offer any additional
degree of freedom to calibrate the strength related to
scenarios different from uniaxial tension, such as the
compressive strength or the shear strength, see the dis-
cussion in (Kumar et al. 2020; De Lorenzis and Mau-
rini 2022). Flexibility in defining the strength envelope
represents the main problem of nucleation under multi-
axial stress states. As for crack propagation under such
stress states, the most relevant problem is to model uni-
lateral contact appropriately (Amor et al. 2009).

There is a wealth of literature proposing phase-
field models for brittle fracture under multi-axial stress
states. The performance of these models is demon-
strated through different examples, making their com-
parison not immediate and their strength and limita-
tions not obvious. Among these contributions, some
preserve a variational nature; the most popular varia-
tional solutions (Amor et al. 2009; Miehe et al. 2010)
are based on elastic energy decompositions. This idea
is adopted in (Freddi and Royer-Carfagni 2010, 2011;

De Lorenzis and Maurini 2022; Navidtehrani et al. 2022), justified through structured deformation theory. Other contributions are inspired by anisotropic materials (van Dijk et al. 2020; He and Shao 2019; Vu et al. 2022), propose cohesive fracture Lorentz (2017), introduce plasticity (Ulloa et al. 2022; You et al. 2020; Fei and Choo 2021) or propose an explicit treatment of the crack direction (Hakimzadeh et al. 2022; Steinke and Kaliske 2019; Storm et al. 2020). Unfortunately, the available variational models implicitly prioritize either nucleation or propagation, and (as we will show in this paper, at least for the most popular ones) none of them can describe both aspects correctly without introducing excessive complexity with respect to the original model. Other models step outside the variational framework; in this way, they more easily achieve the needed flexibility to handle nucleation and propagation, but only at the cost of giving up the theoretical and practical advantages of the variational structure (Kumar et al. 2020; Miehe et al. 2015; Zhang et al. 2017; Feng and Li 2022; Wang et al. 2019; Abrari Vajari et al. 2023). The variational principle naturally comes with a simple stability concept and it allows for the use of the mathematical tools of calculus of variations to discuss the existence of solutions and to study the asymptotic Γ-convergence behavior. Additionally, the finite element tangent stiffness matrix stemming from a variational formulation is automatically symmetric, which gives important advantages in numerics. Moreover, for constitutive models with one scalar damage internal variable, the Drucker-Ilyushin postulate is satisfied only if the damage criterion is derived from a variational formulation where the strain work is a state function (Marigo 1989) and the thermodynamic consistency of non-variational models is not granted.

To date, a study that highlights advantages and drawbacks of different energy decompositions based on a consistent set of criteria and common examples is still lacking. In this contribution, we first define these criteria; subsequently, we perform a systematic review of some available models, highlighting their performance in relation to the defined criteria. We find out that none of the examined energy decompositions is able to satisfy all the proposed criteria, as each of these decompositions performs well with respect to either nucleation or propagation of cracks. As a remedy, we propose a new model, which we denote as *star-convex model*; this model is still based on an energy decomposition, but it is specifically designed to satisfy the desired
Phase-field models are naturally able to predict the evolution of damage from a pristine material, which can be regarded as the capability to predict *nucleation* of cracks. From a theoretical perspective, this ability can be justified by treating the auxiliary regularization variable as a *damage* variable and interpreting the model as a gradient damage model endowed with a finite internal length *ℓ* (Pham et al. 2011b). Gradient damage modeling in the quasi-static rate-independent setting defines the system evolution through the three principles of *stability*, *irreversibility* and *energy balance* (Marigo et al. 2016). In particular, the stability requirement is equivalent to *local* energy minimization, in contrast to the starting sharp-crack problem which instead calls for *global* energy minimization. Evolution following local energy minima is not only numerically convenient when dealing with large-scale problems, but also more physically reasonable as it does not entail crossing energy barriers (Bourdin et al. 2008; Negri 2010; Baldelli and Maurini 2021).

Following the three principles of gradient dam-
age model evolution, localization of damage natu-
rally occurs upon the loss of stability of the homoge-
neous solution whose threshold depends on the length
ℓ (Pham et al. 2011b). Thus, the internal length can
be tuned to calibrate the uniaxial tensile strength of
the material (Tanné et al. 2018; Wu et al. 2017). On
the other hand, once the elastic parameters, the frac-
ture toughness of the material and ℓ are fixed, the stan-
dard phase-field model does not offer any additional
degree of freedom to calibrate the strength related to
scenarios different from uniaxial tension, such as the
compressive strength or the shear strength, see the dis-
cussion in (Kumar et al. 2020; De Lorenzis and Mau-
rini 2022). Flexibility in defining the strength envelope
represents the main problem of nucleation under multi-
axial stress states. As for crack propagation under such
stress states, the most relevant problem is to model uni-
lateral contact appropriately (Amor et al. 2009).

There is a wealth of literature proposing phase-
field models for brittle fracture under multi-axial stress
states. The performance of these models is demon-
strated through different examples, making their com-
parison not immediate and their strength and limita-
tions not obvious. Among these contributions, some
preserve a variational nature; the most popular varia-
tional solutions (Amor et al. 2009; Miehe et al. 2010)
are based on elastic energy decompositions. This idea
is adopted in (Freddi and Royer-Carfagni 2010, 2011;

De Lorenzis and Maurini 2022; Navidtehrani et al. 2022), justified through structured deformation theory. Other contributions are inspired by anisotropic materials (van Dijk et al. 2020; He and Shao 2019; Vu et al. 2022), propose cohesive fracture Lorentz (2017), introduce plasticity (Ulloa et al. 2022; You et al. 2020; Fei and Choo 2021) or propose an explicit treatment of the crack direction (Hakimzadeh et al. 2022; Steinke and Kaliske 2019; Storm et al. 2020). Unfortunately, the available variational models implicitly prioritize either nucleation or propagation, and (as we will show in this paper, at least for the most popular ones) none of them can describe both aspects correctly without introducing excessive complexity with respect to the original model. Other models step outside the variational framework; in this way, they more easily achieve the needed flexibility to handle nucleation and propagation, but only at the cost of giving up the theoretical and practical advantages of the variational structure (Kumar et al. 2020; Miehe et al. 2015; Zhang et al. 2017; Feng and Li 2022; Wang et al. 2019; Abrari Vajari et al. 2023). The variational principle naturally comes with a simple stability concept and it allows for the use of the mathematical tools of calculus of variations to discuss the existence of solutions and to study the asymptotic Γ-convergence behavior. Additionally, the finite element tangent stiffness matrix stemming from a variational formulation is automatically symmetric, which gives important advantages in numerics. Moreover, for constitutive models with one scalar damage internal variable, the Drucker-Ilyushin postulate is satisfied only if the damage criterion is derived from a variational formulation where the strain work is a state function (Marigo 1989) and the thermodynamic consistency of non-variational models is not granted.

To date, a study that highlights advantages and drawbacks of different energy decompositions based on a consistent set of criteria and common examples is still lacking. In this contribution, we first define these criteria; subsequently, we perform a systematic review of some available models, highlighting their performance in relation to the defined criteria. We find out that none of the examined energy decompositions is able to satisfy all the proposed criteria, as each of these decompositions performs well with respect to either nucleation or propagation of cracks. As a remedy, we propose a new model, which we denote as *star-convex model*; this model is still based on an energy decomposition, but it is specifically designed to satisfy the desired
requirements for both nucleation and propagation. All the models that we consider, including our newly proposed energy decomposition, are of phenomenological nature. In other words, they are not grounded on microscopically motivated considerations and do not aim at reproducing fine-scale features of the fracture behavior Marigo and Kazymyrenko (2019), which may well be specific of single materials. Rather, we aim at consistently reproducing the key macroscopic features of the fracture behavior under multiaxial stress states, namely the strength surface and the transmission of stresses across fully developed cracks, that are common to a wide class of brittle materials.

The paper is structured as follows. Section 2 offers a brief review on the basic ingredients of standard phase-field modeling of brittle fracture, shows its limitations and formulates the requirements for an ideal model to describe fracture under multi-axial stress states. In light of the defined requirements, Sect. 3 reports a review of some of the available models, all based on decompositions of the elastic strain energy density. In Sect. 4, the novel star-convex model is introduced. Section 5 is dedicated to numerical experiments that showcase the advantages and limitations of the models previously illustrated, including the new one. The main conclusions are drawn in Sect. 6.

As follows, we summarize the notation and some useful relations. Vectors and second-order tensors are both denoted by boldface fonts, e.g., **u** and **σ** for the displacement vector and stress tensor, respectively. For the standard orthogonal decomposition of second-order tensors in volumetric and deviatoric parts we adopt the following notation (exemplified on **σ**)

$$ 
\begin{aligned}
\boldsymbol{\sigma} &= \boldsymbol{\sigma}_{\text{vol}} + \boldsymbol{\sigma}_{\text{dev}}, \quad \boldsymbol{\sigma}_{\text{vol}} = \frac{\text{tr}(\boldsymbol{\sigma})}{n} \mathbf{I}, \\
\boldsymbol{\sigma}_{\text{dev}} &= \boldsymbol{\sigma} - \frac{\text{tr}(\boldsymbol{\sigma})}{n} \mathbf{I}, \quad \boldsymbol{\sigma}_{\text{vol}} \cdot \boldsymbol{\sigma}_{\text{dev}} = 0,
\end{aligned}
\quad (1) $$

where **I** is the second-order identity tensor, *n* is the number of space dimensions, and the dot denotes the inner product. Note that for problems in plane-strain/stress conditions it is *n* = 3. For an isotropic elastic undamaged material with Young's modulus *E*₀ and Poisson's ratio *v*₀, we denote by (λ₀, μ₀, κ₀) the Lamé constants and the bulk modulus given by:

$$ 
\begin{aligned}
\lambda_0 &= \frac{E_0 v_0}{(1 + v_0)[1 - (n - 1)v_0]}, \quad \mu_0 = \frac{E_0}{2(1 + v_0)}, \\
\kappa_0 &= \lambda_0 + \frac{2\mu_0}{n} = \frac{E_0}{n[1 - (n - 1)v_0]}.
\end{aligned}
\quad (2) $$

Given a scalar valued function: $f: x \to f(x) \in \mathbb{R}$, we define its positive part and negative part as:

$$ 
\begin{aligned}
\langle f(x) \rangle_+ &= \frac{f(x)}{2} + \frac{|f(x)|}{2}, \\
\langle f(x) \rangle_- &= \frac{f(x)}{2} - \frac{|f(x)|}{2}.
\end{aligned}
\quad (3) $$

## 2 Standard model and ideal model requirements

In this section, we first introduce the general formulation of a phase-field model for brittle fracture. We then proceed to specialize this framework to the case of the standard model. Subsequently, we clarify the limitations of this model and leverage these insights to define the ideal requirements for a model able to handle fracture under multi-axial stress states.

### 2.1 General formulation

Let us consider a homogeneous body occupying the domain $\Omega \subset \mathbb{R}^n$. Its current state at point $\mathbf{x} \in \Omega$ is described by the displacement field $\mathbf{u}(\mathbf{x})$ and the irreversible scalar damage field $\alpha(\mathbf{x}) \in [0, 1]$, with $\alpha = 0$ and $\alpha = 1$ denoting a pristine and a fully damaged material, respectively. The strain energy density is a differentiable function of the strain, the damage and the damage gradient

$$ W(\boldsymbol{\varepsilon}, \alpha, \nabla\alpha) := \varphi(\boldsymbol{\varepsilon}, \alpha) + w_1 (w(\alpha) + \ell^2 |\nabla\alpha|^2), \quad (4) $$

where the first term $\varphi$ is the elastic strain energy density and the second term is the dissipated energy density. We assume the elastic strain energy density to be positively homogeneous of degree 2 with respect to $\boldsymbol{\varepsilon}$ at fixed $\alpha$, i.e., $\varphi(r\boldsymbol{\varepsilon}, \alpha) = r^2\varphi(\boldsymbol{\varepsilon}, \alpha) \forall r \ge 0$, and to decrease with respect to $\alpha$ at fixed $\boldsymbol{\varepsilon}$. The dissipated energy density is composed of a local and a non-local part. The local term $w_1w(\alpha)$ corresponds to the amount of energy dissipated per unit volume to damage homogeneously a pristine material; the dissipation function $w(\alpha)$ is a non-negative increasing function of $\alpha$ such that $w(0) = 0$ and $w(1) = 1$, and we refer to $w_1$ as the specific fracture energy. The non-local term is assumed to be a quadratic function of the gradient of the damage, whereby $\ell$ is an internal length.

The total energy at time *t*, $E_t(\mathbf{u}, \alpha)$, is the sum of the strain energy and the potential energy of the external forces:
requirements for both nucleation and propagation. All the models that we consider, including our newly proposed energy decomposition, are of phenomenological nature. In other words, they are not grounded on microscopically motivated considerations and do not aim at reproducing fine-scale features of the fracture behavior Marigo and Kazymyrenko (2019), which may well be specific of single materials. Rather, we aim at consistently reproducing the key macroscopic features of the fracture behavior under multiaxial stress states, namely the strength surface and the transmission of stresses across fully developed cracks, that are common to a wide class of brittle materials.

The paper is structured as follows. Section 2 offers a brief review on the basic ingredients of standard phase-field modeling of brittle fracture, shows its limitations and formulates the requirements for an ideal model to describe fracture under multi-axial stress states. In light of the defined requirements, Sect. 3 reports a review of some of the available models, all based on decompositions of the elastic strain energy density. In Sect. 4, the novel star-convex model is introduced. Section 5 is dedicated to numerical experiments that showcase the advantages and limitations of the models previously illustrated, including the new one. The main conclusions are drawn in Sect. 6.

As follows, we summarize the notation and some useful relations. Vectors and second-order tensors are both denoted by boldface fonts, e.g., **u** and **σ** for the displacement vector and stress tensor, respectively. For the standard orthogonal decomposition of second-order tensors in volumetric and deviatoric parts we adopt the following notation (exemplified on **σ**)

$$ 
\begin{aligned}
\sigma &= \sigma_{\text{vol}} + \sigma_{\text{dev}}, \quad \sigma_{\text{vol}} = \frac{\text{tr}(\sigma)}{n} \mathbf{I}, \\
\sigma_{\text{dev}} &= \sigma - \frac{\text{tr}(\sigma)}{n} \mathbf{I}, \quad \sigma_{\text{vol}} \cdot \sigma_{\text{dev}} = 0,
\end{aligned}
\quad (1) $$

where **I** is the second-order identity tensor, *n* is the number of space dimensions, and the dot denotes the inner product. Note that for problems in plane-strain/stress conditions it is *n* = 3. For an isotropic elastic undamaged material with Young's modulus *E*₀ and Poisson's ratio *v*₀, we denote by (λ₀, μ₀, κ₀) the Lamé constants and the bulk modulus given by:

$$ 
\begin{aligned}
\lambda_0 &= \frac{E_0 v_0}{(1 + v_0)[1 - (n - 1)v_0]}, \quad \mu_0 = \frac{E_0}{2(1 + v_0)}, \\
\kappa_0 &= \lambda_0 + \frac{2\mu_0}{n} = \frac{E_0}{n[1 - (n - 1)v_0]}.
\end{aligned}
\quad (2) $$

Given a scalar valued function: $f: x \to f(x) \in \mathbb{R}$, we define its positive part and negative part as:

$$ 
\begin{aligned}
\langle f(x) \rangle_+ &= \frac{f(x)}{2} + \frac{|f(x)|}{2}, \\
\langle f(x) \rangle_- &= \frac{f(x)}{2} - \frac{|f(x)|}{2}.
\end{aligned}
\quad (3) $$

## 2 Standard model and ideal model requirements

In this section, we first introduce the general formulation of a phase-field model for brittle fracture. We then proceed to specialize this framework to the case of the standard model. Subsequently, we clarify the limitations of this model and leverage these insights to define the ideal requirements for a model able to handle fracture under multi-axial stress states.

### 2.1 General formulation

Let us consider a homogeneous body occupying the domain $\Omega \subset \mathbb{R}^n$. Its current state at point $\mathbf{x} \in \Omega$ is described by the displacement field $\mathbf{u}(\mathbf{x})$ and the irreversible scalar damage field $\alpha(\mathbf{x}) \in [0, 1]$, with $\alpha = 0$ and $\alpha = 1$ denoting a pristine and a fully damaged material, respectively. The strain energy density is a differentiable function of the strain, the damage and the damage gradient

$$ W(\boldsymbol{\epsilon}, \alpha, \nabla\alpha) := \varphi(\boldsymbol{\epsilon}, \alpha) + w_1 (w(\alpha) + \ell^2 |\nabla\alpha|^2), \quad (4) $$

where the first term $\varphi$ is the elastic strain energy density and the second term is the dissipated energy density. We assume the elastic strain energy density to be positively homogeneous of degree 2 with respect to $\boldsymbol{\epsilon}$ at fixed $\alpha$, i.e., $\varphi(r\boldsymbol{\epsilon}, \alpha) = r^2\varphi(\boldsymbol{\epsilon}, \alpha) \forall r \ge 0$, and to decrease with respect to $\alpha$ at fixed $\boldsymbol{\epsilon}$. The dissipated energy density is composed of a local and a non-local part. The local term $w_1w(\alpha)$ corresponds to the amount of energy dissipated per unit volume to damage homogeneously a pristine material; the dissipation function $w(\alpha)$ is a non-negative increasing function of $\alpha$ such that $w(0) = 0$ and $w(1) = 1$, and we refer to $w_1$ as the specific fracture energy. The non-local term is assumed to be a quadratic function of the gradient of the damage, whereby $\ell$ is an internal length.

The total energy at time *t*, $E_t(\mathbf{u}, \alpha)$, is the sum of the strain energy and the potential energy of the external forces:
$$ \mathcal{E}_t(\mathbf{u}, \alpha) := \int_{\Omega} W(\boldsymbol{\epsilon}(\mathbf{u}), \alpha, \nabla\alpha) d\Omega - \int_{\Omega} \mathbf{b}_t \cdot \mathbf{u} d\Omega - \int_{\partial N \Omega} f_t \cdot \mathbf{u} dS, \quad (5) $$

where $\mathbf{b}_t$ is the body force defined in $\Omega$ and $f_t$ is the surface traction applied on $\partial_N \Omega$, both at time $t$, and $\boldsymbol{\epsilon}(\mathbf{u}) := \text{sym}(\nabla \mathbf{u})$ is the linear strain tensor. The applied displacement $\bar{\mathbf{u}}_t$ is applied on the complementary part of the boundary $\partial_D \Omega$.

In the time-discrete setting of the evolution problem, given $\alpha_p$ (the damage at the previous time $t_p$), the displacement and the damage field at time $t = t_p + \Delta t$ are found by solving the minimization problem:

$$ (\mathbf{u}, \alpha) = \arg \min_{(\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p)} \mathcal{E}_t(\hat{\mathbf{u}}, \hat{\alpha}), \quad (6) $$

where

$$ C_t := \{ \mathbf{u} \in H^1(\Omega; \mathbb{R}^n) : \mathbf{u} = \bar{\mathbf{u}}_t \text{ on } \partial_D \Omega \}, \quad (7) $$

$$ \mathcal{D}(\alpha_p) := \{ \alpha \in H^1(\Omega) : \alpha \ge \alpha_p \} $$

are the spaces of the admissible displacement and damage fields at time $t$ from the previous state with damage $\alpha_p$. Equation (6) requires $(\mathbf{u}, \alpha)$ to satisfy:

$$ \forall (\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p), \quad \exists \bar{h} > 0 : \forall h \in [0, \bar{h}] \\ \mathcal{E}_t(\mathbf{u} + h(\hat{\mathbf{u}} - \mathbf{u}), \alpha + h(\hat{\alpha} - \alpha)) - \mathcal{E}_t(\mathbf{u}, \alpha) \ge 0. \quad (8) $$

A necessary condition for this local constrained minimization is found taking into account only the first-order expansion of the energy increment:

$$ \mathcal{E}'_t(\mathbf{u}, \alpha)(\hat{\mathbf{u}} - \mathbf{u}, \hat{\alpha} - \alpha) \ge 0 \quad \forall (\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p), \quad (9) $$

where

$$ \mathcal{E}'_t(\mathbf{u}, \alpha)(\mathbf{v}, \beta) := \frac{d}{dh} \mathcal{E}_t(\mathbf{u} + h\mathbf{v}, \alpha + h\beta) \Big|_{h=0} \quad (10) $$

is the *Gateaux* derivative of the functional $\mathcal{E}_t(\mathbf{u}, \alpha)$ in the direction $(\mathbf{v}, \beta)$. For smooth solutions, we can show with standard arguments of calculus of variations that the first-order optimality condition (9) is equivalent to the *equilibrium equation* and boundary conditions

$$ \operatorname{div} \sigma(\boldsymbol{\epsilon}, \alpha) + \mathbf{b}_t = \mathbf{0} \text{ in } \Omega, \quad \sigma(\boldsymbol{\epsilon}, \alpha) \mathbf{n} = f_t \text{ on } \partial_N \Omega \quad (11) $$

and to the Karush-Kuhn-Tucker (KKT) conditions and boundary conditions

$$ \begin{align} & -Y(\boldsymbol{\epsilon}, \alpha) + w_1 w'(\alpha) - 2\ell^2 w_1 \Delta \alpha \ge 0, \quad \alpha - \alpha_p \ge 0, \tag{12} \\ & (-Y(\boldsymbol{\epsilon}, \alpha) + w_1 w'(\alpha) - 2\ell^2 w_1 \Delta \alpha)(\alpha - \alpha_p) = 0 \text{ on } \Omega, \nonumber \\ & \nabla \alpha \cdot \mathbf{n} \ge 0, \quad \alpha - \alpha_p \ge 0, \quad (\nabla \alpha \cdot \mathbf{n})(\alpha - \alpha_p) = 0 \text{ on } \partial \Omega, \nonumber \end{align} $$

where $\mathbf{n}$ is the outer unit normal to the boundary. We denote the KKT conditions as *damage criterion, irreversibility and loading-unloading condition*, respectively. The conjugate quantities

$$ \boldsymbol{\sigma}(\boldsymbol{\epsilon}, \alpha) := \frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \boldsymbol{\epsilon}}, \quad Y(\boldsymbol{\epsilon}, \alpha) := -\frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \alpha} \quad (13) $$

are the *stress tensor* and the *damage energy release rate*, respectively.

A crucial notion for the following analysis of damage under multi-axial stress states is that of *elastic domains*. In the context of local damage modeling, these are the sets in which stresses and strains must remain in order for the material to follow a linearly elastic behavior without damage evolution, i.e., $\alpha = \alpha_p$. In our non-local context, their boundaries define the elastic limits for materials with homogeneous damage distribution ($\Delta \alpha = 0$). The elastic domains in the strain space $\mathcal{R}(\alpha)$ and in the stress space $\mathcal{R}^*(\alpha)$ are defined as the sets

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \text{Sym} : -\frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \alpha} \le w_1 w'(\alpha) \right\}, \quad (14) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\partial \varphi^*(\boldsymbol{\sigma}, \alpha)}{\partial \alpha} \le w_1 w'(\alpha) \right\}, \quad (15) $$

where Sym denotes the space of symmetric tensors and $\varphi^*(\boldsymbol{\sigma}, \alpha)$ is the *complementary energy density* defined as

$$ \varphi^*(\boldsymbol{\sigma}, \alpha) := \sup_{\hat{\boldsymbol{\epsilon}} \in \text{Sym}} \boldsymbol{\sigma} \cdot \hat{\boldsymbol{\epsilon}} - \varphi(\hat{\boldsymbol{\epsilon}}, \alpha). \quad (16) $$

At a given value of $\alpha \in [0, 1)$ a damage model enjoys the *strain-hardening* property if, $\forall \beta > \alpha$, $\mathcal{R}(\beta) \supset \mathcal{R}(\alpha)$ and the *stress-softening* property if $\mathcal{R}^*(\beta) \subset \mathcal{R}^*(\alpha)$, see Pham and Marigo (2010). The strain-hardening property is important to ensure the uniqueness of the solution for the damage at a given strain, while the stress-softening property is fundamental to allow for damage localization Pham and Marigo (2010).

## 2.2 The standard phase-field model

For the model that we refer to in the following as the *standard* phase-field model, the elastic strain energy density is defined as

$$ \begin{align} & \varphi(\boldsymbol{\epsilon}, \alpha) = a(\alpha)\varphi_0(\boldsymbol{\epsilon}) \quad \text{with} \nonumber \\ & \varphi_0(\boldsymbol{\epsilon}) = \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\epsilon}) + \mu_0 |\boldsymbol{\epsilon}_{\text{dev}}|^2, \quad (17) \end{align} $$
$$ \mathcal{E}_t(\mathbf{u}, \alpha) := \int_{\Omega} W(\boldsymbol{\epsilon}(\mathbf{u}), \alpha, \nabla\alpha) d\Omega - \int_{\Omega} \mathbf{b}_t \cdot \mathbf{u} d\Omega - \int_{\partial N \Omega} f_t \cdot \mathbf{u} dS, \quad (5) $$

where $\mathbf{b}_t$ is the body force defined in $\Omega$ and $f_t$ is the surface traction applied on $\partial_N \Omega$, both at time $t$, and $\boldsymbol{\epsilon}(\mathbf{u}) := \text{sym}(\nabla \mathbf{u})$ is the linear strain tensor. The applied displacement $\bar{\mathbf{u}}_t$ is applied on the complementary part of the boundary $\partial_D \Omega$.

In the time-discrete setting of the evolution problem, given $\alpha_p$ (the damage at the previous time $t_p$), the displacement and the damage field at time $t = t_p + \Delta t$ are found by solving the minimization problem:

$$ (\mathbf{u}, \alpha) = \arg\min_{(\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p)} \mathcal{E}_t(\hat{\mathbf{u}}, \hat{\alpha}), \quad (6) $$

where

$$ \begin{aligned} C_t &:= \{ \mathbf{u} \in H^1(\Omega; \mathbb{R}^n) : \mathbf{u} = \bar{\mathbf{u}}_t \text{ on } \partial_D \Omega \}, \\ \mathcal{D}(\alpha_p) &:= \{ \alpha \in H^1(\Omega) : \alpha \ge \alpha_p \} \end{aligned} \quad (7) $$

are the spaces of the admissible displacement and damage fields at time $t$ from the previous state with damage $\alpha_p$. Equation (6) requires $(\mathbf{u}, \alpha)$ to satisfy:

$$ \begin{aligned} \forall (\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p), \quad \exists \bar{h} > 0 : \forall h \in [0, \bar{h}] \\ \mathcal{E}_t(\mathbf{u} + h(\hat{\mathbf{u}} - \mathbf{u}), \alpha + h(\hat{\alpha} - \alpha)) - \mathcal{E}_t(\mathbf{u}, \alpha) \ge 0. \end{aligned} \quad (8) $$

A necessary condition for this local constrained minimization is found taking into account only the first-order expansion of the energy increment:

$$ \mathcal{E}'_t(\mathbf{u}, \alpha)(\hat{\mathbf{u}} - \mathbf{u}, \hat{\alpha} - \alpha) \ge 0 \quad \forall (\hat{\mathbf{u}}, \hat{\alpha}) \in C_t \times \mathcal{D}(\alpha_p), \quad (9) $$

where

$$ \mathcal{E}'_t(\mathbf{u}, \alpha)(\mathbf{v}, \beta) := \frac{d}{dh} \mathcal{E}_t(\mathbf{u} + h\mathbf{v}, \alpha + h\beta) \Big|_{h=0} \quad (10) $$

is the *Gateaux* derivative of the functional $\mathcal{E}_t(\mathbf{u}, \alpha)$ in the direction $(\mathbf{v}, \beta)$. For smooth solutions, we can show with standard arguments of calculus of variations that the first-order optimality condition (9) is equivalent to the *equilibrium equation* and boundary conditions

$$ \operatorname{div}\sigma(\boldsymbol{\epsilon}, \alpha) + \mathbf{b}_t = \mathbf{0} \text{ in } \Omega, \quad \sigma(\boldsymbol{\epsilon}, \alpha) \mathbf{n} = f_t \text{ on } \partial_N \Omega \quad (11) $$

and to the Karush-Kuhn-Tucker (KKT) conditions and boundary conditions

$$ \begin{aligned} & -Y(\boldsymbol{\epsilon}, \alpha) + w_1 w'(\alpha) - 2\ell^2 w_1 \Delta \alpha \ge 0, \quad \alpha - \alpha_p \ge 0, \\ & (-Y(\boldsymbol{\epsilon}, \alpha) + w_1 w'(\alpha) - 2\ell^2 w_1 \Delta \alpha)(\alpha - \alpha_p) = 0 \text{ on } \Omega, \quad (12) \\ & \nabla \alpha \cdot \mathbf{n} \ge 0, \quad \alpha - \alpha_p \ge 0, \quad (\nabla \alpha \cdot \mathbf{n})(\alpha - \alpha_p) = 0 \text{ on } \partial \Omega, \end{aligned} $$

where $\mathbf{n}$ is the outer unit normal to the boundary. We denote the KKT conditions as *damage criterion, irreversibility and loading-unloading condition*, respectively. The conjugate quantities

$$ \boldsymbol{\sigma}(\boldsymbol{\epsilon}, \alpha) := \frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \boldsymbol{\epsilon}}, \quad Y(\boldsymbol{\epsilon}, \alpha) := -\frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \alpha} \quad (13) $$

are the *stress tensor* and the *damage energy release rate*, respectively.

A crucial notion for the following analysis of damage under multi-axial stress states is that of *elastic domains*. In the context of local damage modeling, these are the sets in which stresses and strains must remain in order for the material to follow a linearly elastic behavior without damage evolution, i.e., $\alpha = \alpha_p$. In our non-local context, their boundaries define the elastic limits for materials with homogeneous damage distribution ($\Delta \alpha = 0$). The elastic domains in the strain space $\mathcal{R}(\alpha)$ and in the stress space $\mathcal{R}^*(\alpha)$ are defined as the sets

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \text{Sym} : -\frac{\partial \varphi(\boldsymbol{\epsilon}, \alpha)}{\partial \alpha} \le w_1 w'(\alpha) \right\}, \quad (14) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\partial \varphi^*(\boldsymbol{\sigma}, \alpha)}{\partial \alpha} \le w_1 w'(\alpha) \right\}, \quad (15) $$

where Sym denotes the space of symmetric tensors and $\varphi^*(\boldsymbol{\sigma}, \alpha)$ is the *complementary energy density* defined as

$$ \varphi^*(\boldsymbol{\sigma}, \alpha) := \sup_{\hat{\boldsymbol{\epsilon}} \in \text{Sym}} \boldsymbol{\sigma} \cdot \hat{\boldsymbol{\epsilon}} - \varphi(\hat{\boldsymbol{\epsilon}}, \alpha). \quad (16) $$

At a given value of $\alpha \in [0, 1)$ a damage model enjoys the *strain-hardening* property if, $\forall \beta > \alpha$, $\mathcal{R}(\beta) \supset \mathcal{R}(\alpha)$ and the *stress-softening* property if $\mathcal{R}^*(\beta) \subset \mathcal{R}^*(\alpha)$, see Pham and Marigo (2010). The strain-hardening property is important to ensure the uniqueness of the solution for the damage at a given strain, while the stress-softening property is fundamental to allow for damage localization Pham and Marigo (2010).

## 2.2 The standard phase-field model

For the model that we refer to in the following as the *standard* phase-field model, the elastic strain energy density is defined as

$$ \begin{aligned} & \varphi(\boldsymbol{\epsilon}, \alpha) = a(\alpha)\varphi_0(\boldsymbol{\epsilon}) \quad \text{with} \\ & \varphi_0(\boldsymbol{\epsilon}) = \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\epsilon}) + \mu_0 |\boldsymbol{\epsilon}_{\text{dev}}|^2, \end{aligned} \quad (17) $$
where $\varphi_0$ is the elastic strain energy density for a homogeneous isotropic linear elastic material and $a(\alpha)$ is the degradation function. This function describes the degradation of the linear elastic properties with damage; it is a decreasing function of $\alpha$ going from $a(0) = 1$ to $a(1) = 0$. By Legendre transformation, the complementary energy density $\varphi^*(\boldsymbol{\sigma}, \alpha)$ is derived as

$$ \varphi^*(\boldsymbol{\sigma}, \alpha) = s(\alpha)\varphi_0^*(\boldsymbol{\sigma}) \quad \text{with} \\ \varphi_0^*(\boldsymbol{\sigma}) = \frac{\text{tr}^2(\boldsymbol{\sigma})}{2n^2\kappa_0} + \frac{|\boldsymbol{\sigma}_{\text{dev}}|^2}{4\mu_0} \quad \text{and} \quad s(\alpha) := \frac{1}{a(\alpha)}. \quad (18) $$

Combining (14-18) we obtain

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \text{Sym} : \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\epsilon}) + \mu_0 |\boldsymbol{\epsilon}_{\text{dev}}|^2 \leq - \frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \quad (19) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\text{tr}^2(\boldsymbol{\sigma})}{2n^2\kappa_0} + \frac{|\boldsymbol{\sigma}_{\text{dev}}|^2}{4\mu_0} \le \frac{w_1 w'(\alpha)}{s'(\alpha)} \right\}. \quad (20) $$

The expressions of the elastic strain energy density in (17) and of the complementary energy density in (18) ensure that the transformation of spaces $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with varying $\alpha$ is a homothety centered in $\boldsymbol{\epsilon} = \mathbf{0}$ and $\boldsymbol{\sigma} = \mathbf{0}$, respectively.

The selection of the degradation function $a(\alpha)$ and of the dissipation function $w(\alpha)$ has a significant impact on the description of the damaging behavior. Two classical expressions are

$$ \text{AT}_1: \quad a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha, \quad \alpha_{\text{peak}} = 0, \quad (21) $$

$$ \text{AT}_2: \quad a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha^2, \quad \alpha_{\text{peak}} = \frac{1}{4}. \quad (22) $$

These two models are both strain-hardening $\forall \alpha \in [0, 1]$, and stress-softening $\forall \alpha \in [\alpha_{\text{peak}}, 1]$, whereas they are stress-hardening $\forall \alpha \in [0, \alpha_{\text{peak}})$. These properties hold regardless of the loading direction due to the homothetic evolution of $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with $\alpha$. The primary difference between AT$_1$ and AT$_2$ is that AT$_1$ displays a linearly elastic behavior up to a non-zero elastic limit stress, whereas AT$_2$ features a zero elastic limit stress, therefore with AT$_2$ an infinitesimal load is sufficient to trigger the onset of damage.

For AT$_1$ and AT$_2$ applied to the 1D bar under tensile loading it is shown in (Pham et al. 2011a, b) that, for sufficiently long bars ($L \gg \ell$), the damage localization in bands (crack nucleation) occurs at the level of damage at which the behavior of the model changes from stress hardening to stress softening, i.e., at $\alpha = \alpha_{\text{peak}}$. For AT$_1$ this transition corresponds to the elastic limit because

$\alpha_{\text{peak}} = 0$. Under multi-axial loading, in (Pham and Marigo 2013) it is also shown that with the AT1 model and for a sufficiently large structure the transition from a stress-hardening to a stress-softening phase is a necessary and sufficient condition for the damage localization in bands (or cracks). For different types of phase-field models, these conclusions do not necessarily hold and a more careful study would be necessary, see (Pham and Marigo 2013; Zolesi and Maurini 2023) for more details.

From now on, we focus only on AT$_1$-like models ($w(\alpha) = \alpha$). Accordingly, we define the strength surface $\mathcal{S}^*$ as the set of stresses at the elastic limit, i.e., all stresses at the boundary of the elastic domain $\mathcal{R}^*(\alpha)$ for $\alpha = \alpha_{\text{peak}} = 0$:

$$ \mathcal{S}^* := \partial \mathcal{R}^*(0) = \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\text{tr}^2(\boldsymbol{\sigma})}{2n^2\kappa_0} + \frac{|\boldsymbol{\sigma}_{\text{dev}}|^2}{4\mu_0} = \frac{w_1 w'(0)}{s'(0)} \right\}. \quad (23) $$

From $\mathcal{S}^*$, one can define the tensile and compressive strengths, $\sigma_e^+$ and $\sigma_e^-$, as the maximum and minimal allowable stress $\sigma$ for the uniaxial stress state $\boldsymbol{\sigma} = \boldsymbol{\sigma} e_1 \otimes e_1$, and the shear strength $\tau_e$ as the maximum allowable stress $\tau$ for stress states of pure shear $\boldsymbol{\sigma} = \tau e_1 \otimes e_2$. For the standard model, these quantities are given by

$$ \sigma_e^+ := \sqrt{\frac{2E_0 w_1 w'(0)}{s'(0)}}, \quad \tau_e := \sqrt{\frac{2\mu_0 w_1 w'(0)}{s'(0)}}, \\ \sigma_e^- := -\sqrt{\frac{2E_0 w_1 w'(0)}{s'(0)}}. \quad (24) $$

In the setting of a 1D bar under tensile loading, the dissipated energy associated to the phase-field smeared representation of a crack is regarded as the dissipation upon rupture and is denoted as fracture toughness $G_c$. In (Pham et al. 2011a, b), the 1D localized damage distribution is derived analytically and the fracture toughness is accordingly expressed as

$$ G_c = 4w_1\ell \int_0^1 \sqrt{w(\beta)}d\beta, \quad (25) $$

which for AT$_1$ gives

$$ G_c = \frac{8}{3}w_1\ell. \quad (26) $$

Thus, from the experimental determination of $G_c$ and $\sigma_e^+$, it is possible to calibrate $w_1$ and $\ell$. In this sense, the regularization length can be viewed as a material property. In Fig. 1, we plot the evolution of the
where $\varphi_0$ is the elastic strain energy density for a homogeneous isotropic linear elastic material and $a(\alpha)$ is the degradation function. This function describes the degradation of the linear elastic properties with damage; it is a decreasing function of $\alpha$ going from $a(0) = 1$ to $a(1) = 0$. By Legendre transformation, the complementary energy density $\varphi^*(\sigma, \alpha)$ is derived as

$$ \varphi^*(\sigma, \alpha) = s(\alpha)\varphi_0^*(\sigma) \quad \text{with} \\ \varphi_0^*(\sigma) = \frac{\text{tr}^2(\sigma)}{2n^2\kappa_0} + \frac{|\sigma_{\text{dev}}|^2}{4\mu_0} \quad \text{and} \quad s(\alpha) := \frac{1}{a(\alpha)}. \quad (18) $$

Combining (14-18) we obtain

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \text{Sym} : \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\epsilon}) + \mu_0 |\boldsymbol{\epsilon}_{\text{dev}}|^2 \leq - \frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \quad (19) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\text{tr}^2(\boldsymbol{\sigma})}{2n^2\kappa_0} + \frac{|\boldsymbol{\sigma}_{\text{dev}}|^2}{4\mu_0} \leq \frac{w_1 w'(\alpha)}{s'(\alpha)} \right\}. \quad (20) $$

The expressions of the elastic strain energy density in (17) and of the complementary energy density in (18) ensure that the transformation of spaces $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with varying $\alpha$ is a homothety centered in $\boldsymbol{\epsilon} = \mathbf{0}$ and $\boldsymbol{\sigma} = \mathbf{0}$, respectively.

The selection of the degradation function $a(\alpha)$ and of the dissipation function $w(\alpha)$ has a significant impact on the description of the damaging behavior. Two classical expressions are

$$ \text{AT}_1: \quad a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha, \quad \alpha_{\text{peak}} = 0, \quad (21) $$

$$ \text{AT}_2: \quad a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha^2, \quad \alpha_{\text{peak}} = \frac{1}{4}. \quad (22) $$

These two models are both strain-hardening $\forall \alpha \in [0, 1]$, and stress-softening $\forall \alpha \in [\alpha_{\text{peak}}, 1]$, whereas they are stress-hardening $\forall \alpha \in [0, \alpha_{\text{peak}})$. These properties hold regardless of the loading direction due to the homothetic evolution of $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with $\alpha$. The primary difference between AT$_1$ and AT$_2$ is that AT$_1$ displays a linearly elastic behavior up to a non-zero elastic limit stress, whereas AT$_2$ features a zero elastic limit stress, therefore with AT$_2$ an infinitesimal load is sufficient to trigger the onset of damage.

For AT$_1$ and AT$_2$ applied to the 1D bar under tensile loading it is shown in (Pham et al. 2011a, b) that, for sufficiently long bars ($L \gg \ell$), the damage localization in bands (crack nucleation) occurs at the level of damage at which the behavior of the model changes from stress hardening to stress softening, i.e., at $\alpha = \alpha_{\text{peak}}$. For AT$_1$ this transition corresponds to the elastic limit because

$\alpha_{\text{peak}} = 0$. Under multi-axial loading, in (Pham and Marigo 2013) it is also shown that with the AT1 model and for a sufficiently large structure the transition from a stress-hardening to a stress-softening phase is a necessary and sufficient condition for the damage localization in bands (or cracks). For different types of phase-field models, these conclusions do not necessarily hold and a more careful study would be necessary, see (Pham and Marigo 2013; Zolesi and Maurini 2023) for more details.

From now on, we focus only on AT$_1$-like models ($w(\alpha) = \alpha$). Accordingly, we define the strength surface $\mathcal{S}^*$ as the set of stresses at the elastic limit, i.e., all stresses at the boundary of the elastic domain $\mathcal{R}^*(\alpha)$ for $\alpha = \alpha_{\text{peak}} = 0$:

$$ \mathcal{S}^* := \partial \mathcal{R}^*(0) = \left\{ \boldsymbol{\sigma} \in \text{Sym} : \frac{\text{tr}^2(\boldsymbol{\sigma})}{2n^2\kappa_0} + \frac{|\boldsymbol{\sigma}_{\text{dev}}|^2}{4\mu_0} = \frac{w_1 w'(0)}{s'(0)} \right\}. \quad (23) $$

From $\mathcal{S}^*$, one can define the tensile and compressive strengths, $\sigma_e^+$ and $\sigma_e^-$, as the maximum and minimal allowable stress $\sigma$ for the uniaxial stress state $\boldsymbol{\sigma} = \boldsymbol{\sigma} e_1 \otimes e_1$, and the shear strength $\tau_e$ as the maximum allowable stress $\tau$ for stress states of pure shear $\boldsymbol{\sigma} = \tau e_1 \otimes e_2$. For the standard model, these quantities are given by

$$ \sigma_e^+ := \sqrt{\frac{2E_0 w_1 w'(0)}{s'(0)}}, \quad \tau_e := \sqrt{\frac{2\mu_0 w_1 w'(0)}{s'(0)}}, \\ \sigma_e^- := -\sqrt{\frac{2E_0 w_1 w'(0)}{s'(0)}}. \quad (24) $$

In the setting of a 1D bar under tensile loading, the dissipated energy associated to the phase-field smeared representation of a crack is regarded as the dissipation upon rupture and is denoted as fracture toughness $G_c$. In (Pham et al. 2011a, b), the 1D localized damage distribution is derived analytically and the fracture toughness is accordingly expressed as

$$ G_c = 4w_1\ell \int_0^1 \sqrt{w(\beta)} d\beta, \quad (25) $$

which for AT$_1$ gives

$$ G_c = \frac{8}{3} w_1 \ell. \quad (26) $$

Thus, from the experimental determination of $G_c$ and $\sigma_e^+$, it is possible to calibrate $w_1$ and $\ell$. In this sense, the regularization length can be viewed as a material property. In Fig. 1, we plot the evolution of the
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Elastic domains for the standard model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$</figcaption></figure>

elastic domains $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with increasing $\alpha$ in both the volumetric-deviatoric plane (Fig. 1a) and the plane spanned by the principal strain/stress components under the plane-strain assumption (Fig. 1b). Given the symmetry of $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with respect to the volumetric axis in Fig. 1a and with respect to the bisector of the first and third quadrants in Fig. 1b, we plot only half of $\mathcal{R}(\alpha)$ (in blue) and half of $\mathcal{R}^*(\alpha)$ (in red). Note that for the elastic domains $\mathcal{R}^*(\alpha)$ in Fig. 1b, depicted in the plane of the principal components $\sigma_{xx} - \sigma_{yy}$, the value of $\sigma_{zz}$ is not constant but related to $\sigma_{xx}$ and $\sigma_{yy}$ through the constitutive law and the plane-strain constraint. The same is true for all the other figures depicting the elastic domains $\mathcal{R}^*(\alpha)$ in the principal stress components.

## 2.3 Limitations of the standard model

The standard phase-field model has major limitations when considering multi-axial stress states. Regarding nucleation, we mention two limits evident from the shape of the elastic domains at $\alpha = 0$ (Figure 1):

* **Tension/compression symmetry:** Both $\mathcal{R}(0)$ and $\mathcal{R}^*(0)$ are ellipsoids in the strain and the stress

(thick lines) are the nucleation domains (or strength surfaces). Volumetric-deviatoric diagram **a**. Principal components diagram under the plane-strain assumption **b**

spaces, respectively. As a consequence, the tensile and compressive strengths have the same magnitude. This is unrealistic for the majority of brittle materials, which typically feature a compressive strength one order of magnitude higher than the tensile strength. Moreover, the shear strength often increases with volumetric compression.

* **Lack of flexibility:** Assuming the elastic constants $\kappa_0$, $\mu_0$, the tensile strength $\sigma_e^+$ and the fracture toughness $G_c$ to be known, it is possible to calibrate $w_1$ and $\ell$ using (24) and (25). In this manner, all the model parameters are determined, thus the shape and size of the ellipsoids $\mathcal{R}(0)$ and $\mathcal{R}^*(0)$ are also fixed and it is not possible to match e.g. experimentally known values of $\sigma_e^-$ and/or $\tau_e$. This aspect is what we term lack of flexibility of the model.

Regarding propagation, we report a major limit related to the local behavior for $\alpha = 1$:

* **Zero residual stress:** In the standard model, the whole elastic strain energy density drives the damage evolution. The consequence is that, once the maximum damage $\alpha = 1$ is locally reached, the
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Elastic domains for the standard model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$</figcaption></figure>

elastic domains $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with increasing $\alpha$ in both the volumetric-deviatoric plane (Fig. 1a) and the plane spanned by the principal strain/stress components under the plane-strain assumption (Fig. 1b). Given the symmetry of $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ with respect to the volumetric axis in Fig. 1a and with respect to the bisector of the first and third quadrants in Fig. 1b, we plot only half of $\mathcal{R}(\alpha)$ (in blue) and half of $\mathcal{R}^*(\alpha)$ (in red). Note that for the elastic domains $\mathcal{R}^*(\alpha)$ in Fig. 1b, depicted in the plane of the principal components $\sigma_{xx} - \sigma_{yy}$, the value of $\sigma_{zz}$ is not constant but related to $\sigma_{xx}$ and $\sigma_{yy}$ through the constitutive law and the plane-strain constraint. The same is true for all the other figures depicting the elastic domains $\mathcal{R}^*(\alpha)$ in the principal stress components.

## 2.3 Limitations of the standard model

The standard phase-field model has major limitations when considering multi-axial stress states. Regarding nucleation, we mention two limits evident from the shape of the elastic domains at $\alpha = 0$ (Figure 1):

* **Tension/compression symmetry:** Both $\mathcal{R}(0)$ and $\mathcal{R}^*(0)$ are ellipsoids in the strain and the stress

(thick lines) are the nucleation domains (or strength surfaces). Volumetric-deviatoric diagram **a**. Principal components diagram under the plane-strain assumption **b**

spaces, respectively. As a consequence, the tensile and compressive strengths have the same magnitude. This is unrealistic for the majority of brittle materials, which typically feature a compressive strength one order of magnitude higher than the tensile strength. Moreover, the shear strength often increases with volumetric compression.

* **Lack of flexibility:** Assuming the elastic constants $\kappa_0, \mu_0$, the tensile strength $\sigma_e^+$ and the fracture toughness $G_c$ to be known, it is possible to calibrate $w_1$ and $\ell$ using (24) and (25). In this manner, all the model parameters are determined, thus the shape and size of the ellipsoids $\mathcal{R}(0)$ and $\mathcal{R}^*(0)$ are also fixed and it is not possible to match e.g. experimentally known values of $\sigma_e^-$ and/or $\tau_e$. This aspect is what we term lack of flexibility of the model.

Regarding propagation, we report a major limit related to the local behavior for $\alpha = 1$:

* **Zero residual stress:** In the standard model, the whole elastic strain energy density drives the damage evolution. The consequence is that, once the maximum damage $\alpha = 1$ is locally reached, the
elastic strain energy density vanishes regardless of the loading direction and the elastic domain in the stress space collapses to the point $\sigma = 0$ (Figure 1). From now on, we denote the stress obtained for $\alpha = 1$, i.e.,

$$ \sigma_R(\boldsymbol{\epsilon}) := \sigma(\boldsymbol{\epsilon}, 1) = \frac{\partial \varphi(\boldsymbol{\epsilon}, 1)}{\partial \boldsymbol{\epsilon}}, \qquad (27) $$

as residual stress. Thus, for the standard model it is $\sigma_R = 0$. This correctly avoids the transmission of tensile tractions across the crack boundary in mode 1 (opening), mode 2 and 3 (shear). However, it also gives zero tractions under compressive loading, thus it cannot properly represent unilateral contact at the crack faces. Determining the threshold that defines which load directions should be associated with zero residual stress and which should not to best model unilateral contact is a nontrivial task. In Chambolle et al. (2018), it is shown that a zero residual stress for $\mathrm{tr}(\boldsymbol{\epsilon}) \ge 0$ and a residual elastic energy density which depends solely on the volumetric energy density associated to the negative part of the strain trace (e.g., $\varphi(\boldsymbol{\epsilon}, 1) = \frac{1}{2}\kappa_0 \langle\mathrm{tr}(\boldsymbol{\epsilon})\rangle^2$) can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence.

## 2.4 Requirements for an ideal model

In light of the previous observations on the standard model, we can now define the requirements for an ideal model as follows:

* **Strain-hardening:** This ensures the uniqueness of the solution for the damage at a given strain.

* **Stress-softening:** It is fundamental to allow for the presence of solutions with localised damage and hence for crack nucleation.

* **Tension/compression asymmetry:** At least we require that $|\sigma_e^-| > \sigma_e^+$. For some materials it would be ideal to have a shear strength that increases with volumetric compression.

* **Flexibility:** The model should contain enough parameters to allow for calibration of $\tau_e$ and $\sigma_e^-$ (or at least one of the two) independently of $\sigma_e^+$ (for given elastic properties and $G_c$).

* **Crack-like residual stress:** The ideal model does not transmit stress through the crack faces in the opening and shear modes but exhibits a compres-

sive residual stress. Hinging on the consistent variational framework in (Chambolle et al. 2018), we require the presence of nonzero residual stresses exclusively for $\mathrm{tr}(\boldsymbol{\epsilon}) < 0$. Accordingly, an elastic energy density such that $\varphi(\boldsymbol{\epsilon}, 1) \propto \langle\mathrm{tr}(\boldsymbol{\epsilon})\rangle^2$ fulfills the requirement.

Several models have been proposed in the literature to solve the limitations of the standard model under multi-axial stress states. In the following, we analyze some of these models in light of the requirements listed above.

## 3 Available variational phase-field models for multi-axial stress states

In this section, we review the main available variational phase-field models for multi-axial stress states. The advantages and limitations of such models in relation to the requirements listed in Sect. 2.4 are then exemplified through numerical experiments in Sect. 5. The models are all based on the following decomposition of the elastic strain energy (also denoted as energy split)

$$ \varphi(\boldsymbol{\epsilon}, \alpha) = a(\alpha)\varphi_D(\boldsymbol{\epsilon}) + \varphi_R(\boldsymbol{\epsilon}) \quad \text{with} \\ \varphi_D(\boldsymbol{\epsilon}) + \varphi_R(\boldsymbol{\epsilon}) = \varphi_0(\boldsymbol{\epsilon}), \qquad (28) $$

where $\varphi_D(\boldsymbol{\epsilon})$ and $\varphi_R(\boldsymbol{\epsilon})$ are respectively the *degradable* and the *residual* components of the elastic strain energy density. These components are non-negative and only vanish for $\boldsymbol{\epsilon} = 0$; their sum yields the elastic energy density of a pristine material $\varphi_0(\boldsymbol{\epsilon})$, defined as in (17). The idea behind this decomposition is that only certain modes of deformation contribute to the driving force for the nucleation and evolution of damage. Residual stresses at full damage are given by

$$ \sigma_R(\boldsymbol{\epsilon}) = \sigma(\boldsymbol{\epsilon}, 1) = \frac{\partial \varphi_R(\boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}}. \qquad (29) $$

Only $\varphi_D(\boldsymbol{\epsilon})$ contributes to the energy release rate $Y(\boldsymbol{\epsilon}, \alpha)$; accordingly, it is the only component which affects the elastic domains

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \mathrm{Sym} : \varphi_D(\boldsymbol{\epsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \qquad (30) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \mathrm{Sym} : \varphi_D(\boldsymbol{\epsilon}(\boldsymbol{\sigma}, \alpha)) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \qquad (31) $$

where $\boldsymbol{\epsilon}(\boldsymbol{\sigma}, \alpha)$ is computed as the inverse of the constitutive law in (13). The choice of the decomposition strongly affects the evolution of the elastic domain and
elastic strain energy density vanishes regardless of the loading direction and the elastic domain in the stress space collapses to the point $\sigma = 0$ (Figure 1). From now on, we denote the stress obtained for $\alpha = 1$, i.e.,

$$ \sigma_R(\boldsymbol{\epsilon}) := \sigma(\boldsymbol{\epsilon}, 1) = \frac{\partial \varphi(\boldsymbol{\epsilon}, 1)}{\partial \boldsymbol{\epsilon}}, \qquad (27) $$

as residual stress. Thus, for the standard model it is $\sigma_R = 0$. This correctly avoids the transmission of tensile tractions across the crack boundary in mode 1 (opening), mode 2 and 3 (shear). However, it also gives zero tractions under compressive loading, thus it cannot properly represent unilateral contact at the crack faces. Determining the threshold that defines which load directions should be associated with zero residual stress and which should not to best model unilateral contact is a nontrivial task. In Chambolle et al. (2018), it is shown that a zero residual stress for $\mathrm{tr}(\boldsymbol{\epsilon}) \ge 0$ and a residual elastic energy density which depends solely on the volumetric energy density associated to the negative part of the strain trace (e.g., $\varphi(\boldsymbol{\epsilon}, 1) = \frac{1}{2} \kappa_0 \langle\mathrm{tr}(\boldsymbol{\epsilon})\rangle^2$) can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence.

## 2.4 Requirements for an ideal model

In light of the previous observations on the standard model, we can now define the requirements for an ideal model as follows:

* **Strain-hardening:** This ensures the uniqueness of the solution for the damage at a given strain.

* **Stress-softening:** It is fundamental to allow for the presence of solutions with localised damage and hence for crack nucleation.

* **Tension/compression asymmetry:** At least we require that $|\sigma_e^-| > \sigma_e^+$. For some materials it would be ideal to have a shear strength that increases with volumetric compression.

* **Flexibility:** The model should contain enough parameters to allow for calibration of $\tau_e$ and $\sigma_e^-$ (or at least one of the two) independently of $\sigma_e^+$ (for given elastic properties and $G_c$).

* **Crack-like residual stress:** The ideal model does not transmit stress through the crack faces in the opening and shear modes but exhibits a compres-

sive residual stress. Hinging on the consistent variational framework in (Chambolle et al. 2018), we require the presence of nonzero residual stresses exclusively for $\mathrm{tr}(\boldsymbol{\epsilon}) < 0$. Accordingly, an elastic energy density such that $\varphi(\boldsymbol{\epsilon}, 1) \propto \langle\mathrm{tr}(\boldsymbol{\epsilon})\rangle^2$ fulfills the requirement.

Several models have been proposed in the literature to solve the limitations of the standard model under multi-axial stress states. In the following, we analyze some of these models in light of the requirements listed above.

## 3 Available variational phase-field models for multi-axial stress states

In this section, we review the main available variational phase-field models for multi-axial stress states. The advantages and limitations of such models in relation to the requirements listed in Sect. 2.4 are then exemplified through numerical experiments in Sect. 5. The models are all based on the following decomposition of the elastic strain energy (also denoted as energy split)

$$ \varphi(\boldsymbol{\epsilon}, \alpha) = a(\alpha)\varphi_D(\boldsymbol{\epsilon}) + \varphi_R(\boldsymbol{\epsilon}) \quad \text{with} \\ \varphi_D(\boldsymbol{\epsilon}) + \varphi_R(\boldsymbol{\epsilon}) = \varphi_0(\boldsymbol{\epsilon}), \qquad (28) $$

where $\varphi_D(\boldsymbol{\epsilon})$ and $\varphi_R(\boldsymbol{\epsilon})$ are respectively the *degradable* and the *residual* components of the elastic strain energy density. These components are non-negative and only vanish for $\boldsymbol{\epsilon} = 0$; their sum yields the elastic energy density of a pristine material $\varphi_0(\boldsymbol{\epsilon})$, defined as in (17). The idea behind this decomposition is that only certain modes of deformation contribute to the driving force for the nucleation and evolution of damage. Residual stresses at full damage are given by

$$ \sigma_R(\boldsymbol{\epsilon}) = \sigma(\boldsymbol{\epsilon}, 1) = \frac{\partial \varphi_R(\boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}}. \qquad (29) $$

Only $\varphi_D(\boldsymbol{\epsilon})$ contributes to the energy release rate $Y(\boldsymbol{\epsilon}, \alpha)$; accordingly, it is the only component which affects the elastic domains

$$ \mathcal{R}(\alpha) := \left\{ \boldsymbol{\epsilon} \in \mathrm{Sym} : \varphi_D(\boldsymbol{\epsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \qquad (30) $$

$$ \mathcal{R}^*(\alpha) := \left\{ \boldsymbol{\sigma} \in \mathrm{Sym} : \varphi_D(\boldsymbol{\epsilon}(\boldsymbol{\sigma}, \alpha)) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)} \right\}, \qquad (31) $$

where $\boldsymbol{\epsilon}(\boldsymbol{\sigma}, \alpha)$ is computed as the inverse of the constitutive law in (13). The choice of the decomposition strongly affects the evolution of the elastic domain and
<figure><img src="image_3.png" /><figcaption>Fig. 2 Elastic domains for the model with volumetric-deviatoric split in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The</figcaption></figure>

grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram a. Principal components diagram under the plane-strain assumption b.

thus the predicted crack nucleation and propagation behavior under multi-axial stress states. Note that formulating the decomposition in the form of (28) ensures that $\mathcal{R}(\alpha)$ evolves with $\alpha$ as a simple homothety centered in the origin. Unlike in the case of the standard model, this is not always guaranteed for $\mathcal{R}^*(\alpha)$.

### 3.1 The volumetric-deviatoric split

The split proposed by Amor et al. (Amor et al. 2009) is based on the decomposition of the elastic strain energy density into a deviatoric and a volumetric part:

$$ 
\begin{aligned}
\varphi_D(\boldsymbol{\epsilon}) &= \frac{1}{2} \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_+^2 + \mu_0 |\boldsymbol{\epsilon}_{\mathrm{dev}}|^2, \\
\varphi_R(\boldsymbol{\epsilon}) &= \frac{1}{2} \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_-^2
\end{aligned} 
\qquad (32) $$

leading to the residual stresses

$$ \boldsymbol{\sigma}_R(\boldsymbol{\epsilon}) = \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_- \mathbf{I}. \qquad (33) $$

This model was constructed to recover unilateral contact conditions under compression. Indeed, in (Cham-

bolle et al. 2018) it is demonstrated that the volumetric-deviatoric decomposition can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence without affecting the tensile and shear behavior in the presence of a crack. The elastic domains are given in Appendix A.1 and illustrated in Fig. 2. In this case, both domains evolve with $\alpha$ as homotheties centered in the origin.

This split introduces the desired asymmetry in tension and compression, as shown in Table 1, without adding residual stresses for $\mathrm{tr}(\boldsymbol{\epsilon}) > 0$. However, it gives no flexibility in the choice of the shear strength, which is the same as in the standard model, nor in the compressive strength.

### 3.2 The spectral split

The decomposition introduced by Miehe et al. (Miehe et al. 2010) is based on the eigenvalues and eigenvectors of the strain tensor
<figure><img src="image_3.png" /><figcaption>Fig. 2 Elastic domains for the model with volumetric-deviatoric split in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The</figcaption></figure>

grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram a. Principal components diagram under the plane-strain assumption b.

thus the predicted crack nucleation and propagation behavior under multi-axial stress states. Note that formulating the decomposition in the form of (28) ensures that $\mathcal{R}(\alpha)$ evolves with $\alpha$ as a simple homothety centered in the origin. Unlike in the case of the standard model, this is not always guaranteed for $\mathcal{R}^*(\alpha)$.

### 3.1 The volumetric-deviatoric split

The split proposed by Amor et al. (Amor et al. 2009) is based on the decomposition of the elastic strain energy density into a deviatoric and a volumetric part:

$$ 
\begin{aligned}
\varphi_D(\boldsymbol{\epsilon}) &= \frac{1}{2} \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_+^2 + \mu_0 |\boldsymbol{\epsilon}_{\mathrm{dev}}|^2, \\
\varphi_R(\boldsymbol{\epsilon}) &= \frac{1}{2} \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_-^2
\end{aligned} 
\qquad (32) 
$$

leading to the residual stresses

$$ \boldsymbol{\sigma}_R(\boldsymbol{\epsilon}) = \kappa_0 \langle \mathrm{tr}(\boldsymbol{\epsilon}) \rangle_- \mathbf{I}. \qquad (33) $$

This model was constructed to recover unilateral contact conditions under compression. Indeed, in (Cham-

bolle et al. 2018) it is demonstrated that the volumetric-deviatoric decomposition can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence without affecting the tensile and shear behavior in the presence of a crack. The elastic domains are given in Appendix A.1 and illustrated in Fig. 2. In this case, both domains evolve with $\alpha$ as homotheties centered in the origin.

This split introduces the desired asymmetry in tension and compression, as shown in Table 1, without adding residual stresses for $\mathrm{tr}(\boldsymbol{\epsilon}) > 0$. However, it gives no flexibility in the choice of the shear strength, which is the same as in the standard model, nor in the compressive strength.

### 3.2 The spectral split

The decomposition introduced by Miehe et al. (Miehe et al. 2010) is based on the eigenvalues and eigenvectors of the strain tensor
<table><caption>**Table 1** Tensile, compressive and shear strengths for all the analyzed models</caption><thead><tr><th>Decomposition</th><th>σ<sub>e</sub><sup>+</sup></th><th>σ<sub>e</sub><sup>-</sup></th><th>τ<sub>e</sub></th></tr></thead><tbody><tr><td>None</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>Vol.-dev.</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(6μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>Spectral</td><td>√(2E<sub>0</sub>(1+v<sub>0</sub>)/1+v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td><td>-√(E<sub>0</sub>(1+v<sub>0</sub>)/v<sub>0</sub><sup>2</sup>/s'(0))</td><td>√(4μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>No-tension</td><td>√(2E<sub>0</sub>(1-v<sub>0</sub>)/1-v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td><td>-∞</td><td>√(2E<sub>0</sub>(1-v<sub>0</sub>)/1-v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td></tr><tr><td>DP-like</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))<br>√(18(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/(√6+y)<sup>2</sup>)<br>for 0<y<√(2/3)μ<sub>0</sub><br>for y≥√(2/3)μ<sub>0</sub></td><td>-√(18(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/(√6-y)<sup>2</sup>)<br>-∞<br>for y≥√6</td><td>√(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/s'(0)</td></tr><tr><td>Star-convex</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(2E<sub>0</sub>(1+(1+y*)<sup>(n-1)(n-1)</sup>/R)<sup>2</sup>/w<sub>1</sub>w'(0)/s'(0))<br>-∞<br>for y*≥n(n-1)/2<sup>κ<sub>0</sub></sup>μ<sub>0</sub><br>for y*<sup>κ<sub>0</sub></sup>≥n(n-1)/2<sup>κ<sub>0</sub></sup>μ<sub>0</sub></td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr></tbody></table>

$$ \varphi_D(\boldsymbol{\epsilon}) = \frac{1}{2}\lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_+^2 + \mu_0\boldsymbol{\epsilon}^+ \cdot \boldsymbol{\epsilon}^+, $$

$$ \varphi_R(\boldsymbol{\epsilon}) = \frac{1}{2}\lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_-^2 + \mu_0\boldsymbol{\epsilon}^- \cdot \boldsymbol{\epsilon}^- \quad (34) $$

with $\boldsymbol{\epsilon}^+ = \sum_i \langle \epsilon_i \rangle_+ \boldsymbol{\epsilon}_i \otimes \boldsymbol{\epsilon}_i$ and $\boldsymbol{\epsilon}^- = \sum_i \langle \epsilon_i \rangle_- \boldsymbol{\epsilon}_i \otimes \boldsymbol{\epsilon}_i$, $\epsilon_i$ being the eigenvalues of the strain tensor and $\boldsymbol{\epsilon}_i$ the corresponding eigenvectors. The residual stress tensor

$$ \sigma_R(\boldsymbol{\epsilon}) = \lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_- \mathbf{I} + 2\mu_0\boldsymbol{\epsilon}^- \quad (35) $$

may be non-zero also for tr(ε) > 0. This implies that, even for positive crack opening, non-zero residual stresses are possible. As we will show through a numerical example in Sect. 5, this undesirable coupling between tension and compression behavior is problematic during damage propagation.

The elastic domains $\mathcal{R}(\alpha)$ and $\mathcal{R}^*(\alpha)$ are detailed in Appendix A.2. The stress-based domain $\mathcal{R}^*(\alpha)$ does not evolve as a simple homothety centered in the origin with respect to $\alpha$, see Fig. 3. We plot the elastic domains only in the space of the principal components, because the model cannot be expressed as a function of the volumetric and deviatoric parts of the strain/stress tensors. Since $\mathcal{R}^*(\alpha)$ is not a homotety, we cannot define a unique damage value $\alpha_{\text{peak}}$ which marks the transition between stress-hardening and stress-softening behavior. With increasing damage, $\mathcal{R}^*(\alpha)$ shrinks only along tensile-dominated stress states. On the other hand, the compression-dominated strength does not necessarily have to decrease with increasing damage. Therefore, even though the evolution of $\mathcal{R}^*(\alpha)$ does not fulfill stress-softening for compression-dominated stress states, it still consistently represents damaging behavior. The asymmetry in tension and com-

pression is gained at the expense of the existence of residual stresses. However, the compressive and shear strengths cannot be calibrated independently of the tensile strength even though they are different from those of the standard model (see Table 1).

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Elastic domains for the model with spectral split in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Principal components diagram under the plane-strain assumption</figcaption></figure>
<table><caption>**Table 1** Tensile, compressive and shear strengths for all the analyzed models</caption><thead><tr><th>Decomposition</th><th>σ<sub>e</sub><sup>+</sup></th><th>σ<sub>e</sub><sup>-</sup></th><th>τ<sub>e</sub></th></tr></thead><tbody><tr><td>None</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>Vol.-dev.</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(6μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>Spectral</td><td>√(2E<sub>0</sub>(1+v<sub>0</sub>)/1+v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td><td>-√(E<sub>0</sub>(1+v<sub>0</sub>)/v<sub>0</sub><sup>2</sup>/s'(0))</td><td>√(4μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr><tr><td>No-tension</td><td>√(2E<sub>0</sub>(1-v<sub>0</sub>)/1-v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td><td>-∞</td><td>√(2E<sub>0</sub>(1-v<sub>0</sub>)/1-v<sub>0</sub>-2v<sub>0</sub><sup>2</sup>/s'(0))</td></tr><tr><td>DP-like</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))<br>√(18(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/(√6+y)<sup>2</sup>)<br>for 0<y<√(2/3)μ<sub>0</sub><br>for y≥√(2/3)μ<sub>0</sub></td><td>-√(18(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/(√6-y)<sup>2</sup>)<br>-∞<br>for 0≤y<√6<br>for y≥√6</td><td>√(k<sub>0</sub>y<sup>2</sup>+2μ<sub>0</sub>)w<sub>1</sub>w'(0)/s'(0)</td></tr><tr><td>Star-convex</td><td>√(2E<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td><td>-√(2E<sub>0</sub>(1+(1+y*)<sup>(n-1)</sup>v<sub>0</sub>-1)/<sup>n</sup>(n-1)<sup>2</sup>μ<sub>0</sub>)<br>-∞<br>for y*≥n(n-1)/2μ<sub>0</sub><br>for -1≤y*<n(n-1)/2</td><td>√(2μ<sub>0</sub>w<sub>1</sub>w'(0)/s'(0))</td></tr></tbody></table>

$$ \varphi_D(\boldsymbol{\epsilon}) = \frac{1}{2}\lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_+^2 + \mu_0\boldsymbol{\epsilon}^+ \cdot \boldsymbol{\epsilon}^+, $$

$$ \varphi_R(\boldsymbol{\epsilon}) = \frac{1}{2}\lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_-^2 + \mu_0\boldsymbol{\epsilon}^- \cdot \boldsymbol{\epsilon}^- \quad (34) $$

with $\boldsymbol{\epsilon}^+ = \sum_i \langle \epsilon_i \rangle_+ \boldsymbol{\epsilon}_i \otimes \boldsymbol{\epsilon}_i$ and $\boldsymbol{\epsilon}^- = \sum_i \langle \epsilon_i \rangle_- \boldsymbol{\epsilon}_i \otimes \boldsymbol{\epsilon}_i$, $\epsilon_i$ being the eigenvalues of the strain tensor and $\boldsymbol{\epsilon}_i$ the corresponding eigenvectors. The residual stress tensor

$$ \boldsymbol{\sigma}_R(\boldsymbol{\epsilon}) = \lambda_0(\mathrm{tr}(\boldsymbol{\epsilon}))_- \mathbf{I} + 2\mu_0\boldsymbol{\epsilon}^- \quad (35) $$

may be non-zero also for tr(ε) > 0. This implies that, even for positive crack opening, non-zero residual stresses are possible. As we will show through a numerical example in Sect. 5, this undesirable coupling between tension and compression behavior is problematic during damage propagation.

The elastic domains $R(\alpha)$ and $R^*(\alpha)$ are detailed in Appendix A.2. The stress-based domain $R^*(\alpha)$ does not evolve as a simple homothety centered in the origin with respect to $\alpha$, see Fig. 3. We plot the elastic domains only in the space of the principal components, because the model cannot be expressed as a function of the volumetric and deviatoric parts of the strain/stress tensors. Since $R^*(\alpha)$ is not a homotethy, we cannot define a unique damage value $\alpha_{peak}$ which marks the transition between stress-hardening and stress-softening behavior. With increasing damage, $R^*(\alpha)$ shrinks only along tensile-dominated stress states. On the other hand, the compression-dominated strength does not necessarily have to decrease with increasing damage. Therefore, even though the evolution of $R^*(\alpha)$ does not fulfill stress-softening for compression-dominated stress states, it still consistently represents damaging behavior. The asymmetry in tension and com-

pression is gained at the expense of the existence of residual stresses. However, the compressive and shear strengths cannot be calibrated independently of the tensile strength even though they are different from those of the standard model (see Table 1).

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Elastic domains for the model with spectral split in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Principal components diagram under the plane-strain assumption</figcaption></figure>
### 3.3 The no-tension model

Freddi and Royer-Carfagni in (Freddi and Royer-Carfagni 2010) use the theory of structured deformations of Del Piero and Owen (Del Piero and Owen 1993) to propose a new elastic energy decomposition. They also demonstrate that certain decompositions already available in the literature (Lancioni and Royer-Carfagni 2009; Amor et al. 2009) and the standard model can be derived by leaning on this theory.

The assumption made is that the existence of micro-cracks leads to a reduction in the elastic energy density of the sound material $\varphi_0(\boldsymbol{\varepsilon})$ due to the presence of inelastic deformations known as structured deformations, denoted as $\boldsymbol{\eta}$. These structured deformations are constrained within a convex set $\mathcal{K}_\varepsilon$, which defines the admissible “structure” of the micro-cracks. With $\varphi_0(\boldsymbol{\varepsilon})$ and $\mathcal{K}_\varepsilon$ determined, the residual elastic energy density is computed by solving the following minimization problem

$$ \bar{\eta}(\boldsymbol{\varepsilon}) := \arg \min_{\boldsymbol{\eta} \in \mathcal{K}_\varepsilon} \varphi_0(\boldsymbol{\varepsilon} - \boldsymbol{\eta}), \\ \varphi_R(\boldsymbol{\varepsilon}) := \min_{\boldsymbol{\eta} \in \mathcal{K}_\varepsilon} \varphi_0(\boldsymbol{\varepsilon} - \boldsymbol{\eta}) = \varphi_0(\boldsymbol{\varepsilon} - \bar{\eta}(\boldsymbol{\varepsilon})). \quad (36) $$

The result of the minimization problem (36), reported in Appendix A.3, is applied to the elastic energy density decomposition to derive $\varphi_D(\boldsymbol{\varepsilon}) = \varphi_0(\boldsymbol{\varepsilon}) - \varphi_R(\boldsymbol{\varepsilon})$. By knowing $\varphi_D$ and $\varphi_R$, the residual stresses can be obtained as follows:

$$ - \text{ if } \varepsilon_3 \ge 0, \quad \boldsymbol{\sigma}_R = \mathbf{0}, $$

$$ - \text{ else if } \varepsilon_2 + v_0\varepsilon_3 \ge 0, \quad \boldsymbol{\sigma}_R = \frac{\mu_0(3\lambda_0+2\mu_0)}{\lambda_0+\mu_0}\varepsilon_3\boldsymbol{e}_3 \otimes \boldsymbol{e}_3, $$

$$ - \text{ else if } \varepsilon_1 + \frac{v_0}{1-v_0}(\varepsilon_2 + \varepsilon_3) \ge 0, $$

$$ \boldsymbol{\sigma}_R = \frac{2\mu_0}{\lambda_0+2\mu_0} \left[ - \frac{\lambda_0(3\lambda_0+2\mu_0)(\varepsilon_2+\varepsilon_3)}{\lambda_0+2\mu_0} \boldsymbol{e}_1 \otimes \boldsymbol{e}_1 + (2(\lambda_0+\mu_0)\varepsilon_3 + \lambda_0\varepsilon_2) \boldsymbol{e}_2 \otimes \boldsymbol{e}_2 + (2(\lambda_0+\mu_0)\varepsilon_2 + \lambda_0\varepsilon_3) \boldsymbol{e}_3 \otimes \boldsymbol{e}_3 \right], $$

$$ - \text{ else, } \boldsymbol{\sigma}_R = 2\mu_0\boldsymbol{\varepsilon} + \lambda_0\text{tr}(\boldsymbol{\varepsilon})\mathbf{I}. $$

The elastic domains are detailed in Appendix A.3 and represented in Fig. 4 only in the space of principal components, because the model cannot be expressed as a function of the volumetric and deviatoric parts of the strain/stress tensors. With the no-tension model, the degree of asymmetry between tension and compression behavior is increased compared to the volumetric-deviatoric and the spectral splits. However, the model has some limitations common to the previous ones. The residual stresses lead to the same problematic coupling obtained with the spectral split. Furthermore, Chambolle et al. (2018) study the behavior of this model as

<figure><img src="image_5.png" /><figcaption>Fig. 4 Elastic domains for the no-tension model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}(v_0 = 0.3)$. The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Principal components diagram under the plane-strain assumption</figcaption></figure>

$\ell \to 0$ and demonstrate that a relative displacement of the crack faces results in infinite energy. This constraint does not physically describe a crack. In addition, there is still no flexibility in the choice of the compressive and shear strengths (see Table 1).

### 3.4 The Drucker-Prager-like (DP-like) model

The introduction of structured deformations potentially leads to extra flexibility, namely in the selection of the convex set $\mathcal{K}_\varepsilon$. De Lorenzis and Maurini (De Lorenzis and Maurini 2022) leverage this flexibility by introducing a new decomposition that depends on an additional material parameter $\gamma$. This parameter allows for elastic domains which partially recover the Drucker-Prager model, commonly used for modeling compressive failure in cohesive-frictional materials like rocks or concrete, and enables the independent calibration of the tensile and the compressive or shear strengths.

The residual energy is determined through the structured deformation problem (36), utilizing a convex cone of admissible structured deformations

$$ \mathcal{K}_\varepsilon := \{\boldsymbol{\eta} \in \text{Sym} : \text{tr}(\boldsymbol{\eta}) \ge \gamma | \boldsymbol{\eta}_{\text{dev}}|\} \quad \text{with } \gamma \ge 0. \quad (37) $$
### 3.3 The no-tension model

Freddi and Royer-Carfagni in (Freddi and Royer-Carfagni 2010) use the theory of structured deformations of Del Piero and Owen (Del Piero and Owen 1993) to propose a new elastic energy decomposition. They also demonstrate that certain decompositions already available in the literature (Lancioni and Royer-Carfagni 2009; Amor et al. 2009) and the standard model can be derived by leaning on this theory.

The assumption made is that the existence of micro-cracks leads to a reduction in the elastic energy density of the sound material $\varphi_0(\boldsymbol{\varepsilon})$ due to the presence of inelastic deformations known as structured deformations, denoted as $\boldsymbol{\eta}$. These structured deformations are constrained within a convex set $\mathcal{K}_\varepsilon$, which defines the admissible “structure” of the micro-cracks. With $\varphi_0(\boldsymbol{\varepsilon})$ and $\mathcal{K}_\varepsilon$ determined, the residual elastic energy density is computed by solving the following minimization problem

$$ \bar{\eta}(\boldsymbol{\varepsilon}) := \arg \min_{\boldsymbol{\eta} \in \mathcal{K}_\varepsilon} \varphi_0(\boldsymbol{\varepsilon} - \boldsymbol{\eta}), \\ \varphi_R(\boldsymbol{\varepsilon}) := \min_{\boldsymbol{\eta} \in \mathcal{K}_\varepsilon} \varphi_0(\boldsymbol{\varepsilon} - \boldsymbol{\eta}) = \varphi_0(\boldsymbol{\varepsilon} - \bar{\eta}(\boldsymbol{\varepsilon})). \quad (36) $$

The result of the minimization problem (36), reported in Appendix A.3, is applied to the elastic energy density decomposition to derive $\varphi_D(\boldsymbol{\varepsilon}) = \varphi_0(\boldsymbol{\varepsilon}) - \varphi_R(\boldsymbol{\varepsilon})$. By knowing $\varphi_D$ and $\varphi_R$, the residual stresses can be obtained as follows:

$$ - \text{ if } \varepsilon_3 \ge 0, \quad \boldsymbol{\sigma}_R = \mathbf{0}, $$

$$ - \text{ else if } \varepsilon_2 + v_0\varepsilon_3 \ge 0, \quad \boldsymbol{\sigma}_R = \frac{\mu_0(3\lambda_0 + 2\mu_0)}{\lambda_0 + \mu_0}\varepsilon_3\boldsymbol{e}_3 \otimes \boldsymbol{e}_3, $$

$$ - \text{ else if } \varepsilon_1 + \frac{v_0}{1-v_0}(\varepsilon_2 + \varepsilon_3) \ge 0, $$

$$ \boldsymbol{\sigma}_R = \frac{2\mu_0}{\lambda_0+2\mu_0} \left[ - \frac{\lambda_0(3\lambda_0+2\mu_0)(\varepsilon_2+\varepsilon_3)}{\lambda_0+2\mu_0} \boldsymbol{e}_1 \otimes \boldsymbol{e}_1 + (2(\lambda_0+\mu_0)\varepsilon_3 + \lambda_0\varepsilon_2) \boldsymbol{e}_2 \otimes \boldsymbol{e}_2 + (2(\lambda_0+\mu_0)\varepsilon_2 + \lambda_0\varepsilon_3) \boldsymbol{e}_3 \otimes \boldsymbol{e}_3 \right], $$

$$ - \text{ else, } \boldsymbol{\sigma}_R = 2\mu_0\boldsymbol{\varepsilon} + \lambda_0\text{tr}(\boldsymbol{\varepsilon})\mathbf{I}. $$

The elastic domains are detailed in Appendix A.3 and represented in Fig. 4 only in the space of principal components, because the model cannot be expressed as a function of the volumetric and deviatoric parts of the strain/stress tensors. With the no-tension model, the degree of asymmetry between tension and compression behavior is increased compared to the volumetric-deviatoric and the spectral splits. However, the model has some limitations common to the previous ones. The residual stresses lead to the same problematic coupling obtained with the spectral split. Furthermore, Chambolle et al. (2018) study the behavior of this model as

<figure><img src="image_5.png" /><figcaption>Fig. 4 Elastic domains for the no-tension model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}(v_0 = 0.3)$. The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Principal components diagram under the plane-strain assumption</figcaption></figure>

$\ell \to 0$ and demonstrate that a relative displacement of the crack faces results in infinite energy. This constraint does not physically describe a crack. In addition, there is still no flexibility in the choice of the compressive and shear strengths (see Table 1).

### 3.4 The Drucker-Prager-like (DP-like) model

The introduction of structured deformations potentially leads to extra flexibility, namely in the selection of the convex set $\mathcal{K}_\varepsilon$. De Lorenzis and Maurini (De Lorenzis and Maurini 2022) leverage this flexibility by introducing a new decomposition that depends on an additional material parameter $\gamma$. This parameter allows for elastic domains which partially recover the Drucker-Prager model, commonly used for modeling compressive failure in cohesive-frictional materials like rocks or concrete, and enables the independent calibration of the tensile and the compressive or shear strengths.

The residual energy is determined through the structured deformation problem (36), utilizing a convex cone of admissible structured deformations

$$ \mathcal{K}_\varepsilon := \{\boldsymbol{\eta} \in \text{Sym} : \text{tr}(\boldsymbol{\eta}) \ge \gamma | \boldsymbol{\eta}_{\text{dev}}|\} \quad \text{with } \gamma \ge 0. \quad (37) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5 Elastic domains for the DP-like model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma = \sqrt{2\mu_0/\kappa_0}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength sur-</figcaption></figure>

faces). The grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram a. Principal components diagram under the plane-strain assumption b.

The elastic strain energy density decomposition, solution of the minimization problem (36), is detailed in Appendix A.4. The residual stresses read as follows:

- if $|\boldsymbol{\epsilon}_{\text{dev}}| < \text{tr}(\boldsymbol{\epsilon})/\gamma$, $\boldsymbol{\sigma}_R = \mathbf{0}$,

- else if $|\boldsymbol{\epsilon}_{\text{dev}}| \ge -\frac{\gamma\kappa_0}{2\mu_0}\text{tr}(\boldsymbol{\epsilon})$,

$$ \boldsymbol{\sigma}_R = \frac{2\kappa_0\mu_0}{\kappa_0\gamma^2+2\mu_0} \left( (\text{tr}\boldsymbol{\epsilon} - \gamma|\boldsymbol{\epsilon}_{\text{dev}}|)\mathbf{I} + \gamma \left( \gamma - \frac{\text{tr}\boldsymbol{\epsilon}}{|\boldsymbol{\epsilon}_{\text{dev}}|} \right) \boldsymbol{\epsilon}_{\text{dev}} \right), $$

- else, $\boldsymbol{\sigma}_R = 2\mu_0\boldsymbol{\epsilon}_{\text{dev}} + \kappa_0\text{tr}(\boldsymbol{\epsilon})\mathbf{I}$.

With this model, the shear or the compressive strength can be calibrated independently of the tensile strength through the new parameter $\gamma$ (see Figs. 5 and 6). However, during damage evolution the residual stresses lead to the same undesired coupling between tension and compression as with the spectral split and the no-tension model.

## 3.5 Summary

In light of the requirements defined in Sect. 2.4, we summarize the merits and disadvantages of the ana-

lyzed models in Table 2. None of the available splits meets all requirements, as they either lack in flexibility (a requisite that is most important for nucleation, since this defines the strength surface of the material) or give rise to residual stresses which contain spurious components (a requisite that is most important for damage evolution and for obtaining a crack-like behavior of the fully developed damage localization bands). With the goal of satisfying the above requirements under both aspects of flexibility and crack-like residual stresses, in the next section we introduce a novel model, which we denote as the star-convex model.

## 4 The star-convex model

In this section, we propose a new model that aims at satisfying the requirements defined in Sect. 2.4. As mentioned earlier, in (Chambolle et al. 2018) it is shown that the volumetric-deviatoric split can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence. This success hinges on a residual elas-
<figure><img src="image_6.png" /><figcaption>Fig. 5 Elastic domains for the DP-like model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma = \sqrt{2\mu_0/\kappa_0}$ ($\nu_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength sur-</figcaption></figure>

faces). The grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram a. Principal components diagram under the plane-strain assumption b

The elastic strain energy density decomposition, solution of the minimization problem (36), is detailed in Appendix A.4. The residual stresses read as follows:

- if $|\boldsymbol{\epsilon}_{\text{dev}}| < \text{tr}(\boldsymbol{\epsilon})/\gamma$, $\boldsymbol{\sigma}_R = \mathbf{0}$,

- else if $|\boldsymbol{\epsilon}_{\text{dev}}| \ge -\frac{\gamma\kappa_0}{2\mu_0}\text{tr}(\boldsymbol{\epsilon})$,

$$ \boldsymbol{\sigma}_R = \frac{2\kappa_0\mu_0}{\kappa_0\gamma^2+2\mu_0} \left( (\text{tr}\boldsymbol{\epsilon} - \gamma|\boldsymbol{\epsilon}_{\text{dev}}|)\mathbf{I} + \gamma \left( \gamma - \frac{\text{tr}\boldsymbol{\epsilon}}{|\boldsymbol{\epsilon}_{\text{dev}}|} \right) \boldsymbol{\epsilon}_{\text{dev}} \right), $$

- else, $\boldsymbol{\sigma}_R = 2\mu_0\boldsymbol{\epsilon}_{\text{dev}} + \kappa_0\text{tr}(\boldsymbol{\epsilon})\mathbf{I}$.

With this model, the shear or the compressive strength can be calibrated independently of the tensile strength through the new parameter $\gamma$ (see Figs. 5 and 6). However, during damage evolution the residual stresses lead to the same undesired coupling between tension and compression as with the spectral split and the no-tension model.

## 3.5 Summary

In light of the requirements defined in Sect. 2.4, we summarize the merits and disadvantages of the ana-

lyzed models in Table 2. None of the available splits meets all requirements, as they either lack in flexibility (a requisite that is most important for nucleation, since this defines the strength surface of the material) or give rise to residual stresses which contain spurious components (a requisite that is most important for damage evolution and for obtaining a crack-like behavior of the fully developed damage localization bands). With the goal of satisfying the above requirements under both aspects of flexibility and crack-like residual stresses, in the next section we introduce a novel model, which we denote as the star-convex model.

## 4 The star-convex model

In this section, we propose a new model that aims at satisfying the requirements defined in Sect. 2.4. As mentioned earlier, in (Chambolle et al. 2018) it is shown that the volumetric-deviatoric split can be used to approximate the non-interpenetration constraint in the sense of $\Gamma$-convergence. This success hinges on a residual elas-
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Elastic domains for the DP-like model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma = \sqrt{10\mu_0/\kappa_0}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength</figcaption></figure>

surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram **a**. Principal components diagram under the plane-strain assumption **b**

<caption>**Table 2** Analyzed models and requirements for an ideal model under multi-axial stress states: blue checks and red crosses denote respectively success and failure in satisfying a requirement. For the sake of compactness the last line of the table includes the</caption>

star-convex model even though it is first introduced in the following section. The star-convex model meets all requirements with partial flexibility

<table><thead><tr><th></th><th>Strain-hardening</th><th>Stress-softening</th><th>Tens./compr. asymmetry</th><th>Flexibility</th><th>Crack-like residual stress</th></tr></thead><tbody><tr><td>Standard</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>Vol.-dev.</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Spectral</td><td>✓</td><td>✗</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>No-tension</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>DP-like</td><td>✓</td><td>✓</td><td>✓</td><td>✓<sup>a</sup></td><td>✗</td></tr><tr><td>Star-convex</td><td>✓</td><td>✓</td><td>✓</td><td>✓<sup>b</sup></td><td>✓</td></tr></tbody></table>

tic strain energy density depending solely on the volumetric energy associated to a negative strain trace, i.e., $\phi_R(\epsilon) \propto \langle\text{tr}(\epsilon)\rangle^2$. All the other decompositions suffer from spurious non-zero tangential stiffness for fully developed cracks. Therefore, the volumetric-deviatoric split represents a sound option to model crack-like

stresses in the phase-field framework, but lacks flexibility in calibrating the compressive strength $\sigma_e^-$ independently of the tensile strength $\sigma_e^+$. To compensate for this disadvantage but keep a residual energy density (and stress) exclusively related to a local volume

<footnote><sup>a</sup> Partial flexibility: out of tensile, compressive and shear strengths, the model allows for the independent calibration of any pair, while the third property is automatically fixed</footnote>

<footnote><sup>b</sup> Partial flexibility: the model allows for the independent calibration of tensile and compressive or shear and compressive strengths, while the third property is automatically fixed</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Elastic domains for the DP-like model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma = \sqrt{10\mu_0/\kappa_0}$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength</figcaption></figure>

surfaces). The grey background indicates the parts of the elastic domains shared with the standard model. Volumetric-deviatoric diagram **a**. Principal components diagram under the plane-strain assumption **b**

<caption>**Table 2** Analyzed models and requirements for an ideal model under multi-axial stress states: blue checks and red crosses denote respectively success and failure in satisfying a requirement. For the sake of compactness the last line of the table includes the</caption>

star-convex model even though it is first introduced in the following section. The star-convex model meets all requirements with partial flexibility

<table><thead><tr><th></th><th>Strain-hardening</th><th>Stress-softening</th><th>Tens./compr. asymmetry</th><th>Flexibility</th><th>Crack-like residual stress</th></tr></thead><tbody><tr><td>Standard</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td><td>✗</td></tr><tr><td>Vol.-dev.</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✓</td></tr><tr><td>Spectral</td><td>✓</td><td>✗</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>No-tension</td><td>✓</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>DP-like</td><td>✓</td><td>✓</td><td>✓</td><td>✓<sup>a</sup></td><td>✗</td></tr><tr><td>Star-convex</td><td>✓</td><td>✓</td><td>✓</td><td>✓<sup>b</sup></td><td>✓</td></tr></tbody></table>

tic strain energy density depending solely on the volumetric energy associated to a negative strain trace, i.e., $\phi_R(\varepsilon) \propto \langle\text{tr}(\varepsilon)\rangle^2$. All the other decompositions suffer from spurious non-zero tangential stiffness for fully developed cracks. Therefore, the volumetric-deviatoric split represents a sound option to model crack-like

stresses in the phase-field framework, but lacks flexibility in calibrating the compressive strength $\sigma_e^-$ independently of the tensile strength $\sigma_e^+$. To compensate for this disadvantage but keep a residual energy density (and stress) exclusively related to a local volume

<footnote><sup>a</sup> Partial flexibility: out of tensile, compressive and shear strengths, the model allows for the independent calibration of any pair, while the third property is automatically fixed</footnote>

<footnote><sup>b</sup> Partial flexibility: the model allows for the independent calibration of tensile and compressive or shear and compressive strengths, while the third property is automatically fixed</footnote>
contraction $\langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_-$, we propose the following energy decomposition

$$ \varphi_D(\boldsymbol{\varepsilon}) = \mu_0 |\boldsymbol{\varepsilon}_{\text{dev}}|^2 + \frac{1}{2} \kappa_0 \left( \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_+^2 - \gamma^* \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_-^2 \right), \quad (38) $$

$$ \varphi_R(\boldsymbol{\varepsilon}) = (1 + \gamma^*) \frac{1}{2} \kappa_0 \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_-^2, $$

where $\gamma^* \ge -1$ is the additional parameter controlling the $\sigma_e^- / \sigma_e^+$ ratio. Hence, the volumetric and deviatoric components of the stress are derived as

$$ \begin{aligned} \sigma_{\text{vol}} &= \kappa_0 \left\{ a(\alpha) \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_+ + \left[ 1 + \gamma^* (1 - a(\alpha)) \right] \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_- \right\} \mathbf{I}, \\ \sigma_{\text{dev}} &= 2a(\alpha) \mu_0 \boldsymbol{\varepsilon}_{\text{dev}}. \end{aligned} \quad (39) $$

For a reason that will become clear later, we denote the corresponding model as star-convex model (or star-convex energy decomposition). In the following subsections, we analyze the elastic domains for the proposed model in the strain and stress spaces. In Appendix B, we provide further insights on the development of the proposed model in comparison to alternative options.

## 4.1 Strain space

Figures 7, 8 and 9 illustrate the elastic domain of the proposed model in the strain space, which is given by the conditions

$$ \left\{ \begin{array}{ll} \mu_0 |\boldsymbol{\varepsilon}_{\text{dev}}|^2 + \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\varepsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)}, & \text{for } \text{tr}(\boldsymbol{\varepsilon}) \ge 0 \\ \mu_0 |\boldsymbol{\varepsilon}_{\text{dev}}|^2 - \gamma^* \frac{\kappa_0}{2} \text{tr}^2(\boldsymbol{\varepsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)}, & \text{for } \text{tr}(\boldsymbol{\varepsilon}) < 0. \end{array} \right. $$

For $\text{tr}(\boldsymbol{\varepsilon}) \ge 0$, the new model shares the same boundary of the elastic domain $\partial \mathcal{R}(\alpha)$ with the standard model and $\partial \mathcal{R}(\alpha)$ is an ellipse in the $\text{tr}(\boldsymbol{\varepsilon}) - |\boldsymbol{\varepsilon}_{\text{dev}}|$ diagram. On the same diagram but for $\text{tr}(\boldsymbol{\varepsilon}) < 0$, the parameter $\gamma^*$ determines the type of conic section (Fig. 7) such that

* for $-1 \le \gamma^* < 0$, $\partial \mathcal{R}(\alpha)$ is an ellipse,

* for $\gamma^* = 0$, $\partial \mathcal{R}(\alpha)$ is a degenerate parabola (horizontal straight line),

* for $\gamma^* > 0$, $\partial \mathcal{R}(\alpha)$ is a hyperbola lying above the asymptote $|\boldsymbol{\varepsilon}_{\text{dev}}| = -\sqrt{\frac{\kappa_0}{2\mu_0}} \gamma^* \text{tr}(\boldsymbol{\varepsilon})$

In particular, we retrieve the standard model when $\gamma^* = -1$ and the volumetric-deviatoric split for $\gamma^* = 0$.

For $-1 \le \gamma^* \le 0$, the elastic domain $\mathcal{R}(\alpha)$ is convex, whereas for $\gamma^* > 0$ $\partial \mathcal{R}(\alpha)$ consists of the ellipse defined over $\text{tr}(\boldsymbol{\varepsilon}) > 0$ smoothly joined at $\text{tr}(\boldsymbol{\varepsilon}) = 0$ with the hyperbola defined over $\text{tr}(\boldsymbol{\varepsilon}) < 0$. Hence, in

this case $\mathcal{R}(\alpha)$ is not convex but rather 0-star-convex, i.e. $\forall \boldsymbol{\varepsilon} \in \mathcal{R}(\alpha)$ and $\forall s \in [0, 1]$, $s\boldsymbol{\varepsilon} \in \mathcal{R}(\alpha)$ Hansen et al. (2020). This is the reason for denoting the new model as star-convex model (or star-convex decomposition).

Note that starting from the Drucker-Ilyushin postulate, which requires the non-negativity of the interior work in a strain cycle, it can be proved that an elasto-plastic material must have a convex elastic domain. However, the same thermodynamic postulate for an elasto-damaging material does not necessarily imply convexity, but at most star-convexity with respect to 0 (Marigo 1989).

The star-convex model enjoys the strain-hardening property $\forall \alpha \in [0, 1]$. Additionally, for $\gamma^* > 0$, the hyperbola for $\text{tr}(\boldsymbol{\varepsilon}) < 0$ evolves but keeps a constant asymptote as this depends only on the elastic constants and on $\gamma^*$. Since $\partial \mathcal{R}(\alpha)$ lies above the asymptote (grey dashed line in Fig. 7), the strain states below the asymptote are unable to produce additional damage and the material behaves as linearly elastic. Accordingly, the compressive strength can be calibrated by modifying the slope of the asymptote by tuning $\gamma^*$, thus providing the desired flexibility in nucleation. At the limit $\alpha = 1$, the star-convex elastic domain becomes $\mathcal{R}(1) = \text{Sym}$, meaning that for a fully broken material the whole strain space is admissible.

## 4.2 Stress space

The elastic domain in the stress space $\mathcal{R}^*(\alpha)$ (Figs. 8, 9) is the set of $\boldsymbol{\sigma} \in \text{Sym}$ such that

$$ \begin{aligned} - & \text{ if } \text{tr}(\boldsymbol{\sigma}) \ge 0, \quad \frac{1}{4\mu_0} |\boldsymbol{\sigma}_{\text{dev}}|^2 + \frac{1}{2\kappa_0 n^2} \text{tr}^2(\boldsymbol{\sigma}) \le \frac{w_1 w'(\alpha)}{s'(\alpha)}, \\ - & \text{ if } \text{tr}(\boldsymbol{\sigma}) < 0, \quad \frac{1}{4\mu_0} |\boldsymbol{\sigma}_{\text{dev}}|^2 + \\ & \qquad - \gamma^* \frac{a(\alpha)^2}{2\kappa_0 n^2 (1+\gamma^*(1-a(\alpha)))^2} \text{tr}^2(\boldsymbol{\sigma}) \le \frac{w_1 w'(\alpha)}{s'(\alpha)}. \end{aligned} $$

Similarly to the strain space, for $\text{tr}(\boldsymbol{\sigma}) \ge 0$ the star-convex split shares with the standard model the same elliptic boundary $\partial \mathcal{R}^*(\alpha)$ in the $\text{tr}(\boldsymbol{\sigma}) - |\boldsymbol{\sigma}_{\text{dev}}|$ diagram, whereas for $\text{tr}(\boldsymbol{\sigma}) < 0$

* for $-1 \le \gamma^* < 0$, $\partial \mathcal{R}^*(\alpha)$ is an ellipse,

* for $\gamma^* = 0$, $\partial \mathcal{R}^*(\alpha)$ is a degenerate parabola (horizontal straight line),

* for $\gamma^* > 0$, $\partial \mathcal{R}^*(\alpha)$ is a hyperbola lying above the asymptote $|\boldsymbol{\sigma}_{\text{dev}}| = -\frac{2a(\alpha)}{n(1+\gamma^*(1-a(\alpha)))} \sqrt{\frac{2\mu_0}{\kappa_0}} \gamma^* \text{tr}(\boldsymbol{\sigma})$.

A three-dimensional representation of $\mathcal{R}^*(0)$ is given for $\gamma^* = 1, 5$ in Fig. 10. For $-1 \le \gamma^* \le 0$, the
contraction $\langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_-$, we propose the following energy decomposition

$$ \varphi_D(\boldsymbol{\varepsilon}) = \mu_0 |\boldsymbol{\varepsilon}_{\mathrm{dev}}|^2 + \frac{1}{2} \kappa_0 \left( \langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_+^2 - \gamma^* \langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_-^2 \right), \quad (38) $$

$$ \varphi_R(\boldsymbol{\varepsilon}) = (1 + \gamma^*) \frac{1}{2} \kappa_0 \langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_-^2, $$

where $\gamma^* \ge -1$ is the additional parameter controlling the $\sigma_e^- / \sigma_e^+$ ratio. Hence, the volumetric and deviatoric components of the stress are derived as

$$ \begin{aligned} \sigma_{\mathrm{vol}} &= \kappa_0 \{ a(\alpha) \langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_+ + [1 + \gamma^* (1 - a(\alpha))] \langle \mathrm{tr}(\boldsymbol{\varepsilon}) \rangle_- \} \mathbf{I}, \\ \sigma_{\mathrm{dev}} &= 2a(\alpha) \mu_0 \boldsymbol{\varepsilon}_{\mathrm{dev}}. \end{aligned} \quad (39) $$

For a reason that will become clear later, we denote the corresponding model as star-convex model (or star-convex energy decomposition). In the following subsections, we analyze the elastic domains for the proposed model in the strain and stress spaces. In Appendix B, we provide further insights on the development of the proposed model in comparison to alternative options.

## 4.1 Strain space

Figures 7, 8 and 9 illustrate the elastic domain of the proposed model in the strain space, which is given by the conditions

$$ \begin{cases} \mu_0 |\boldsymbol{\varepsilon}_{\mathrm{dev}}|^2 + \frac{\kappa_0}{2} \mathrm{tr}^2(\boldsymbol{\varepsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)}, & \text{for } \mathrm{tr}(\boldsymbol{\varepsilon}) \ge 0 \\ \mu_0 |\boldsymbol{\varepsilon}_{\mathrm{dev}}|^2 - \gamma^* \frac{\kappa_0}{2} \mathrm{tr}^2(\boldsymbol{\varepsilon}) \le -\frac{w_1 w'(\alpha)}{a'(\alpha)}, & \text{for } \mathrm{tr}(\boldsymbol{\varepsilon}) < 0. \end{cases} $$

For $\mathrm{tr}(\boldsymbol{\varepsilon}) \ge 0$, the new model shares the same boundary of the elastic domain $\partial \mathcal{R}(\alpha)$ with the standard model and $\partial \mathcal{R}(\alpha)$ is an ellipse in the $\mathrm{tr}(\boldsymbol{\varepsilon}) - |\boldsymbol{\varepsilon}_{\mathrm{dev}}|$ diagram. On the same diagram but for $\mathrm{tr}(\boldsymbol{\varepsilon}) < 0$, the parameter $\gamma^*$ determines the type of conic section (Fig. 7) such that

* for $-1 \le \gamma^* < 0$, $\partial \mathcal{R}(\alpha)$ is an ellipse,

* for $\gamma^* = 0$, $\partial \mathcal{R}(\alpha)$ is a degenerate parabola (horizontal straight line),

* for $\gamma^* > 0$, $\partial \mathcal{R}(\alpha)$ is a hyperbola lying above the asymptote $|\boldsymbol{\varepsilon}_{\mathrm{dev}}| = -\sqrt{\frac{\kappa_0}{2\mu_0}} \gamma^* \mathrm{tr}(\boldsymbol{\varepsilon})$

In particular, we retrieve the standard model when $\gamma^* = -1$ and the volumetric-deviatoric split for $\gamma^* = 0$.

For $-1 \le \gamma^* \le 0$, the elastic domain $\mathcal{R}(\alpha)$ is convex, whereas for $\gamma^* > 0$ $\partial \mathcal{R}(\alpha)$ consists of the ellipse defined over $\mathrm{tr}(\boldsymbol{\varepsilon}) > 0$ smoothly joined at $\mathrm{tr}(\boldsymbol{\varepsilon}) = 0$ with the hyperbola defined over $\mathrm{tr}(\boldsymbol{\varepsilon}) < 0$. Hence, in

this case $\mathcal{R}(\alpha)$ is not convex but rather 0-star-convex, i.e. $\forall \boldsymbol{\varepsilon} \in \mathcal{R}(\alpha)$ and $\forall s \in [0, 1]$, $s\boldsymbol{\varepsilon} \in \mathcal{R}(\alpha)$ Hansen et al. (2020). This is the reason for denoting the new model as star-convex model (or star-convex decomposition).

Note that starting from the Drucker-Ilyushin postulate, which requires the non-negativity of the interior work in a strain cycle, it can be proved that an elasto-plastic material must have a convex elastic domain. However, the same thermodynamic postulate for an elasto-damaging material does not necessarily imply convexity, but at most star-convexity with respect to 0 (Marigo 1989).

The star-convex model enjoys the strain-hardening property $\forall \alpha \in [0, 1]$. Additionally, for $\gamma^* > 0$, the hyperbola for $\mathrm{tr}(\boldsymbol{\varepsilon}) < 0$ evolves but keeps a constant asymptote as this depends only on the elastic constants and on $\gamma^*$. Since $\partial \mathcal{R}(\alpha)$ lies above the asymptote (grey dashed line in Fig. 7), the strain states below the asymptote are unable to produce additional damage and the material behaves as linearly elastic. Accordingly, the compressive strength can be calibrated by modifying the slope of the asymptote by tuning $\gamma^*$, thus providing the desired flexibility in nucleation. At the limit $\alpha = 1$, the star-convex elastic domain becomes $\mathcal{R}(1) = \mathrm{Sym}$, meaning that for a fully broken material the whole strain space is admissible.

## 4.2 Stress space

The elastic domain in the stress space $\mathcal{R}^*(\alpha)$ (Figs. 8, 9) is the set of $\boldsymbol{\sigma} \in \mathrm{Sym}$ such that

$$ \begin{aligned} & - \text{if } \mathrm{tr}(\boldsymbol{\sigma}) \ge 0, \frac{1}{4\mu_0} |\boldsymbol{\sigma}_{\mathrm{dev}}|^2 + \frac{1}{2\kappa_0 n^2} \mathrm{tr}^2(\boldsymbol{\sigma}) \le \frac{w_1 w'(\alpha)}{s'(\alpha)}, \\ & - \text{if } \mathrm{tr}(\boldsymbol{\sigma}) < 0, \frac{1}{4\mu_0} |\boldsymbol{\sigma}_{\mathrm{dev}}|^2 + \\ & \qquad - \gamma^* \frac{a(\alpha)^2}{2\kappa_0 n^2 (1+\gamma^*(1-a(\alpha)))^2} \mathrm{tr}^2(\boldsymbol{\sigma}) \le \frac{w_1 w'(\alpha)}{s'(\alpha)}. \end{aligned} $$

Similarly to the strain space, for $\mathrm{tr}(\boldsymbol{\sigma}) \ge 0$ the star-convex split shares with the standard model the same elliptic boundary $\partial \mathcal{R}^*(\alpha)$ in the $\mathrm{tr}(\boldsymbol{\sigma}) - |\boldsymbol{\sigma}_{\mathrm{dev}}|$ diagram, whereas for $\mathrm{tr}(\boldsymbol{\sigma}) < 0$

* for $-1 \le \gamma^* < 0$, $\partial \mathcal{R}^*(\alpha)$ is an ellipse,

* for $\gamma^* = 0$, $\partial \mathcal{R}^*(\alpha)$ is a degenerate parabola (horizontal straight line),

* for $\gamma^* > 0$, $\partial \mathcal{R}^*(\alpha)$ is a hyperbola lying above the asymptote $|\boldsymbol{\sigma}_{\mathrm{dev}}| = -\frac{2a(\alpha)}{n(1+\gamma^*(1-a(\alpha)))} \sqrt{\frac{2\mu_0}{\kappa_0}} \gamma^* \mathrm{tr}(\boldsymbol{\sigma})$.

A three-dimensional representation of $\mathcal{R}^*(0)$ is given for $\gamma^* = 1, 5$ in Fig. 10. For $-1 \le \gamma^* \le 0$, the
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Nucleation domain in the strain space for the star-convex model with varying $\gamma^*$ ($v_0 = 0.3$). For $\text{tr}(\boldsymbol{\varepsilon}) < 0$, the boundary is an ellipse if $-1 \le \gamma^* < 0$, a degenerate parabola if $\gamma^* = 0$ and a hyperbola (with the dashed grey line as asymptote) if $\gamma^* > 0$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Elastic domains for the star-convex model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ for $\gamma^* = 1$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared</figcaption></figure>

domain $\mathcal{R}^*(\alpha)$ is convex, whereas for $\gamma^* > 0$ $\mathcal{R}^*(\alpha)$ is not convex but **0**-star-convex, i.e. $\forall \sigma \in \mathcal{R}^*(\alpha)$ and $\forall s \in [0, 1]$, $s\sigma \in \mathcal{R}^*(\alpha)$ Hansen et al. (2020).

The star-convex model enjoys the stress-softening property $\forall \alpha \in [0, 1)$. At the limit $\alpha = 1$, the star-convex elastic domain becomes strictly convex and collapses to the negative hydrostatic pressure half-line, meaning that non-zero deviatoric and positive volumetric stresses are not admissible in a fully broken material. In particular, according to (27), the stress tensor when $\alpha = 1$ is

with the standard model. Volumetric-deviatoric diagram **a**: damage cannot evolve for loading directions below the black dashed line in the strain space which is also the constant asymptote. Principal components diagram under the plane-strain assumption **b**

$$ \sigma_R = (1 + \gamma^*) \kappa_0 \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_{-} \mathbf{I}, \quad (40) $$

thus, the mapping between admissible strains and admissible stresses for $\alpha = 1$ is such that only negative volumetric strains can produce work through negative volumetric stresses. In this manner, the model fulfils the crack-like stress requirement.

Equation (40) shows that, for negative volumetric strains, the volumetric stiffness for $\alpha = 1$ is $\kappa_0(1 + \gamma^*)$. On the other hand, according to (39), it is equal to $\kappa_0$ when $\alpha = 0$. Therefore, for negative volumetric strains, the volumetric stiffness of a fully broken material is larger than the one of the pristine material when $\gamma^* > 0$. This increase can be explained by consider-
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Nucleation domain in the strain space for the star-convex model with varying $\gamma^*$ ($v_0 = 0.3$). For $\text{tr}(\boldsymbol{\varepsilon}) < 0$, the boundary is an ellipse if $-1 \le \gamma^* < 0$, a degenerate parabola if $\gamma^* = 0$ and a hyperbola (with the dashed grey line as asymptote) if $\gamma^* > 0$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8** Elastic domains for the star-convex model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ for $\gamma^* = 1$ ($v_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared</figcaption></figure>

domain $\mathcal{R}^*(\alpha)$ is convex, whereas for $\gamma^* > 0$ $\mathcal{R}^*(\alpha)$ is not convex but **0**-star-convex, i.e. $\forall \sigma \in \mathcal{R}^*(\alpha)$ and $\forall s \in [0, 1]$, $s\sigma \in \mathcal{R}^*(\alpha)$ Hansen et al. (2020).

The star-convex model enjoys the stress-softening property $\forall \alpha \in [0, 1)$. At the limit $\alpha = 1$, the star-convex elastic domain becomes strictly convex and collapses to the negative hydrostatic pressure half-line, meaning that non-zero deviatoric and positive volumetric stresses are not admissible in a fully broken material. In particular, according to (27), the stress tensor when $\alpha = 1$ is

with the standard model. Volumetric-deviatoric diagram **a**: damage cannot evolve for loading directions below the black dashed line in the strain space which is also the constant asymptote. Principal components diagram under the plane-strain assumption **b**

$$ \sigma_R = (1 + \gamma^*) \kappa_0 \langle \text{tr}(\boldsymbol{\varepsilon}) \rangle_{-} \mathbf{I}, \quad (40) $$

thus, the mapping between admissible strains and admissible stresses for $\alpha = 1$ is such that only negative volumetric strains can produce work through negative volumetric stresses. In this manner, the model fulfils the crack-like stress requirement.

Equation (40) shows that, for negative volumetric strains, the volumetric stiffness for $\alpha = 1$ is $\kappa_0(1 + \gamma^*)$. On the other hand, according to (39), it is equal to $\kappa_0$ when $\alpha = 0$. Therefore, for negative volumetric strains, the volumetric stiffness of a fully broken material is larger than the one of the pristine material when $\gamma^* > 0$. This increase can be explained by consider-
<figure><img src="image_10.png" /><figcaption>Fig. 9 Elastic domains for the star-convex model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma^* = 5$ ($\nu_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared</figcaption></figure>

with the standard model. Volumetric-deviatoric diagram a: damage cannot evolve for loading directions below the black dashed line in the strain space which is also the constant asymptote. Principal components diagram under the plane-strain assumption b

ing that the damage evolution occurs only above the asymptote, where the deviatoric contribution exceeds the volumetric one. Damage evolution redistributes the stiffness such that the deviatoric stiffness decreases to zero and the volumetric stiffness increases, while the total elastic energy density (which is the sum of the volumetric and deviatoric contributions) always decreases.

Defining as usual the strength surface as the boundary of the elastic domain in the stress space for $\alpha = \alpha_{peak} = 0$, i.e. $\partial \mathcal{R}^*(0)$, we can derive the tensile strength, the shear strength and the compressive strength as

$$ 
\begin{aligned}
\sigma_e^+ &= \sqrt{\frac{2 E_0 w_1 w'(0)}{s'(0)}}, & \tau_e &= \sqrt{\frac{2 \mu_0 w_1 w'(0)}{s'(0)}}, \\
\sigma_e^- &= \begin{cases} -\sqrt{\frac{1}{1+(1+\gamma^*)^{\frac{(n-1)(n-1)}{n}} \frac{2 E_0 w_1 w'(0)}{s'(0)}}} & \text{for } \gamma^* < \frac{n(n-1)}{2} \frac{\kappa_0}{\mu_0}, \\ -\infty & \text{otherwise.} \end{cases}
\end{aligned} 
\quad (41) 
$$

From (41) we see that the shear strength $\tau_e$ is independent of $\gamma^*$. Thus, the flexibility offered by the star-convex model is exclusively related to the compressive strength $\sigma_e^-$.

## 5 Numerical tests

As follows, we demonstrate the advantages and limitations of the existing models outlined in Sect. 3 and of the novel star-convex model introduced in Sect. 4 through numerical experiments. These are carried out on a bi-axially loaded disk, a plate with hole under compression and two blocks in relative sliding.

Computational solvers for phase-field fracture models based on the finite element method seek a quasi-static solution by solving the time- and space-discretized weak problem. At time step $t$, we look for the solution of the space-discretized weak form of (11) and (12) by using an alternate minimization scheme implemented in FEniCSx Scroggs et al. (2022a, b). In Bourdin et al. (2008), such scheme was introduced in the context of the standard model taking advantage of the separate directional convexity of functions $u \to \mathcal{E}_t(u, \cdot)$ and $\alpha \to \mathcal{E}_t(\cdot, \alpha)$ (see also Ambati et al. (2015); Gerasimov and De Lorenzis (2016) for a detailed discussion on monolithic and alternate minimization schemes). With this scheme, using as initial guess $(u^{k-1}, \alpha^{k-1})$, we first solve the problem in $u$ using a Newton-based nonlinear solver with line search and a maximum of 100
<figure><img src="image_10.png" /><figcaption>Fig. 9 Elastic domains for the star-convex model in the strain (blue dash-dotted line) and stress (red solid line) space for $\alpha = \{0, 0.25, 0.5, 0.75, 1\}$ and $\gamma^* = 5$ ($\nu_0 = 0.3$). The domains for $\alpha = 0$ are the nucleation domains (or strength surfaces). The grey background indicates the parts of the elastic domains shared</figcaption></figure>

with the standard model. Volumetric-deviatoric diagram a: damage cannot evolve for loading directions below the black dashed line in the strain space which is also the constant asymptote. Principal components diagram under the plane-strain assumption b

ing that the damage evolution occurs only above the asymptote, where the deviatoric contribution exceeds the volumetric one. Damage evolution redistributes the stiffness such that the deviatoric stiffness decreases to zero and the volumetric stiffness increases, while the total elastic energy density (which is the sum of the volumetric and deviatoric contributions) always decreases.

Defining as usual the strength surface as the boundary of the elastic domain in the stress space for $\alpha = \alpha_{peak} = 0$, i.e. $\partial \mathcal{R}^*(0)$, we can derive the tensile strength, the shear strength and the compressive strength as

$$ \sigma_e^+ = \sqrt{\frac{2 E_0 w_1 w'(0)}{s'(0)}}, \quad \tau_e = \sqrt{\frac{2 \mu_0 w_1 w'(0)}{s'(0)}}, \\ \sigma_e^- = \begin{cases} -\sqrt{\frac{1}{1+(1+\gamma^*)^{\frac{\alpha-1}{n}}} \frac{2 E_0 w_1 w'(0)}{s'(0)}}, & \text{for } \gamma^* < \frac{n(n-1)}{2} \frac{\kappa_0}{\mu_0}, \\ -\infty, & \text{otherwise.} \end{cases} \quad (41) $$

From (41) we see that the shear strength $\tau_e$ is independent of $\gamma^*$. Thus, the flexibility offered by the star-convex model is exclusively related to the compressive strength $\sigma_e^-$.

## 5 Numerical tests

As follows, we demonstrate the advantages and limitations of the existing models outlined in Sect. 3 and of the novel star-convex model introduced in Sect. 4 through numerical experiments. These are carried out on a bi-axially loaded disk, a plate with hole under compression and two blocks in relative sliding.

Computational solvers for phase-field fracture models based on the finite element method seek a quasi-static solution by solving the time- and space-discretized weak problem. At time step $t$, we look for the solution of the space-discretized weak form of (11) and (12) by using an alternate minimization scheme implemented in FEniCSx Scroggs et al. (2022a, b). In Bourdin et al. (2008), such scheme was introduced in the context of the standard model taking advantage of the separate directional convexity of functions $u \to \mathcal{E}_t(u, \cdot)$ and $\alpha \to \mathcal{E}_t(\cdot, \alpha)$ (see also Ambati et al. (2015); Gerasimov and De Lorenzis (2016) for a detailed discussion on monolithic and alternate minimization schemes). With this scheme, using as initial guess $(u^{k-1}, \alpha^{k-1})$, we first solve the problem in $u$ using a Newton-based nonlinear solver with line search and a maximum of 100
<figure><img src="image_11.png" /><figcaption>**Fig. 10** Three-dimensional nucleation domain in the stress space for the star-convex model with $\gamma^* = 1$ **a** and $\gamma^* = 5$ **b** ($v_0 = 0.3$)</figcaption></figure>

iterations, which gives us the output $\mathbf{u}^k$. Then, starting from $(\mathbf{u}^k, \alpha^{k-1})$, we solve the constrained nonlinear problem in $\alpha$ using a reduced-space active set solver for variational inequalities based on Newton's method with a maximum of 100 iterations, which returns $\alpha^k$. Consequently, we calculate $R_u$, the $L_2$-norm of the residual for the displacement problem for $(\mathbf{u}^k, \alpha^k)$. The convergence of the alternate minimization scheme is achieved if $R_u$ is below the tolerance level ($t \circ 1 = 10^{-6}$).

For $\text{tr}(\boldsymbol{\epsilon}) < 0$, the DP-like model and the star-convex model display a linearly elastic behavior with no damage for $|\boldsymbol{\epsilon}_{\text{dev}}| < -\frac{\gamma^{\kappa_0}}{2\mu_0} \text{tr}(\boldsymbol{\epsilon})$ and $|\boldsymbol{\epsilon}_{\text{dev}}| < -\sqrt{\frac{\kappa_0}{2\mu_0}} \gamma^* \text{tr}(\boldsymbol{\epsilon})$, respectively. By overlapping these two lines, we obtain the relationship $\gamma = \sqrt{2\gamma^*\mu_0/\kappa_0}$ among the parameters $\gamma$ and $\gamma^*$. Accordingly, we adopt $\gamma^* = \{1, 5\}$ for the star-convex model and $\gamma = \{\sqrt{2\mu_0/\kappa_0}, \sqrt{10\mu_0/\kappa_0}\}$ for the DP-like model to facilitate their comparison.

## 5.1 Bi-axially loaded disk test

We first focus on nucleation and consider a bi-axially loaded disk of diameter $D$ under the plane-strain assumption. The geometry and boundary conditions are illustrated in Fig. 11. The center of the disk coincides with the origin of the cartesian reference system $x-y$ and the out-of-plane direction is denoted with $z$. We set $D = 1$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1.5$ and $l = 0.04$. At $t = 0$, the material is assumed to be intact, i.e., $\alpha = 0$ everywhere. The goal of the test is to numerically determine the elastic domains of the various models for $\alpha = 0$ in the $\varepsilon_{xx}-\varepsilon_{yy}$ and $\sigma_{xx}-\sigma_{yy}$ planes.

<figure><img src="image_12.png" /><figcaption>**Fig. 11** Geometry and loading for the bi-axially loaded disk</figcaption></figure>

Defining on the $\varepsilon_{xx}-\varepsilon_{yy}$ plane the angle $\theta \in [0, 2\pi)$, for a given $\theta$ we prescribe a displacement $\bar{\mathbf{u}}_t = t(x \cos(\theta), y \sin(\theta))$ on the disk boundary. In the linear elastic regime, this Dirichlet boundary condition ensures that the disk is subjected to a homogeneous strain along the $\theta$-direction, i.e.,

$$ \frac{\varepsilon_{yy}}{\varepsilon_{xx}} = \tan(\theta). \qquad (42) $$

Maintaining $\theta$ fixed, we increment the value of $t$ until we reach the elastic limit, while always imposing the Dirichlet boundary condition $\alpha = 0$ on the whole boundary of the disk.

We discretize the problem with standard linear triangular finite elements using an unstructured mesh with average size $h = l/3$.

Figures 12 and 13 illustrate the contour of the elastic domain in the strain and stress planes at $\alpha = 0$ for all
<figure><img src="image_11.png" /><figcaption>Fig. 10 Three-dimensional nucleation domain in the stress space for the star-convex model with $\gamma^* = 1$ a and $\gamma^* = 5$ b ($v_0 = 0.3$)</figcaption></figure>

iterations, which gives us the output $\mathbf{u}^k$. Then, starting from $(\mathbf{u}^k, \alpha^{k-1})$, we solve the constrained nonlinear problem in $\alpha$ using a reduced-space active set solver for variational inequalities based on Newton's method with a maximum of 100 iterations, which returns $\alpha^k$. Consequently, we calculate $R_u$, the $L_2$-norm of the residual for the displacement problem for $(\mathbf{u}^k, \alpha^k)$. The convergence of the alternate minimization scheme is achieved if $R_u$ is below the tolerance level ($t \circ 1 = 10^{-6}$).

For $\text{tr}(\boldsymbol{\varepsilon}) < 0$, the DP-like model and the star-convex model display a linearly elastic behavior with no damage for $|\boldsymbol{\varepsilon}_{\text{dev}}| < -\frac{\gamma^{\kappa_0}}{2\mu_0} \text{tr}(\boldsymbol{\varepsilon})$ and $|\boldsymbol{\varepsilon}_{\text{dev}}| < -\sqrt{\frac{\kappa_0}{2\mu_0}} \gamma^* \text{tr}(\boldsymbol{\varepsilon})$, respectively. By overlapping these two lines, we obtain the relationship $\gamma = \sqrt{2\gamma^*\mu_0/\kappa_0}$ among the parameters $\gamma$ and $\gamma^*$. Accordingly, we adopt $\gamma^* = \{1, 5\}$ for the star-convex model and $\gamma = \{\sqrt{2\mu_0/\kappa_0}, \sqrt{10\mu_0/\kappa_0}\}$ for the DP-like model to facilitate their comparison.

## 5.1 Bi-axially loaded disk test

We first focus on nucleation and consider a bi-axially loaded disk of diameter $D$ under the plane-strain assumption. The geometry and boundary conditions are illustrated in Fig. 11. The center of the disk coincides with the origin of the cartesian reference system $x-y$ and the out-of-plane direction is denoted with $z$. We set $D = 1$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1.5$ and $l = 0.04$. At $t = 0$, the material is assumed to be intact, i.e., $\alpha = 0$ everywhere. The goal of the test is to numerically determine the elastic domains of the various models for $\alpha = 0$ in the $\varepsilon_{xx}-\varepsilon_{yy}$ and $\sigma_{xx}-\sigma_{yy}$ planes.

<figure><img src="image_12.png" /><figcaption>Fig. 11 Geometry and loading for the bi-axially loaded disk</figcaption></figure>

Defining on the $\varepsilon_{xx}-\varepsilon_{yy}$ plane the angle $\theta \in [0, 2\pi)$, for a given $\theta$ we prescribe a displacement $\tilde{\mathbf{u}}_t = t(x \cos(\theta), y \sin(\theta))$ on the disk boundary. In the linear elastic regime, this Dirichlet boundary condition ensures that the disk is subjected to a homogeneous strain along the $\theta$-direction, i.e.,

$$ \frac{\varepsilon_{yy}}{\varepsilon_{xx}} = \tan(\theta). \qquad (42) $$

Maintaining $\theta$ fixed, we increment the value of $t$ until we reach the elastic limit, while always imposing the Dirichlet boundary condition $\alpha = 0$ on the whole boundary of the disk.

We discretize the problem with standard linear triangular finite elements using an unstructured mesh with average size $h = l/3$.

Figures 12 and 13 illustrate the contour of the elastic domain in the strain and stress planes at $\alpha = 0$ for all
models. The blue and red lines correspond to the theo-
retical strain and stress domains, respectively, while the
dots of the respective colors represent the numerically
obtained strain and stress states right before damage
nucleation. In some cases, the dots are absent because
certain loading directions do not reach the limit or reach
it after a very large number of time steps (outside the
range of the figure). For each load direction we plot the
damage field at the time step when nucleation occurs, if
it occurs; for directions at which no nucleation occurs
we plot a zero damage field for completeness. Hence,
the test serves not only as a nucleation test but also pro-
vides valuable insights into the possible post-critical
crack patterns under various loading conditions.

Figure 12a presents the simulation results for the standard model, which displays elliptical elastic domains and maintains symmetry of behavior in both tension and compression. Accordingly, the elastic domains are symmetric not only with respect to the bisector of the first and third quadrants but also with respect to the bisector of the second and fourth quadrants.

Figure 12b illustrates the outcomes of the test using the volumetric-deviatoric model. The “tensile” part (tr(e) ≥ 0) behaves as in the standard model, while an asymmetry in the compressive behavior is introduced by the split. Due to the plane-strain assumption, there is no purely volumetric direction; as a result, the elastic domain is closed, leading to damage nucleation for every loading direction. In the third quadrant of Fig. 12b nucleation occurs in all cases, although in three of them the damage level is so low to be hardly visible. Understanding why these cases nucleate but do not localize into a crack would require a second-order stability analysis for multi-axial stress states Pham and Marigo (2013); Baldelli and Maurini (2021), a task going beyond the scope of the present work and currently in progress (Zolesi and Maurini 2023).

Figure 12c and 12d provide the results for the spec-
tral and the no-tension models, respectively. In these
cases, the domains are open since for some directions
(some values of θ) nucleation cannot occur. The extent
of the no-nucleation region varies for the spectral and
the no-tension models but cannot be adjusted based
on experimental results. As discussed in Sect. 3.2, the
evolution of the elastic domain in the stress plane for
the spectral split is a complex non-homothetic trans-
formation with respect to the origin that does not enjoy
the stress-softening property. This can explain the pres-

ence in Fig. 12c of some localized damage fields not
reaching $\alpha = 1$.

Figure 12e and 12f depict the results of the test using the DP-like model with $\gamma = \sqrt{2\mu_0/\kappa_0}$ and $\gamma = \sqrt{10\mu_0/\kappa_0}$, respectively. Similar to the spectral and no-tension models, the domains remain open, hence for some values of $\theta$ no nucleation occurs. However, the opening angle of the no-nucleation region can now be adjusted using the additional parameter $\gamma$, namely, it becomes wider with increasing $\gamma$. At the same time the region of $\partial\mathcal{R}(\alpha)$ and $\partial\mathcal{R}^*(\alpha)$ shared with the standard model becomes narrower as $\gamma$ increases. While on this shared region the nucleation behavior is the same as in the standard model (the dots have the same locations), the propagation behavior (which can be appreciated looking at the plotted localized damage fields after nucleation) is different. For instance, for $\gamma = \sqrt{10\mu_0/\kappa_0}$ in Fig. 12f, the cracks exhibit abnormal thickening with respect to those in Fig. 12a even in the first quadrant where the shared portion lies. Hence, while this model exhibits the desired flexibility in nucleation, it leads to a problematic behavior during propagation related to the non-zero residual stresses for $\text{tr}(e) > 0$.

Finally, Fig. 13a and 13b illustrate the results for the star-convex split with $\gamma^* = \{1, 5\}$. Similar to the DP-like model, the star-convex split also exhibits a certain flexibility in the opening angle of the no-nucleation region. As shown in the figures, for increasing values of $\gamma^*$, the opening angle widens. At the same time, the star-convex split shows a better propagation behavior, without the abnormal thickening of cracks observed with the DP-like model. In fact, the residual stresses are entirely absent for $\text{tr}(e) \ge 0$ and proportional to the negative volumetric strain.

5.2 Plate with hole test

We now discuss the compression of a plate with hole
under plane-strain conditions. Assuming the problem
to be symmetric, the test is performed on a quarter of
the plate as shown in Fig. 14.

Chao Correas et al. (2023) performed recently a sim-
ilar nucleation test under multiaxial loading, comparing
the standard model and the no-tension model to the pre-
diction of cohesive-zone models and the finite-fracture
mechanics approach.
models. The blue and red lines correspond to the theo-
retical strain and stress domains, respectively, while the
dots of the respective colors represent the numerically
obtained strain and stress states right before damage
nucleation. In some cases, the dots are absent because
certain loading directions do not reach the limit or reach
it after a very large number of time steps (outside the
range of the figure). For each load direction we plot the
damage field at the time step when nucleation occurs, if
it occurs; for directions at which no nucleation occurs
we plot a zero damage field for completeness. Hence,
the test serves not only as a nucleation test but also pro-
vides valuable insights into the possible post-critical
crack patterns under various loading conditions.

Figure 12a presents the simulation results for the standard model, which displays elliptical elastic domains and maintains symmetry of behavior in both tension and compression. Accordingly, the elastic domains are symmetric not only with respect to the bisector of the first and third quadrants but also with respect to the bisector of the second and fourth quadrants.

Figure 12b illustrates the outcomes of the test using the volumetric-deviatoric model. The “tensile” part (tr(e) ≥ 0) behaves as in the standard model, while an asymmetry in the compressive behavior is introduced by the split. Due to the plane-strain assumption, there is no purely volumetric direction; as a result, the elastic domain is closed, leading to damage nucleation for every loading direction. In the third quadrant of Fig. 12b nucleation occurs in all cases, although in three of them the damage level is so low to be hardly visible. Understanding why these cases nucleate but do not localize into a crack would require a second-order stability analysis for multi-axial stress states Pham and Marigo (2013); Baldelli and Maurini (2021), a task going beyond the scope of the present work and currently in progress (Zolesi and Maurini 2023).

Figure 12c and 12d provide the results for the spec-
tral and the no-tension models, respectively. In these
cases, the domains are open since for some directions
(some values of θ) nucleation cannot occur. The extent
of the no-nucleation region varies for the spectral and
the no-tension models but cannot be adjusted based
on experimental results. As discussed in Sect. 3.2, the
evolution of the elastic domain in the stress plane for
the spectral split is a complex non-homothetic trans-
formation with respect to the origin that does not enjoy
the stress-softening property. This can explain the pres-

ence in Fig. 12c of some localized damage fields not
reaching $\alpha = 1$.

Figure 12e and 12f depict the results of the test using the DP-like model with $\gamma = \sqrt{2\mu_0/\kappa_0}$ and $\gamma = \sqrt{10\mu_0/\kappa_0}$, respectively. Similar to the spectral and no-tension models, the domains remain open, hence for some values of $\theta$ no nucleation occurs. However, the opening angle of the no-nucleation region can now be adjusted using the additional parameter $\gamma$, namely, it becomes wider with increasing $\gamma$. At the same time the region of $\partial\mathcal{R}(\alpha)$ and $\partial\mathcal{R}^*(\alpha)$ shared with the standard model becomes narrower as $\gamma$ increases. While on this shared region the nucleation behavior is the same as in the standard model (the dots have the same locations), the propagation behavior (which can be appreciated looking at the plotted localized damage fields after nucleation) is different. For instance, for $\gamma = \sqrt{10\mu_0/\kappa_0}$ in Fig. 12f, the cracks exhibit abnormal thickening with respect to those in Fig. 12a even in the first quadrant where the shared portion lies. Hence, while this model exhibits the desired flexibility in nucleation, it leads to a problematic behavior during propagation related to the non-zero residual stresses for $\text{tr}(e) > 0$.

Finally, Fig. 13a and 13b illustrate the results for the star-convex split with $\gamma^* = \{1, 5\}$. Similar to the DP-like model, the star-convex split also exhibits a certain flexibility in the opening angle of the no-nucleation region. As shown in the figures, for increasing values of $\gamma^*$, the opening angle widens. At the same time, the star-convex split shows a better propagation behavior, without the abnormal thickening of cracks observed with the DP-like model. In fact, the residual stresses are entirely absent for $\text{tr}(e) \ge 0$ and proportional to the negative volumetric strain.

5.2 Plate with hole test

We now discuss the compression of a plate with hole
under plane-strain conditions. Assuming the problem
to be symmetric, the test is performed on a quarter of
the plate as shown in Fig. 14.

Chao Correas et al. (2023) performed recently a sim-
ilar nucleation test under multiaxial loading, comparing
the standard model and the no-tension model to the pre-
diction of cohesive-zone models and the finite-fracture
mechanics approach.
<figure><img src="image_13.png" /><figcaption>**Fig. 12** Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the available models</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 12 Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the available models</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the star-convex model</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Geometry and loading for the plate with hole</figcaption></figure>

The plate is a square with side length $2L$ and a central hole with radius $R$. The hole center coincides with the origin of both the cartesian reference system $x-y$ and the polar one $r-\theta$. The out-of-plane direction is denoted with $z$. Rigid body motions are prevented by the symmetry constraints $u_y(x, y = 0) = 0$ and $u_x(x = 0, y) = 0$ and compression is prescribed through the upper vertical displacement $u_y(x, y = L) = -tL$ with $t$ going from 0 up to 0.3. At $t = 0$, the material is assumed to be intact, i. e., $\alpha = 0$ everywhere. We adopt $L = 1$, $R = 0.3$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1$ and $\ell = 0.02$. We discretize the problem with standard bilinear quadrilateral finite elements using an unstructured mesh with average size $h = \ell/3$.

During the elastic regime, the highest stress concentration is expected at the boundary of the hole ($r = R$), along which the radial stress $\sigma_{rr}$ and the tangential stress $\sigma_{r\theta}$ are both zero ($\sigma_{rr} = \sigma_{r\theta} = 0$) and the hoop stress $\sigma_{\theta\theta}$ and the out-of-plane stress $\sigma_{zz}$ are related by the plane-strain condition $\sigma_{zz} = -v_0\sigma_{\theta\theta}$. Along the boundary of the hole, point A ($r = R, \theta = \pi/2$) undergoes the maximum tensile hoop stress $\sigma_{\theta\theta}^A \ge 0$ whereas point B ($r = R, \theta = 0$) is subjected to the lowest compressive hoop stress $\sigma_{\theta\theta}^B \le -\sigma_{\theta\theta}^A$.

We now wonder whether the crack will originate from A or B. Experiments with this setup for many brittle materials produce a so-called “splitting” crack that starts at point A and proceeds parallel to the loading direction upwards (Romani et al. 2015). On the other hand, nucleation as predicted by a phase-field model will start from the point at which the stress state first touches the limit curve $\partial R^*(0)$. For first nucleation in A, the stress state in A must reach the intersection between $\partial R^*(0)$ and the half-line $\sigma_{zz} = -v_0\sigma_{\theta\theta} \le 0$. We denote with $\sigma_{\theta\theta_e}^+$ the value of the hoop stress at this intersection point (this concept is exemplified in Fig. 15 for the star-convex model). Similarly, for first nucleation in B, the stress state in B must reach the intersection between $\partial R^*(0)$ and the half-line $\sigma_{zz} = -v_0\sigma_{\theta\theta} \ge 0$. We denote with $\sigma_{\theta\theta_e}^-$ the value of the hoop stress at this intersection point (Fig. 15). The $|\sigma_{\theta\theta_e}^- / \sigma_{\theta\theta_e}^+$ ratio increases with the compressive to tensile strength ratio $|\sigma_e^- / \sigma_e^+|$ with a trend which is specific to each model. Since when $L \gg R$, $\sigma_{\theta\theta}^B = -3\sigma_{\theta\theta}^A$ (Kirsch 1898), we
<figure><img src="image_14.png" /><figcaption>**Fig. 13** Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the star-convex model</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 14** Geometry and loading for the plate with hole</figcaption></figure>

The plate is a square with side length $2L$ and a central hole with radius $R$. The hole center coincides with the origin of both the cartesian reference system $x-y$ and the polar one $r-\theta$. The out-of-plane direction is denoted with $z$. Rigid body motions are prevented by the symmetry constraints $u_y(x, y = 0) = 0$ and $u_x(x = 0, y) = 0$ and compression is prescribed through the upper vertical displacement $u_y(x, y = L) = -tL$ with $t$ going from 0 up to 0.3. At $t = 0$, the material is assumed to be intact, i. e., $\alpha = 0$ everywhere. We adopt $L = 1$, $R = 0.3$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1$ and $\ell = 0.02$. We discretize the problem with standard bilinear quadrilateral finite elements using an unstructured mesh with average size $h = \ell/3$.

During the elastic regime, the highest stress concentration is expected at the boundary of the hole ($r = R$), along which the radial stress $\sigma_{rr}$ and the tangential stress $\sigma_{r\theta}$ are both zero ($\sigma_{rr} = \sigma_{r\theta} = 0$) and the hoop stress $\sigma_{\theta\theta}$ and the out-of-plane stress $\sigma_{zz}$ are related by the plane-strain condition $\sigma_{zz} = -v_0\sigma_{\theta\theta}$. Along the boundary of the hole, point A ($r = R, \theta = \pi/2$) undergoes the maximum tensile hoop stress $\sigma_{\theta\theta}^A \ge 0$ whereas point B ($r = R, \theta = 0$) is subjected to the lowest compressive hoop stress $\sigma_{\theta\theta}^B \le -\sigma_{\theta\theta}^A$.

We now wonder whether the crack will originate from A or B. Experiments with this setup for many brittle materials produce a so-called “splitting” crack that starts at point A and proceeds parallel to the loading direction upwards (Romani et al. 2015). On the other hand, nucleation as predicted by a phase-field model will start from the point at which the stress state first touches the limit curve $\partial R^*(0)$. For first nucleation in A, the stress state in A must reach the intersection between $\partial R^*(0)$ and the half-line $\sigma_{zz} = -v_0\sigma_{\theta\theta} \le 0$. We denote with $\sigma_{\theta\theta_e}^+$ the value of the hoop stress at this intersection point (this concept is exemplified in Fig. 15 for the star-convex model). Similarly, for first nucleation in B, the stress state in B must reach the intersection between $\partial R^*(0)$ and the half-line $\sigma_{zz} = -v_0\sigma_{\theta\theta} \ge 0$. We denote with $\sigma_{\theta\theta_e}^-$ the value of the hoop stress at this intersection point (Fig. 15). The $|\sigma_{\theta\theta_e}^- / \sigma_{\theta\theta_e}^+$ ratio increases with the compressive to tensile strength ratio $|\sigma_e^- / \sigma_e^+|$ with a trend which is specific to each model. Since when $L \gg R$, $\sigma_{\theta\theta}^B = -3\sigma_{\theta\theta}^A$ (Kirsch 1898), we
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Effect of $\gamma^*$ on the nucleation domain in the plane spanned by hoop stress and out-of-plane stress for the plate with hole test. Elastic limit hoop stresses (a): the red solid line and the black dashed line represent the contour $\partial\mathcal{R}^*(0)$ for $\gamma^* = \{-1, 0, 1, 2, 5\}$ and the stress state under the plane-strain</figcaption></figure>

assumption, respectively. Their intersection determines the compressive limit $\sigma_{\bar{\theta}\theta_e}^-$ (blue on the left) and the tensile limit $\sigma_{\theta\theta_e}^+$ (blue on the right). Compressive-tensile hoop limits ratio (b): trend of the limits ratio vs. $\gamma^*$ (red solid line). For $\gamma^* > \gamma_3^*$ the ratio is greater than 3, for $\gamma^* \ge \gamma_\infty^*$ the ratio is equal to $\infty$

<table><caption>**Table 3** $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio for the plate with hole. With reference to Fig. 14, first nucleation is expected in B when the ratio is well below 3. Nucleation is expected in A when the ratio is well above 3</caption><thead><tr><th>Model</th><th>|σ<sub>θθ<sub>e</sub></sub><sup>-</sup>/σ<sub>θθ<sub>e</sub></sub><sup>+</sup>|</th><th>v<sub>0</sub> = 0.3</th></tr></thead><tbody><tr><td>Standard</td><td>1</td><td>= 1</td></tr><tr><td>Vol.-dev.</td><td>√(3(1-v<sub>0</sub><sup>2</sup>)/(2(1+v<sub>0</sub><sup>3</sup>)))</td><td>≈ 1.15</td></tr><tr><td>Spectral</td><td>√(1-v<sub>0</sub><sup>2</sup>)/(v<sub>0</sub>)</td><td>≈ 2.60</td></tr><tr><td>No-tension</td><td>∞</td><td>= ∞</td></tr><tr><td>DP-like</td><td><div style="display: flex; align-items: center; justify-content: space-around; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">√(3(1-v<sub>0</sub>)(γ<sub>∞</sub>ζ + γ<sup>2</sup>))</span><span style="text-align: center; line-height: 1.2;">√(1-2v<sub>0</sub>)(1+v<sub>0</sub>)(γ<sub>∞</sub>-γ)</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub>+γ</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub>-γ</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">∞</span><span style="text-align: center; line-height: 1.2;">∞</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">for</span><span style="text-align: center; line-height: 1.2;">γ</span><span style="text-align: center; line-height: 1.2;">&lt;</span><span style="text-align: center; line-height: 1.2;">ζ</span><span style="text-align: center; line-height: 1.2;">= (3(1-2v<sub>0</sub>)/(γ<sub>∞</sub>(1+v<sub>0</sub>)))</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">for</span><span style="text-align: center; line-height: 1.2;">ζ</span><span style="text-align: center; line-height: 1.2;">&le;</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub></span><span style="text-align: center; line-height: 1.2;">= √(6(1-v<sub>0</sub>+v<sub>0</sub><sup>2</sup>)/(1+v<sub>0</sub>))</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">otherwise</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub></span><span style="text-align: center; line-height: 1.2;">&approx;</span><span style="text-align: center; line-height: 1.2;">1.67</span></div></td><td><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">ζ</span><span style="text-align: center; line-height: 1.2;">&approx;</span><span style="text-align: center; line-height: 1.2;">0.55</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub></span><span style="text-align: center; line-height: 1.2;">&approx;</span><span style="text-align: center; line-height: 1.2;">1.67</span></div></td></tr><tr><td>Star-convex</td><td><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">√(γ<sub>∞</sub><sup>2</sup>+1)</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub><sup>2</sup>-γ<sup>*</sup></span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">∞</span><span style="text-align: center; line-height: 1.2;">∞</span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">for</span><span style="text-align: center; line-height: 1.2;">γ<sup>*</sup></span><span style="text-align: center; line-height: 1.2;">&lt;</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub><sup>*</sup></span></div><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">otherwise</span><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub><sup>*</sup></span><span style="text-align: center; line-height: 1.2;">&approx;</span><span style="text-align: center; line-height: 1.2;">3.04</span></div></td><td><div style="display: flex; align-items: center; justify-content: space-between; margin-bottom: 5px;"><span style="text-align: center; line-height: 1.2;">γ<sub>∞</sub><sup>*</sup></span><span style="text-align: center; line-height: 1.2;">&approx;</span><span style="text-align: center; line-height: 1.2;">3.04</span></div></td></tr></tbody></table>

can infer that if the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is well below 3, nucleation will occur in B, and if it is well above 3, it will occur in A (Table 3).

For the standard phase-field model, which has an elastic domain symmetric in tension and compression, damage starts at point B, as shown in Fig. 16b. Additionally, the asymmetry alone is not sufficient to prevent crack initiation in B, as visible in Figs. 16d, 16f and 17f as for these models the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is well below 3 (Table 3). Unlike the other models, the DP-like and the star-convex models contain the additional parameters $\gamma$ and $\gamma^*$ that allow for the flexible calibration of the

strength ratio $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$. For the star-convex model, Fig. 15 shows how the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is affected by the variation of $\gamma^*$.

Figures 16 and 17 present the results of the plate with hole test. For each model, we show the plot of the vertical stress reaction per unit thickness $R_y$, equal to the integral of $\sigma_{yy}$ over the upper boundary, against the pseudo-time $t$, as well as the damage field at the last time step. Some models do not reach iterative convergence of the numerical solver for all time steps; for these models, we report the results up to the last converged time step (a brief discussion on iterative con-
<figure><img src="image_16.png" /><figcaption>**Fig. 15** Effect of $\gamma^*$ on the nucleation domain in the plane spanned by hoop stress and out-of-plane stress for the plate with hole test. Elastic limit hoop stresses (a): the red solid line and the black dashed line represent the contour $\partial\mathcal{R}^*(0)$ for $\gamma^* = \{-1, 0, 1, 2, 5\}$ and the stress state under the plane-strain</figcaption></figure>

assumption, respectively. Their intersection determines the compressive limit $\sigma_{\bar{\theta}\theta_e}^-$ (blue on the left) and the tensile limit $\sigma_{\theta\theta_e}^+$ (blue on the right). Compressive-tensile hoop limits ratio (b): trend of the limits ratio vs. $\gamma^*$ (red solid line). For $\gamma^* > \gamma_3^*$ the ratio is greater than 3, for $\gamma^* \ge \gamma_\infty^*$ the ratio is equal to $\infty$

<table><caption>**Table 3** $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio for the plate with hole. With reference to Fig. 14, first nucleation is expected in B when the ratio is well below 3. Nucleation is expected in A when the ratio is well above 3</caption><thead><tr><th>Model</th><th>|σ<sub>θθ<sub>e</sub></sub><sup>-</sup>/σ<sub>θθ<sub>e</sub></sub><sup>+</sup>|</th><th>v<sub>0</sub> = 0.3</th></tr></thead><tbody><tr><td>Standard</td><td>1</td><td>= 1</td></tr><tr><td>Vol.-dev.</td><td>√(3(1-v<sub>0</sub><sup>2</sup>)/(2(1+v<sub>0</sub><sup>3</sup>)))</td><td>≈ 1.15</td></tr><tr><td>Spectral</td><td>√(1-v<sub>0</sub><sup>2</sup>) / v<sub>0</sub></td><td>≈ 2.60</td></tr><tr><td>No-tension</td><td>∞</td><td>= ∞</td></tr><tr><td>DP-like</td><td>√(3(1-v<sub>0</sub>)(γ<sub>∞</sub>ζ + γ<sup>2</sup>)) / (√(1-2v<sub>0</sub>)(1+v<sub>0</sub>)(γ<sub>∞</sub>-γ))</td><td>ζ ≈ 0.55</td></tr><tr><td>Star-convex</td><td>√(γ<sub>∞</sub><sup>2</sup>+1) / (γ<sub>∞</sub><sup>2</sup>-γ<sup>*</sup>)</td><td>γ<sub>∞</sub> ≈ 1.67</td></tr></tbody></table>

can infer that if the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is well below 3, nucleation will occur in B, and if it is well above 3, it will occur in A (Table 3).

For the standard phase-field model, which has an elastic domain symmetric in tension and compression, damage starts at point B, as shown in Fig. 16b. Additionally, the asymmetry alone is not sufficient to prevent crack initiation in B, as visible in Figs. 16d, 16f and 17f as for these models the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is well below 3 (Table 3). Unlike the other models, the DP-like and the star-convex models contain the additional parameters $\gamma$ and $\gamma^*$ that allow for the flexible calibration of the

strength ratio $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$. For the star-convex model, Fig. 15 shows how the $|\sigma_{\bar{\theta}\theta_e}/\sigma_{\theta\theta_e}^+|$ ratio is affected by the variation of $\gamma^*$.

Figures 16 and 17 present the results of the plate with hole test. For each model, we show the plot of the vertical stress reaction per unit thickness $R_y$, equal to the integral of $\sigma_{yy}$ over the upper boundary, against the pseudo-time $t$, as well as the damage field at the last time step. Some models do not reach iterative convergence of the numerical solver for all time steps; for these models, we report the results up to the last converged time step (a brief discussion on iterative con-
<figure><img src="image_17.png" /><figcaption>**Fig. 16** Results of the plate with hole test for the standard, volumetric-deviatoric, spectral and no-tension models. Vertical reaction against pseudo-time *t* a, c, e, g: the black dash-dotted line corresponds to the elastic phase, the *e* point to damage nucleation and the red solid line to the post-nucleation phase. Damage fields **b**, **d**, **f**, **h**</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Fig. 16** Results of the plate with hole test for the standard, volumetric-deviatoric, spectral and no-tension models. Vertical reaction against pseudo-time *t* a, c, e, g: the black dash-dotted line corresponds to the elastic phase, the *e* point to damage nucleation and the red solid line to the post-nucleation phase. Damage fields **b**, **d**, **f**, **h**</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 17** Results of the plate with hole test for the DP-like ($\gamma = \{\sqrt{2\mu_0/\kappa_0}, \sqrt{10\mu_0/\kappa_0}\}$) and star-convex models ($\gamma^* = \{1, 5\}$). Vertical reaction against pseudo-time $t$ **a**, **c**, **e**, **g**: the black dash-dotted line corresponds to the elastic phase, the **e** point to damage nucleation and the red solid line to the post-nucleation phase. Damage fields **b**, **d**, **f**, **h**</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 17** Results of the plate with hole test for the DP-like ($\gamma = \{\sqrt{2\mu_0}/\kappa_0, \sqrt{10\mu_0}/\kappa_0\}$) and star-convex models ($\gamma^* = \{1, 5\}$). Vertical reaction against pseudo-time $t$ **a**, **c**, **e**, **g**: the black dash-dotted line corresponds to the elastic phase, the **e** point to damage nucleation and the red solid line to the post-nucleation phase. Damage fields **b**, **d**, **f**, **h**</figcaption></figure>
vergence issues is postponed to Sect. 5.4). In the $R_y$ vs. $t$ plots, the elastic limit point $e$ marks the transition between the elastic regime (black dash-dotted line) and the stage after damage nucleation (red solid line). From Figs. 16a, c, e and 17e we observe that, in models for which nucleation occurs at point B, the elastic limit $e$ is met at different time steps. While the volumetric-deviatoric, spectral, and star-convex split with $\gamma^* = 1$ introduce tension/compression asymmetry, this is not sufficient to obtain damage onset at point A. In the case of the spectral split in Fig. 16f, the nucleation of the diffuse damage in B occurs before vertical localization at point A. Consequently, the nucleation point $e$ in Fig. 16e corresponds to the initiation of damage in B. Thus, the only models that exhibit a crack propagation consistent with the experimental evidence Romani et al. (2015) are those with a certain level of flexibility, adjustable with the parameter $\gamma$ (DP-like model) or $\gamma^*$ (star-convex split), along with the no-tension model. However, in the case of the DP-like model (Fig. 17d), for $\gamma = \sqrt{10\mu_0/\kappa_0}$, the damage exhibits a diffuse pattern rather than being fully localized. This is attributed to the phenomenon of crack thickening as $\gamma$ increases, as previously highlighted in the bi-axially loaded disk test.

## 5.3 Sliding test

The last test concerns a square plate of side $L$ with a crack in the center, as shown in Fig. 18, under plane-strain conditions. The bottom left corner of the plate coincides with the origin of the cartesian reference system $x-y$ and the out-of-plane direction is denoted with $z$. In order to emulate a crack using a localized damage field, we impose the Dirichlet boundary condition $\alpha(x, y = L/2) = 1$ and we solve the minimization problem (6) while maintaining the displacement field homogeneously equal to zero. Starting from the obtained damage field, we apply the Dirichlet boundary conditions in Fig. 18 to the displacement field: the lower part of the plate is blocked in both directions and a horizontal displacement $u_x(x, y = L) = tL$ is applied to the upper part. Furthermore, the two sides are vertically blocked in their lower part $u_y(x = 0, y \le L/2) = u_y(x = L, y \le L/2) = 0$ and we apply a constant (i.e., independent of $t$) vertical displacement on the upper part $u_y(x = 0, y \ge L/2) = u_y(x = L, y \ge L/2) = \bar{u}L$, where $\bar{u} = 0.1$. This test aims to emphasize the

<figure><img src="image_19.png" /><figcaption>Fig. 18 Geometry and loading for the sliding test</figcaption></figure>

importance of having zero residual shear stresses at full damage. Indeed, if these stresses are not zero, even though the two faces of the block are separate, physically unrealistic stresses are transmitted during sliding. This test is inspired by the one illustrated in Strobl and Seelig (2016) with additional separation of the two blocks. This excludes the possibility that the transmitted stresses are attributable to friction. We set $L = 1$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1$ and $\ell = 0.05$. We discretize the problem with standard bilinear quadrilateral finite elements using a structured mesh with $h = \ell/5$.

Figure 19 depicts the deformed configuration at $t = 0.2$ for each model, except for the spectral and no-tension models, both of which experience convergence issues and thus stop at $t = 0.06$. The no-tension model and the DP-like model for $\gamma = \sqrt{10\mu_0/\kappa_0}$ in Fig. 19d and f, respectively, display unrealistic deformations. Figure 20 illustrates for each model the horizontal stress reaction per unit thickness $R_x$, given by the integral of $\sigma_{xy}$ over the upper side of the plate, against $t$, compared to the horizontal elastic stress reaction of an undamaged plate (green dotted line). Given the presence of the crack and the initial separation, we expect the horizontal reaction force to be identically zero. However, this occurs only for the standard model, which has no residual stresses, and for the volumetric-deviatoric and the star-convex models, both of which have residual stresses proportional to $\langle\text{tr}(\boldsymbol{e})\rangle$. The spectral and the no-tension models, as well as the DP-like model for $\gamma = \sqrt{10\mu_0/\kappa_0}$, exhibit a significant reaction force, slightly lower than half of that observed in an undamaged material. The DP-like model for $\gamma = \sqrt{2\mu_0/\kappa_0}$ initially displays no reaction and a deformation as shown in Fig. 19, which is acceptable. However, it subsequently shows an increase in the
vergence issues is postponed to Sect. 5.4). In the $R_y$ vs. $t$ plots, the elastic limit point $e$ marks the transition between the elastic regime (black dash-dotted line) and the stage after damage nucleation (red solid line). From Figs. 16a, c, e and 17e we observe that, in models for which nucleation occurs at point B, the elastic limit $e$ is met at different time steps. While the volumetric-deviatoric, spectral, and star-convex split with $\gamma^* = 1$ introduce tension/compression asymmetry, this is not sufficient to obtain damage onset at point A. In the case of the spectral split in Fig. 16f, the nucleation of the diffuse damage in B occurs before vertical localization at point A. Consequently, the nucleation point $e$ in Fig. 16e corresponds to the initiation of damage in B. Thus, the only models that exhibit a crack propagation consistent with the experimental evidence Romani et al. (2015) are those with a certain level of flexibility, adjustable with the parameter $\gamma$ (DP-like model) or $\gamma^*$ (star-convex split), along with the no-tension model. However, in the case of the DP-like model (Fig. 17d), for $\gamma = \sqrt{10\mu_0/\kappa_0}$, the damage exhibits a diffuse pattern rather than being fully localized. This is attributed to the phenomenon of crack thickening as $\gamma$ increases, as previously highlighted in the bi-axially loaded disk test.

## 5.3 Sliding test

The last test concerns a square plate of side $L$ with a crack in the center, as shown in Fig. 18, under plane-strain conditions. The bottom left corner of the plate coincides with the origin of the cartesian reference system x-y and the out-of-plane direction is denoted with $z$. In order to emulate a crack using a localized damage field, we impose the Dirichlet boundary condition $\alpha(x, y = L/2) = 1$ and we solve the minimization problem (6) while maintaining the displacement field homogeneously equal to zero. Starting from the obtained damage field, we apply the Dirichlet boundary conditions in Fig. 18 to the displacement field: the lower part of the plate is blocked in both directions and a horizontal displacement $u_x(x, y = L) = tL$ is applied to the upper part. Furthermore, the two sides are vertically blocked in their lower part $u_y(x = 0, y \le L/2) = u_y(x = L, y \le L/2) = 0$ and we apply a constant (i.e., independent of $t$) vertical displacement on the upper part $u_y(x = 0, y \ge L/2) = u_y(x = L, y \ge L/2) = \bar{u}L$, where $\bar{u} = 0.1$. This test aims to emphasize the

importance of having zero residual shear stresses at full damage. Indeed, if these stresses are not zero, even though the two faces of the block are separate, physically unrealistic stresses are transmitted during sliding. This test is inspired by the one illustrated in Strobl and Seelig (2016) with additional separation of the two blocks. This excludes the possibility that the transmitted stresses are attributable to friction. We set $L = 1$, $E_0 = 100$, $v_0 = 0.3$, $w_1 = 1$ and $\ell = 0.05$. We discretize the problem with standard bilinear quadrilateral finite elements using a structured mesh with $h = \ell/5$.

Figure 19 depicts the deformed configuration at $t = 0.2$ for each model, except for the spectral and no-tension models, both of which experience convergence issues and thus stop at $t = 0.06$. The no-tension model and the DP-like model for $\gamma = \sqrt{10\mu_0/\kappa_0}$ in Fig. 19d and f, respectively, display unrealistic deformations. Figure 20 illustrates for each model the horizontal stress reaction per unit thickness $R_x$, given by the integral of $\sigma_{xy}$ over the upper side of the plate, against $t$, compared to the horizontal elastic stress reaction of an undamaged plate (green dotted line). Given the presence of the crack and the initial separation, we expect the horizontal reaction force to be identically zero. However, this occurs only for the standard model, which has no residual stresses, and for the volumetric-deviatoric and the star-convex models, both of which have residual stresses proportional to $\langle\text{tr}(\boldsymbol{e})\rangle$. The spectral and the no-tension models, as well as the DP-like model for $\gamma = \sqrt{10\mu_0/\kappa_0}$, exhibit a significant reaction force, slightly lower than half of that observed in an undamaged material. The DP-like model for $\gamma = \sqrt{2\mu_0/\kappa_0}$ initially displays no reaction and a deformation as shown in Fig. 19, which is acceptable. However, it subsequently shows an increase in the

<figure><img src="image_19.png" /><figcaption>Fig. 18 Geometry and loading for the sliding test</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19 Damage fields on the deformed configurations for all models</figcaption></figure>

reaction force. This is because a lower $\gamma$ value delays the undesired residual stress effect.

## 5.4 Iterative convergence issues

As described in Sect. 5, the iterative convergence criterion of the alternate minimization is met when $R_u < tol$. Given the convexity of the functional in the separate problems of displacement and damage, the standard model is robust in terms of iterative convergence. However, due to the high non-linearity introduced by the energy decompositions, for the other models iterative convergence is often problematic. In fact, we already showed in Sects. 5.2 and 5.3 that, for each model, we had to stop specific numerical tests at the first non-converging time step. The presented results exclusively relate to converged simulations. On the other hand, the same tests with modified parameters can lead to no convergence. In Fig. 21, we illustrate the exemplary case of an angle $\theta$ in the bi-axially loaded disk test using the star-convex model, for which iterative convergence is lost when varying the mesh size $h$. This issue calls for a thorough numerical investigation addressing the specific non-linearities introduced by the splits. It is a crucial issue to be resolved for the proper numerical implementation of phase-field modeling of fracture under multi-axial stress states, but it goes beyond the scope of this work.

<figure><img src="image_21.png" /><figcaption>Fig. 20 Horizontal reaction vs. pseudo-time $t$ obtained from different models: linear elasticity (green dotted line) and phase-field models with spectral split (black solid line with empty circles), no-tension model (purple dashed line), DP-like model with $\gamma = \sqrt{2\mu_0/\kappa_0}$ (red dash-dotted line), DP-like model with $\gamma = \sqrt{10\mu_0/\kappa_0}$ (blue solid line with point marker). The grey solid line corresponds to the star-convex split with $\gamma^* = \{1, 5\}$, to the volumetric-deviatoric split ($\gamma^* = 0$) and to the standard model ($\gamma^* = -1$)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 21 Bi-axially loaded disk with the star-convex model ($\gamma^* = 1$, $\theta = \frac{3}{4}\pi$): $R_u$ against iterations $k$ of the alternate minimization scheme for the time step at which the elastic limit is reached. All the parameters are as in Sect. 5.1 for the solid blue line ($h = \frac{\ell}{3}$). For the solid red line, only the mesh size differs ($h = \frac{\ell}{5}$)</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 19 Damage fields on the deformed configurations for all models</figcaption></figure>

reaction force. This is because a lower $\gamma$ value delays the undesired residual stress effect.

## 5.4 Iterative convergence issues

As described in Sect. 5, the iterative convergence criterion of the alternate minimization is met when $R_u < tol$. Given the convexity of the functional in the separate problems of displacement and damage, the standard model is robust in terms of iterative convergence. However, due to the high non-linearity introduced by the energy decompositions, for the other models iterative convergence is often problematic. In fact, we already showed in Sects. 5.2 and 5.3 that, for each model, we had to stop specific numerical tests at the first non-converging time step. The presented results exclusively relate to converged simulations. On the other hand, the same tests with modified parameters can lead to no convergence. In Fig. 21, we illustrate the exemplary case of an angle $\theta$ in the bi-axially loaded disk test using the star-convex model, for which iterative convergence is lost when varying the mesh size $h$. This issue calls for a thorough numerical investigation addressing the specific non-linearities introduced by the splits. It is a crucial issue to be resolved for the proper numerical implementation of phase-field modeling of fracture under multi-axial stress states, but it goes beyond the scope of this work.

<figure><img src="image_21.png" /><figcaption>Fig. 20 Horizontal reaction vs. pseudo-time $t$ obtained from different models: linear elasticity (green dotted line) and phase-field models with spectral split (black solid line with empty circles), no-tension model (purple dashed line), DP-like model with $\gamma = \sqrt{2\mu_0/\kappa_0}$ (red dash-dotted line), DP-like model with $\gamma = \sqrt{10\mu_0/\kappa_0}$ (blue solid line with point marker). The grey solid line corresponds to the star-convex split with $\gamma^* = \{1, 5\}$, to the volumetric-deviatoric split ($\gamma^* = 0$) and to the standard model ($\gamma^* = -1$)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 21 Bi-axially loaded disk with the star-convex model ($\gamma^* = 1$, $\theta = \frac{3}{4}\pi$): $R_u$ against iterations $k$ of the alternate minimization scheme for the time step at which the elastic limit is reached. All the parameters are as in Sect. 5.1 for the solid blue line ($h = \frac{\ell}{3}$). For the solid red line, only the mesh size differs ($h = \frac{\ell}{5}$)</figcaption></figure>
# 6 Conclusions

In this contribution, we focused on variational phase-field models for brittle fracture under multi-axial stress states based on energy decomposition. We first reviewed some available models of this type, namely the volumetric-deviatoric (Amor et al. 2009), the spectral (Miehe et al. 2010), the no-tension (Freddi and Royer-Carfagni 2010) and the DP-like models (De Lorenzis and Maurini 2022). We then proposed a new model that we denoted as star-convex model. All these models are of phenomenological nature and focus on reproducing macroscopic behavioral features. The major contents and findings of the paper can be summarized as follows:

* We defined essential requirements for a phase-field model of brittle fracture dealing with multi-axial stress states (Sect. 2.4): strain-hardening, stress-softening, tension/compression asymmetry, flexibility (i.e., the ability to independently calibrate not only the uniaxial tensile strength but also the uniaxial compressive strength and the shear strength), and crack-like residual stress.

* In light of these requirements we discussed the advantages and limitations of the available models. As summarized in Table 2, none of the analyzed existing decompositions was found to meet all the requirements.

* Our newly proposed star-convex model, based on a minimal modification of the volumetric-deviatoric decomposition, is equipped with a $\gamma^*$ parameter that allows for the independent calibration of compressive and tensile strengths. Such partial flexibility can be extended to the shear strength by modifying the softening laws as analyzed in (Zolesi and Maurini 2023). Additionally, the model satisfies all other requirements. Thus, it represents a very simple but effective step forward towards the realistic prediction of brittle fracture mechanisms under multiaxial stress states.

The issue of iterative convergence, which is crucial for the robustness of the simulations with all the presented models, remains open and may well be the topic of further investigations.

**Acknowledgements** This research has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 861061 – NEWFRAC Project. LDL also acknowledges funding

by the Swiss National Science Foundation (SNF) through grant N. 200021-219407 "Phase-field modeling of fracture and fatigue: from rigorous theory to fast predictive simulations". This work was performed using HPC resources from GENCI-IDRIS (Grant 2023-AD010913451R1)

**Author contributions** FV & CZ: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – Review & Editing, Visualization. PC: Conceptualization, Methodology, Writing – review & editing. CM & LDL: Conceptualization, Methodology, Writing – review & editing, Supervision, Funding acquisition.

# Declarations

**Conflict of interest** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

# References

Abrari Vajari S, Neuner M, Arunachala PK, Linder C (2023) Investigation of driving forces in a phase field approach to mixed mode fracture of concrete. Comput Methods Appl Mech Eng 417:116404

Ambati M, Gerasimov T, De Lorenzis L (2015) A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech 55(2):383–405

Ambrosio L (1992) On the approximation of free discontinuity problems. Boll Un Mat Ital B 7:105–123

Amor H, Marigo J-J, Maurini C (2009) Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. J Mech Phys Solids 57(8):1209–1229

Baldelli AAL, Maurini C (2021) Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. J Mech Phys Solids 152:104424

Bourdin B, Francfort GA, Marigo J-J (2000) Numerical experiments in revisited brittle fracture. J Mech Phys Solids 48(4):797–826

Bourdin B, Francfort GA, Marigo J-J (2008) The variational approach to fracture. J Elast 91:5–148

Chambolle A, Conti S, Francfort GA (2018) Approximation of a brittle fracture energy with a constraint of non-interpenetration. Arch Ration Mech Anal 228(3):867–889
# 6 Conclusions

In this contribution, we focused on variational phase-field models for brittle fracture under multi-axial stress states based on energy decomposition. We first reviewed some available models of this type, namely the volumetric-deviatoric (Amor et al. 2009), the spectral (Miehe et al. 2010), the no-tension (Freddi and Royer-Carfagni 2010) and the DP-like models (De Lorenzis and Maurini 2022). We then proposed a new model that we denoted as star-convex model. All these models are of phenomenological nature and focus on reproducing macroscopic behavioral features. The major contents and findings of the paper can be summarized as follows:

* We defined essential requirements for a phase-field model of brittle fracture dealing with multi-axial stress states (Sect. 2.4): strain-hardening, stress-softening, tension/compression asymmetry, flexibility (i.e., the ability to independently calibrate not only the uniaxial tensile strength but also the uniaxial compressive strength and the shear strength), and crack-like residual stress.

* In light of these requirements we discussed the advantages and limitations of the available models. As summarized in Table 2, none of the analyzed existing decompositions was found to meet all the requirements.

* Our newly proposed star-convex model, based on a minimal modification of the volumetric-deviatoric decomposition, is equipped with a $\gamma^*$ parameter that allows for the independent calibration of compressive and tensile strengths. Such partial flexibility can be extended to the shear strength by modifying the softening laws as analyzed in (Zolesi and Maurini 2023). Additionally, the model satisfies all other requirements. Thus, it represents a very simple but effective step forward towards the realistic prediction of brittle fracture mechanisms under multiaxial stress states.

The issue of iterative convergence, which is crucial for the robustness of the simulations with all the presented models, remains open and may well be the topic of further investigations.

**Acknowledgements** This research has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 861061 – NEWFRAC Project. LDL also acknowledges funding

by the Swiss National Science Foundation (SNF) through grant N. 200021-219407 "Phase-field modeling of fracture and fatigue: from rigorous theory to fast predictive simulations". This work was performed using HPC resources from GENCI-IDRIS (Grant 2023-AD010913451R1)

**Author contributions** FV & CZ: Conceptualization, Methodology, Software, Formal analysis, Investigation, Data curation, Writing – original draft, Writing – Review & Editing, Visualization. PC: Conceptualization, Methodology, Writing – review & editing. CM & LDL: Conceptualization, Methodology, Writing – review & editing, Supervision, Funding acquisition.

# Declarations

**Conflict of interest** The authors declare no competing interests.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

# References

Abrari Vajari S, Neuner M, Arunachala PK, Linder C (2023) Investigation of driving forces in a phase field approach to mixed mode fracture of concrete. Comput Methods Appl Mech Eng 417:116404

Ambati M, Gerasimov T, De Lorenzis L (2015) A review on phase-field models of brittle fracture and a new fast hybrid formulation. Comput Mech 55(2):383–405

Ambrosio L (1992) On the approximation of free discontinuity problems. Boll Un Mat Ital B 7:105–123

Amor H, Marigo J-J, Maurini C (2009) Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments. J Mech Phys Solids 57(8):1209–1229

Baldelli AAL, Maurini C (2021) Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. J Mech Phys Solids 152:104424

Bourdin B, Francfort GA, Marigo J-J (2000) Numerical experiments in revisited brittle fracture. J Mech Phys Solids 48(4):797–826

Bourdin B, Francfort GA, Marigo J-J (2008) The variational approach to fracture. J Elast 91:5–148

Chambolle A, Conti S, Francfort GA (2018) Approximation of a brittle fracture energy with a constraint of non-interpenetration. Arch Ration Mech Anal 228(3):867–889
Chao Correas A, Sapora A, Reinoso J, Corrado M, Cornetti P (2023) Coupled versus energetic nonlocal failure criteria: A case study on the crack onset from circular holes under biaxial loadings. Eur J Mech A/Solids 101:105037

De Lorenzis L, Maurini C (2022) Nucleation under multi-axial loading in variational phase-field models of brittle fracture. Int J Fract 237(1-2):61–81

Del Piero G, Owen DR (1993) Structured deformations of continua. Arch Ration Mech Anal 124:99–99

Fei F, Choo J (2021) Double-phase-field formulation for mixed-mode fracture in rocks. Comput Methods Appl Mech Eng 376:113655

Feng Y, Li J (2022) Phase-field method with additional dissipation force for mixed-mode cohesive fracture. J Mech Phys Solids 159:104693

Francfort GA, Marigo J-J (1998) Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids 46(8):1319–1342

Freddi F, Royer-Carfagni G (2010) Regularized variational theories of fracture: a unified approach. J Mech Phys Solids 58(8):1154–1174

Freddi F, Royer-Carfagni G (2011) Variational fracture mechanics to model compressive splitting of masonry-like materials. Annal Solid Struct Mech 2:57–67

Gerasimov T, De Lorenzis L (2016) A line search assisted monolithic approach for phase-field computing of brittle fracture. Comput Methods Appl Mech Eng 312:276–303

Hakimzadeh M, Agrawal V, Dayal K, Mora-Corral C (2022) Phase-field finite deformation fracture with an effective energy for regularized crack face contact. J Mech Phys Solids 167:104994

Hansen G, Herburt I, Martini H, Moszyńska M (2020) Starshaped sets. Aequ Math 94:1001–1092

He Q-C, Shao Q (2019) Closed-form coordinate-free decompositions of the two-dimensional strain and stress for modeling tension-compression dissymmetry. J Appl Mech 86(3):031007

Kirsch EG (1898) Die Theorie der Elastizität und die Bedürfnisse der Festigkeitslehre. Zeitshrift des Vereines deutscher Ingenieure 42:797–807

Kumar A, Bourdin B, Francfort GA, Lopez-Pamies O (2020) Revisiting nucleation in the phase-field approach to brittle fracture. J Mech Phys Solids 142:104027

Lancioni G, Royer-Carfagni G (2009) The variational approach to fracture mechanics. a practical application to the french panthéon in paris. J Elast 95:1–30

Lorentz E (2017) A nonlocal damage model for plain concrete consistent with cohesive fracture. Int J Fract 207(2):123–159

Marigo J-J (1989) Constitutive relations in plasticity, damage and fracture mechanics based on a work property. Nucl Eng Des 114(3):249–272

Marigo J-J, Kazymyrenko K (2019) A micromechanical inspired model for the coupled to damage elasto-plastic behavior of geomaterials under compression. Mech Indus. https://doi.org/10.1051/meca/2018043

Marigo J-J, Maurini C, Pham K (2016) An overview of the modelling of fracture by gradient damage models. Meccanica 51(12):3107–3128. https://doi.org/10.1007/s11012-016-0538-4

Miehe C, Welschinger F, Hofacker M (2010) Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations. Int J Numer Methods Eng 83(10):1273–1311

Miehe C, Schaenzel L-M, Ulmer H (2015) Phase field modeling of fracture in multi-physics problems. part i. balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Comput Methods Appl Mech Eng 294:449–485

Navidtehrani Y, Betegon C, Martinez-Paneda E (2022) A general framework for decomposing the phase field fracture driving force, particularised to a drucker-prager failure surface. Theor Appl Fract Mech 121:103555

Negri M (2010) A comparative analysis on variational models for quasi-static brittle crack propagation

Pham K, Marigo J-J (2010) Approche variationnelle de l'endommagement?: I. Les concepts fondamentaux. Compt Rend Mécanique 338(4):191–198. https://doi.org/10.1016/j.cme.2010.03.009

Pham K, Marigo J-J (2013) Stability of homogeneous states with gradient damage models: size effects and shape effects in the three-dimensional setting. J Elast 110:63–93

Pham K, Marigo J-J, Maurini C (2011a) The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models. J Mech Phys Solids 59(6):1163–1190

Pham K, Amor H, Marigo J-J, Maurini C (2011b) Gradient damage models and their use to approximate brittle fracture. Int J Damage Mech 20(4):618–652

Romani R, Bornert M, Leguillon D, Le Roy R, Sab K (2015) Detection of crack onset in double cleavage drilled specimens of plaster under compression by digital image correlation-theoretical predictions based on a coupled criterion. Eur J Mech A/Solids 51:172–182

Sacco E (1990) Modellazione e calcolo di strutture in materiale non resistente a trazione. Atti della Accademia Nazionale dei Lincei. Classe di Scienze Fisiche, Matematiche e Naturali. Rend Lincei-Mat Appl 1(3):235–258

Scroggs MW, Baratta IA, Richardson CN, Wells GN (2022a) Basix: a runtime finite element basis evaluation library. J Open Source Softw 7(73):3982

Scroggs MW, Dokken JS, Richardson CN, Wells GN (2022b) Construction of arbitrary order finite element degree-of-freedom maps on polygonal anfd polyhedral cell meshes. ACM Trans Math Softw 48(2):18:1-18:23

Steinke C, Kaliske M (2019) A phase-field crack model based on directional stress decomposition. Comput Mech 63:1019–1046

Storm J, Supriatna D, Kaliske M (2020) The concept of representative crack elements for phase-field fracture: Anisotropic elasticity and thermo-elasticity. Int J Numer Methods Eng 121(5):779–805

Strobl M, Seelig T (2016) On constitutive assumptions in phase field approaches to brittle fracture. Procedia Struct Integr 2:3705–3712