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

Phase-ﬁeld 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-ﬁeld variable. In light of the identiﬁed requirements, we review some of the existing energy decompositions, showcasing their merits and limitations, and conclude that none of them is able to fulﬁl all requirements. As a partial remedy to this outcome, we propose a new energy decomposition, denoted as star-convex model, which involves a minimal modiﬁcation of the volumetric-deviatoric


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 sharpcrack) model is limited by the difficulty to handle displacement jumps in the numerical setting.Phasefield 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 , which defines the typical width of the localization bands representing the approximation of the cracks.
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 damage model evolution, localization of damage naturally occurs upon the loss of stability of the homogeneous 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 fracture toughness of the material and are fixed, the standard 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 discussion in (Kumar et al. 2020;De Lorenzis and Maurini 2022).Flexibility in defining the strength envelope represents the main problem of nucleation under multiaxial stress states.As for crack propagation under such stress states, the most relevant problem is to model unilateral contact appropriately (Amor et al. 2009).
There is a wealth of literature proposing phasefield models for brittle fracture under multi-axial stress states.The performance of these models is demonstrated through different examples, making their comparison not immediate and their strength and limitations not obvious.Among these contributions, some preserve a variational nature; the most popular variational solutions (Amor et al. 2009;Miehe et al. 2010) are based on elastic energy decompositions.This idea is adopted in (Freddi andRoyer-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 phasefield 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 σ ) 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 0 and Poisson's ratio ν 0 , we denote by (λ 0 , μ 0 , κ 0 ) the Lamé constants and the bulk modulus given by: Given a scalar valued function: we define its positive part and negative part as: (3)

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.

General formulation
Let us consider a homogeneous body occupying the domain ⊂ R n .Its current state at point x ∈ is described by the displacement field u(x) and the irreversible scalar damage field α(x) ∈ [0, 1], with α = 0 and α = 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 where the first term ϕ 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 ε at fixed α, i.e., ϕ(r ε, α) = r 2 ϕ(ε, α) ∀r ≥ 0, and to decrease with respect to α at fixed ε.The dissipated energy density is composed of a local and a non-local part.The local term w 1 w(α) corresponds to the amount of energy dissipated per unit volume to damage homogeneously a pristine material; the dissipation function w(α) is a non-negative increasing function of α 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 is an internal length.
The total energy at time t, E t (u, α), is the sum of the strain energy and the potential energy of the external forces: where b t is the body force defined in and f t is the surface traction applied on ∂ N , both at time t, and ε(u) := sym(∇u) is the linear strain tensor.The applied displacement ūt is applied on the complementary part of the boundary ∂ D .
In the time-discrete setting of the evolution problem, given α p (the damage at the previous time t p ), the displacement and the damage field at time t = t p + t are found by solving the minimization problem: where are the spaces of the admissible displacement and damage fields at time t from the previous state with damage α p .Equation ( 6) requires (u, α) to satisfy: A necessary condition for this local constrained minimization is found taking into account only the firstorder expansion of the energy increment: where is the Gateaux derivative of the functional E t (u, α) in the direction (v, β).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 divσ (ε, α) and to the Karush-Kuhn-Tucker (KKT) conditions and boundary conditions where 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 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., α = α p .In our non-local context, their boundaries define the elastic limits for materials with homogeneous damage distribution ( α = 0).The elastic domains in the strain space R(α) and in the stress space R * (α) are defined as the sets (15) where Sym denotes the space of symmetric tensors and ϕ * (σ , α) is the complementary energy density defined as ϕ * (σ , α) At a given value of α ∈ [0, 1) a damage model enjoys the strain-hardening property if, ∀β > α, R(β) ⊃ R(α) and the stress-softening property if R * (β) ⊂ R * (α), 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).

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 where ϕ 0 is the elastic strain energy density for a homogeneous isotropic linear elastic material and a(α) is the degradation function.This function describes the degradation of the linear elastic properties with damage; it is a decreasing function of α going from a(0) = 1 to a(1) = 0.By Legendre transformation, the complementary energy density ϕ * (σ , α) is derived as . ( 18) Combining (14-18) we obtain The expressions of the elastic strain energy density in (17) and of the complementary energy density in (18) ensure that the transformation of spaces R(α) and R * (α) with varying α is a homothety centered in ε = 0 and σ = 0, respectively.The selection of the degradation function a(α) and of the dissipation function w(α) has a significant impact on the description of the damaging behavior.Two classical expressions are These two models are both strain-hardening ∀α ∈ [0, 1], and stress-softening ∀α ∈ [α peak , 1], whereas they are stress-hardening ∀α ∈ [0, α peak ).These properties hold regardless of the loading direction due to the homothetic evolution of R(α) and R * (α) with α.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 ), 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 α = α peak .For AT 1 this transition corresponds to the elastic limit because α 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 phasefield 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(α) = α).Accordingly, we define the strength surface S * as the set of stresses at the elastic limit, i.e., all stresses at the boundary of the elastic domain R * (α) for α = α peak = 0: From S * , one can define the tensile and compressive strengths, σ + e and σ − e , as the maximum and minimal allowable stress σ for the uniaxial stress state σ = σ e 1 ⊗ e 1 , and the shear strength τ e as the maximum allowable stress τ for stress states of pure shear σ = τ e 1 ⊗ e 2 .For the standard model, these quantities are given by 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 which for AT 1 gives Thus, from the experimental determination of G c and σ + e , it is possible to calibrate w 1 and .In this sense, the regularization length can be viewed as a material property.In Fig. 1, we plot the evolution of the elastic domains R(α) and R * (α) with increasing α 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 R(α) and R * (α) 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 R(α) (in blue) and half of R * (α) (in red).Note that for the elastic domains R * (α) in Fig. 1b, depicted in the plane of the principal components σ x x − σ yy , the value of σ zz is not constant but related to σ x x and σ yy through the constitutive law and the plane-strain constraint.The same is true for all the other figures depicting the elastic domains R * (α) in the principal stress components.

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 α = 0 (Figure 1): • Tension/compression symmetry: Both R(0) and R * (0) are ellipsoids in the strain and the stress 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 κ 0 , μ 0 , the tensile strength σ + e and the fracture toughness G c to be known, it is possible to calibrate w 1 and using ( 24) and ( 25).In this manner, all the model parameters are determined, thus the shape and size of the ellipsoids R(0) and R * (0) are also fixed and it is not possible to match e.g.experimentally known values of σ − e and/or τ 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 α = 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 α = 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 σ = 0 (Figure 1).From now on, we denote the stress obtained for as residual stress.Thus, for the standard model it is σ 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 tr(ε) ≥ 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., ϕ(ε, 1) = 1 2 κ 0 tr(ε) 2 − ) can be used to approximate the non-interpenetration constraint in the sense of -convergence.

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 |σ − e | > σ + 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 τ e and σ − e (or at least one of the two) independently of σ + 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 tr(ε) < 0. Accordingly, an elastic energy density such that ϕ(ε, 1) ∝ tr(ε) 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.

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) where ϕ D (ε) and ϕ R (ε) are respectively the degradable and the residual components of the elastic strain energy density.These components are non-negative and only vanish for ε = 0; their sum yields the elastic energy density of a pristine material ϕ 0 (ε), 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 Only ϕ D (ε) contributes to the energy release rate Y (ε, α); accordingly, it is the only component which affects the elastic domains where ε(σ , α) 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 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 R(α) evolves with α as a simple homothety centered in the origin.Unlike in the case of the standard model, this is not always guaranteed for R * (α).

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: This model was constructed to recover unilateral contact conditions under compression.Indeed, in (Cham-bolle et al. 2018) it is demonstrated that the volumetricdeviatoric decomposition can be used to approximate the non-interpenetration constraint in the sense ofconvergence 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 α 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 tr(ε) > 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.

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 − 2 E0 w1 w ( 0) Vol.-dev.
2 E0 w1 w (0) s (0) − 6 μ0 w1 w ( 0) ε i being the eigenvalues of the strain tensor and e i the corresponding eigenvectors.The residual stress tensor 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(α) and R * (α) are detailed in Appendix A.2.The stress-based domain R * (α) does not evolve as a simple homothety centered in the origin with respect to α, 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 * (α) is not a homotethy, we cannot define a unique damage value α peak which marks the transition between stress-hardening and stress-softening behavior.With increasing damage, R * (α) 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 * (α) 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).

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 microcracks leads to a reduction in the elastic energy density of the sound material ϕ 0 (ε) due to the presence of inelastic deformations known as structured deformations, denoted as η.These structured deformations are constrained within a convex set K ε , which defines the admissible "structure" of the micro-cracks.With ϕ 0 (ε) and K ε determined, the residual elastic energy density is computed by solving the following minimization problem η(ε) := arg min The result of the minimization problem (36), reported in Appendix A.3, is applied to the elastic energy density decomposition to derive ϕ D (ε) = ϕ 0 (ε) − ϕ R (ε).By knowing ϕ D and ϕ R , the residual stresses can be obtained as follows: 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 volumetricdeviatoric 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 → 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).

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 K ε .De Lorenzis and Maurini (De Lorenzis and Maurini 2022) leverage this flexibility by introducing a new decomposition that depends on an additional material parameter γ .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 (37) 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 With this model, the shear or the compressive strength can be calibrated independently of the tensile strength through the new parameter γ (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 notension model.

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.

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 -convergence.This success hinges on a residual elas-  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 star-convex model even though it is first introduced in the following section.The star-convex model meets all requirements with partial flexibility

Strain-hardening
Stress-softening Tens./compr.asymmetry Flexibility Crack-like residual stress a 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 b 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 tic strain energy density depending solely on the volumetric energy associated to a negative strain trace, i.e., ϕ R (ε) ∝ tr(ε) 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 σ − e independently of the tensile strength σ + e .To compensate for this disadvantage but keep a residual energy density (and stress) exclusively related to a local volume contraction tr(ε) − , we propose the following energy decomposition where γ ≥ −1 is the additional parameter controlling the σ − e /σ + e ratio.Hence, the volumetric and deviatoric components of the stress are derived as For a reason that will become clear later, we denote the corresponding model as star-convex model (or starconvex 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.

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 For tr(ε) ≥ 0, the new model shares the same boundary of the elastic domain ∂R(α) with the standard model and ∂R(α) is an ellipse in the tr(ε)−|ε dev | diagram.On the same diagram but for tr(ε) < 0, the parameter γ determines the type of conic section (Fig. 7) such that In particular, we retrieve the standard model when γ = −1 and the volumetric-deviatoric split for γ = 0.
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 elastoplastic 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 ∀α ∈ [0, 1).Additionally, for γ > 0, the hyperbola for tr(ε) < 0 evolves but keeps a constant asymptote as this depends only on the elastic constants and on γ .Since ∂R(α) 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 γ , thus providing the desired flexibility in nucleation.At the limit α = 1, the star-convex elastic domain becomes R(1) = Sym, meaning that for a fully broken material the whole strain space is admissible.
The star-convex model enjoys the stress-softening property ∀α ∈ [0, 1).At the limit α = 1, the starconvex 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 thus, the mapping between admissible strains and admissible stresses for α = 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 α = 1 is κ 0 (1+γ ).On the other hand, according to (39), it is equal to κ 0 when α = 0. Therefore, for negative volumetric strains, the volumetric stiffness of a fully broken material is larger than the one of the pristine material when γ > 0. This increase can be explained by consider- 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 α = α peak = 0, i.e. ∂R * (0), we can derive the tensile strength, the shear strength and the compressive strength as (41) From ( 41) we see that the shear strength τ e is independent of γ .Thus, the flexibility offered by the starconvex model is exclusively related to the compressive stength σ − e .

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 quasistatic 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 → E t (u, •) and α → E t (•, α) (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 , α k−1 ), we first solve the problem in u using a Newton-based nonlinear solver with line search and a maximum of 100 iterations, which gives us the output u k .Then, starting from (u k , α k−1 ), we solve the constrained nonlinear problem in α using a reduced-space active set solver for variational inequalities based on Newton's method with a maximum of 100 iterations, which returns α k .Consequently, we calculate R u , the L 2 -norm of the residual for the displacement problem for (u k , α k ).The convergence of the alternate minimization scheme is achieved if R u is below the tolerance level (tol =10 −6 ) .

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, ν 0 = 0.3, w 1 = 1.5 and = 0.04.At t = 0, the material is assumed to be intact, i.e., α = 0 everywhere.The goal of the test is to numerically determine the elastic domains of the various models for α = 0 in the ε x x -ε yy and σ x x -σ yy planes.Defining on the ε x x -ε yy plane the angle θ ∈ [0, 2π), for a given θ we prescribe a displacement ūt = t (x cos(θ ), y sin(θ )) on the disk boundary.In the linearly elastic regime, this Dirichlet boundary condition ensures that the disk is subjected to a homogeneous strain along the θ -direction, i.e., Maintaining θ fixed, we increment the value of t until we reach the elastic limit, while always imposing the Dirichlet boundary condition α = 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 = /3.
Figures 12 and 13 illustrate the contour of the elastic domain in the strain and stress planes at α = 0 for all models.The blue and red lines correspond to the theoretical 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 provides 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(ε) ≥ 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 spectral 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 transformation 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 α = 1.
Figure 12e and 12f depict the results of the test using the DP-like model with γ = √ 2μ 0 /κ 0 and γ = √ 10μ 0 /κ 0 , respectively.Similar to the spectral and no-tension models, the domains remain open, hence for some values of θ no nucleation occurs.However, the opening angle of the no-nucleation region can now be adjusted using the additional parameter γ , namely, it becomes wider with increasing γ .At the same time the region of ∂R(α) and ∂R * (α) shared with the standard model becomes narrower as γ 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 γ = √ 10μ 0 /κ 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 tr(ε) > 0.
Finally, Fig. 13a and 13b illustrate the results for the star-convex split with γ = {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 γ , 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 tr(ε) ≥ 0 and proportional to the negative volumetric strain.

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 similar nucleation test under multiaxial loading, comparing the standard model and the no-tension model to the prediction of cohesive-zone models and the finite-fracture mechanics approach.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 -θ .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) = −t L with t going from 0 up to 0.3.At t = 0, the material is assumed to be intact, i. e., α = 0 everywhere.We adopt L = 1, R = 0.3, E 0 = 100, ν 0 = 0.3, w 1 = 1 and = 0.02.We discretize the problem with standard bilinear quadrilateral finite elements using an unstructured mesh with average size h = /3.
During the elastic regime, the highest stress concentration is expected at the boundary of the hole (r = R), along which the radial stress σ rr and the tangential stress σ r θ are both zero (σ rr = σ r θ = 0) and the hoop stress σ θθ and the out-of-plane stress σ zz are related by the plane-strain condition σ zz = −ν 0 σ θθ .Along the boundary of the hole, point A (r = R, θ = π/2) undergoes the maximum tensile hoop stress σ A θθ ≥ 0 whereas point B (r = R, θ = 0) is subjected to the lowest compressive hoop stress σ B θθ ≤ −σ 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 ∂R * (0).For first nucleation in A, the stress state in A must reach the intersection between ∂R * (0) and the half-line σ zz = −ν 0 σ θθ ≤ 0. We denote with σ + θθ 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 ∂R * (0) and the half-line σ zz = −ν 0 σ θθ ≥ 0. We denote with σ − θθ e the value of the hoop stress at this intersection point (Fig. 15).The  assumption, respectively.Their intersection determines the compressive limit σ − θθe (blue on the left) and the tensile limit σ + θθe (blue on the right).Compressive-tensile hoop limits ratio (b): trend of the limits ratio vs. γ (red solid line).For γ > γ 3 the ratio is greater than 3, for γ ≥ γ ∞ the ratio is equal to ∞ can infer that if the |σ − θθ e /σ + θθ 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 |σ − θθ e /σ + θθ e | ratio is well below 3 (Table 3).Unlike the other models, the DP-like and the star-convex models contain the additional parameters γ and γ that allow for the flexible calibration of the strength ratio |σ − θθ e /σ + θθ e |.For the star-convex model, Fig. 15 shows how the |σ − θθ e /σ + θθ e | ratio is affected by the variation of γ .
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 σ 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- 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 volumetricdeviatoric, spectral, and star-convex split with γ = 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 γ (DP-like model) or γ (star-convex split), along with the no-tension model.However, in the case of the DP-like model (Fig. 17d), for γ = √ 10μ 0 /κ 0 , the damage exhibits a diffuse pattern rather than being fully localized.This is attributed to the phenomenon of crack thickening as γ increases, as previously highlighted in the bi-axially loaded disk test.

Sliding test
The last test concerns a square plate of side L with a crack in the center, as shown in Fig. 18, under planestrain 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 α(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) = t L is applied to the upper part.Furthermore, the two sides are vertically blocked in their lower part u y (x = 0, y ≤ L/2) = u y (x = L , y ≤ L/2) = 0 and we apply a constant (i.e., independent of t) vertical displacement on the upper part u y (x = 0, y ≥ L/2) = u y (x = L , y ≥ L/2) = ūL, where ū = 0.1.This test aims to emphasize the Fig. 18 Geometry and loading for the sliding test 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, ν 0 = 0.3, w 1 = 1 and = 0.05.We discretize the problem with standard bilinear quadrilateral finite elements using a structured mesh with h = /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 γ = √ 10μ 0 /κ 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 σ 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 tr(ε) − .The spectral and the no-tension models, as well as the DP-like model for γ = √ 10μ 0 /κ 0 , exhibit a significant reaction force, slightly lower than half of that observed in an undamaged material.The DP-like model for γ = √ 2μ 0 /κ 0 initially displays no reaction and a deformation as shown in Fig. 19, which is acceptable.However, it subsequently shows an increase in the reaction force.This is because a lower γ value delays the undesired residual stress effect.

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 θ 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.4 π : 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 = 3 .For the solid red line, only the mesh size differs h = 5

Conclusions
In this contribution, we focused on variational phasefield 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, stresssoftening, 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 γ 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.

Fig. 3
Fig. 3 Elastic domains for the model with spectral split in the strain (blue dash-dotted line) and stress (red solid line) space for α = {0, 0.25, 0.5, 0.75, 1} (ν 0 = 0.3).The domains for α = 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

Fig. 4
Fig. 4 Elastic domains for the no-tension model in the strain (blue dash-dotted line) and stress (red solid line) space for α = {0, 0.25, 0.5, 0.75, 1}(ν 0 = 0.3).The domains for α = 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 planestrain assumption

Fig. 11
Fig. 11 Geometry and loading for the bi-axially loaded disk

Fig. 12
Fig. 12 Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the available models

Fig. 13
Fig. 13 Analytical (lines) versus numerical (dots) nucleation domains and damage fields after nucleation for the star-convex model |σ − θθ e /σ + θθ e | ratio increases with the compressive to tensile strength ratio |σ − e /σ + e | with a trend which is specific to each model.Since when L R, σ B θθ = −3σ A θθ (Kirsch 1898), we

Fig. 15
Fig. 15Effect of γ 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 ∂R * (0) for γ = {−1, 0, 1, 2, 5} and the stress state under the plane-strain

Fig. 16 Fig. 17
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

Fig. 19
Fig. 19 Damage fields on the deformed configurations for all models

Table 1
Tensile, compressive and shear strengths for all the analyzed models