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

HAL
open science

# On the strain energy decomposition in phase field brittle fracture: established models and novel cleavage plane-based techniques

Nour Habib, Saber El Arem, Amine Ammar

► To cite this version:

Nour Habib, Saber El Arem, Amine Ammar. On the strain energy decomposition in phase field brittle fracture: established models and novel cleavage plane-based techniques. Advanced Modeling and Simulation in Engineering Sciences, 2025, 12, 10.1186/s40323-025-00295-5. http://hal-05093095

HAL Id: hal-05093095

https://hal.science/hal-05093095v1

Submitted on 2 Jun 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# On the strain energy decomposition in phase field brittle fracture: established models and novel cleavage plane-based techniques

Nour Habib ID, Saber El Arem and Amine Ammar

LAMPA Arts et Metiers Institute
of Technology, 2 boulevard du
Ronceray, BP 93525, 49035
Angers cedex 01, France

## Abstract

This work offers a detailed examination of the phase field approach for modeling brittle fracture, emphasizing its theoretical foundations, mathematical descriptions, and computational strategies. Central to our discussion is an in-depth analysis of strain energy decomposition methods integral to phase field models. We introduce an innovative technique using a cleavage plane based degradation that has shown promising results under various loading scenarios. We meticulously evaluate each method's inherent limitations and challenges to highlight their respective advantages and drawbacks across different loading scenarios. This review aims not only to catalog existing knowledge but also to pave the way for future research directions in the application of phase field approach to fracture analysis.

**Keywords:** Phase field method, Brittle fracture, Fracture mechanics, Damage, Strain energy formulation

## Introduction

Continuum damage mechanics (CDM) has emerged as a foundational framework for predicting the behavior of materials that undergo progressive degradation. It is a powerful tool for understanding and preventing failure in engineering materials as fracture and crack propagation are fundamental aspects of material behavior, and understanding and predicting crack evolution has led to significant advancements in structural design, safety, and reliability across various industries. In CDM the phenomenology of fracture and crack growth is considered through the use of a continuous damage variable that will represent the density of the cracks and voids, and the loss of stiffness or strength in a material due to the accumulation of strain. In the context of Finite Element Method (FEM), solving continuum damage mechanics holds several drawbacks that can compromise the accuracy and robustness of the analyses. These limitations encompass mesh dependency as the evolution of damage and the resultant stress distribution are highly sensitive to the discretization of the problem. Furthermore, issues such as loss of ellipticity in static scenarios and challenges in ensuring hyperbolicity in dynamic scenarios under-

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

<footnote>*Correspondence:
nour.habib@ensam.eu</footnote>
score the complexities in achieving unique and stable numerical solutions. Addressing these challenges, the Extended Finite Element Method (XFEM) [21] offers a robust alternative. By employing level sets and enrichment functions, XFEM circumvents traditional mesh constraints, enabling dynamic simulation of crack propagation without necessitating mesh modifications. Despite its advantages, XFEM's integration into existing codes bears complexity, requiring additional functions (enrichment functions) and managing extra degrees of freedom. Moreover, the enriched solution space escalates computational demands, posing challenges in large-scale or multifaceted crack scenarios.

Parallelly, non-local damage models, advocated by pioneering studies [4], incorporate the influence of neighboring zone damage, addressing mesh sensitivity through interactions defined over a physical length scale. This approach not only captures complex interactions within materials but also excels in scenarios where cracks are distributed rather than localized, as observed in quasi-brittle materials like concrete. However, the broader consideration of damage interactions necessitates significant computational resources and intricate calibration of model parameters, often relying on extensive experimental data.

The phase-field model presents an innovative mathematical framework for interface evolution between material phases. It offers a distinct advantage by blurring out boundaries and introducing a continuous development of field variables across the interfacial zone. This smearing of the interface enables a more realistic representation of crack propagation and other interfacial phenomena. The inspiration for the phase field method can be traced back to the Ginzburg-Landau equation developed to determine surface tension at the boundary between distinct phases. They introduced the concept of an order parameter and its gradient to track the values of superconductivity. The key feature of the phase field is the use of the gradient of the field variable that is added to the free energy functional to penalize the formation of sharp interfaces, ensuring that the phase-field variable remains smooth over a defined characteristic width. These advantages have led to their widespread application in simulating solidification processes, grain growth, fracture mechanics, and a variety of other phenomena. In the context of brittle fracture, the phase field approach provides a versatile framework for modeling crack propagation. Rather than explicitly tracking crack boundaries, the phase field model incorporates an order parameter. In the context of fracture, this parameter is the damage variable (d), which describes the presence and extent of cracks within the material. The idea of phase field modeling would be representing the damage variable as continuous everywhere, with a process zone of a characteristic length where it would rapidly change value to transition from the undamaged area to the cracked area. By introducing the gradient of the damage (∇d) and a crack width parameter (or characteristic length, $l_c$), gradient crack evolution can be accurately captured within the numerical model. Another alternative method for modeling diffuse crack behavior is the gradient damage model proposed by [8], drawing inspiration from a variational approach rooted in Griffith's seminal work [10]. This approach addresses the traditional methods' limitations by facilitating crack nucleation and determining crack propagation paths. Their numerical model for brittle fracture involves minimizing a weak functional, similar to Mumford-Shah's model for image segmentation and the volume integral approximation done by Ambrosio and Tortorelli, and incorporates $\Gamma$-convergence regularization for numerical implementation [5]. The Francfort-Marigo model serves as the foundation for many variational approaches
and has paved the way for phase field models of fracture. [18] also outlined a phase field model for crack propagation based on these works but from a continuum mechanics and thermodynamics outline with a gradient-type regularized surface energy. All of these non-local gradient damage models share a similar numerical structure, using regularization parameters that represent the material's internal length scale and employing gradient damage models. The distinguishing factor lies in the interpretation of the formulations and their specific decomposition approaches.

The primary focus of the various decomposition approaches lies in identifying the key factors driving damage evolution. These approaches typically involve separating the influence of active and inactive components on damage progression. The active component facilitates damage evolution, while the inactive component inhibits further damage growth. This decomposition is largely based on distinguishing tensile and compressive loading modes. Tensile loads, associated with positive strains and stresses, contribute significantly to damage progression, whereas compressive loads, linked to negative strains and stresses, are generally not considered to exacerbate damage. There exists in the literature various split approaches for the elastic strain energy to account for the unilateral aspect of damage. The most prevalent works being the preliminary decomposition by [2, 20]. Both of these rely on the split of the strain tensor. Several alternative models have been proposed in the literature, exemplifying diverse approaches to structural analysis. Notably, the model by [22] utilizes a split based on the orthogonality between traction and compression components. Additionally, the work of [26] introduces a modified perspective on the volumetric-deviatoric split, while the study presented in [9] features a reformulation of the compliance tensor. Other significant contributions include models that depend on the decomposition of the stress tensor, as exemplified in [28], and the directional split of the stress tensor, as detailed in [24]. The latter model forms part of a broader category of directional split methods, which also encompasses the studies by [7, 25]. Furthermore, we introduce a new decomposition in the subsequent sections of this work, aligning with the directional split methods.

In this work, we present a comprehensive study of phase field modeling for brittle fracture, focusing on different approaches adopted in existing models.

* We explore five established methods for the decomposition of the elastic strain energy,

* We introduce a sixth technique, a novel approach intended for anisotropic materials that involves considering a weakest cleavage plane inherent to the material,

* We compare the results of the different decomposition methods in the case of tension/compression and pure shear testing,

* We elucidate the advantages and limitations of each approach for the different loading cases.

By revisiting the remarkable strides made in the last decades, this work aims to both honor the legacy of earlier research and contribute to the ongoing development of phase field modeling, fostering a deeper comprehension and refining the methodologies employed in the study of material degradation and fracture mechanics.
<figure><img src="image_5.png" /><figcaption>Fig. 1 Diffuse fracture modeling</figcaption></figure>

## Phase field model for brittle fracture

In the phase field model for brittle fracture, the damage variable, denoted as $d$, quantifies the extent of material damage. This variable ranges from 0, representing a sound state, to 1, indicating complete material failure. The initiation and progression of cracks are governed by a critical threshold related to the surface energy, as conceptualized by Griffith. The fracture toughness, $G_c$, signifies the energy dissipation required for crack extension per unit area. Consequently, the total energy of a body, comprising volume $\Omega$ and cracked surface $\Gamma$, integrates both elastic and surface energies as follows:

$$ \mathcal{E}(\boldsymbol{\varepsilon}(\mathbf{u}), d) = \int_{\Omega} \psi(\boldsymbol{\varepsilon}, d) \, d\Omega + \int_{\Gamma} G_c d\Gamma \quad (1) $$

The elastic energy is stored within the body, while fracture energy is accounted for by integrating over the cracked surface. To incorporate the regularization of the crack surface via the phase field parameter, represented by damage in brittle fracture contexts, and to address the interface smearing across a characteristic width, $l_c$, the energy formulation is adapted using [18]'s crack surface density function $\gamma$:

$$ \mathcal{E}(\boldsymbol{\varepsilon}(\mathbf{u}), d, \nabla d) = \int_{\Omega} \psi(\boldsymbol{\varepsilon}, d) d\Omega + \int_{\Omega} G_c \gamma(d, \nabla d) d\Omega \quad (2) $$

This reformulation transforms the surface integral into a volume integral, simplifying the computational approach by eliminating the need to explicitly track the evolving crack surface. Instead, the integration spans the entire volume, treating the crack representation as an inherent part of the material's continuum description.

## The crack surface density function

The phase field model characterizes the crack surface area, $A_\Gamma$, through the crack surface density function per unit volume, $\gamma$, effectively translating the traditional concept of a discrete crack into a continuous field within the material volume $\Omega$:

$$ A_\Gamma = \int_{\Gamma} d\Gamma = \int_{\Omega} \gamma(d, \nabla d) \, d\Omega \quad (3) $$

This approach delineates the phase field model's non-local nature, with $\gamma$ as a function of both the damage variable $d$ and its gradient $\nabla d$. The characteristic length scale, $l_c$, signifies the transition zone's width from undamaged to damaged material, facilitating a smooth damage gradient rather than a sharp discontinuity (Fig. 1). As $l_c$ narrows, the model approaches the traditional sharp crack representation. A prevalent formulation for $\gamma$ is:

$$ \gamma(d, \nabla d) = \frac{1}{c_0} \left( \frac{1}{l_c} \alpha(d) + l_c \nabla d \cdot \nabla d \right) \quad (4) $$
with $c_0$ normalizing the regularization to ensure fidelity to the physical crack surface, defined as:

$$c_0 = 4 \int_{0}^{1} \sqrt{\alpha(\beta)} d\beta \qquad (5)$$

The function $\alpha(d)$ correlates the damage variable to crack geometry. Its common expression is:

$$\alpha(d) = d\xi_1 + d^2\xi_2 + d^3\xi_3 + \dots + d^m\xi_m \qquad (6)$$

with the most used form being:

$$\alpha(d) = \xi d + (1 - \xi) d^2 \qquad (7)$$

Here, $\xi$ modulates the function's shape, aligning $\alpha(d)$ and $d$ within the [0, 1] interval, ensuring $\alpha(0) = 0$ and $\alpha(1) = 1$. The literature primarily recognizes quadratic ($\xi = 0$, yielding $\alpha(d) = d^2$) and linear [$\xi = 1$, yielding $\alpha(d) = d$] forms.

The quadratic expression [18]:

$$\gamma(d, \nabla d) = \frac{1}{2l_c} (d^2 + l_c^2 \nabla d \cdot \nabla d) \qquad (8)$$

lacks an initial pure elastic response, with damage onset coinciding with loading application. In contrast, the linear form [23]:

$$\gamma(d, \nabla d) = \frac{3}{8l_c} (d + l_c^2 \nabla d \cdot \nabla d) \qquad (9)$$

introduces a threshold, delaying damage initiation until a critical elastic strain energy density is reached.

## Damage evolution equation

The progression of damage within the material, adheres to the Allen-Cahn equation, first introduced in [1], also recognized as the Ginzburg-Landau time-dependent evolution equation:

$$\dot{d} = -L \frac{\delta F}{\delta d} \qquad (10)$$

where $L$ represents the mobility factor, inversely proportional to the viscosity parameter, and $F$ denotes the Helmholtz free energy of the system that encompasses both the volumetric energy and the crack surface energy contributions:

$$F = \int_{\Omega} \psi(\boldsymbol{\epsilon}, d) d\Omega + \int_{\Omega} G_c \gamma(d, \nabla d) d\Omega \qquad (11)$$

To facilitate analysis, the Helmholtz free energy density, $\Psi$, is defined as the sum of the elastic energy density and the energy density associated with crack formation:

$$\Psi = \psi(\boldsymbol{\epsilon}, d) + G_c \gamma(d, \nabla d) \qquad (12)$$

## The degradation function

The phase field model accounts for the influence of damage on the stored energy by dividing it into two distinct components:

$$\psi(\boldsymbol{\epsilon}, d) = g(d) \psi_{0}^{+}(\boldsymbol{\epsilon}) + \psi_{0}^{-}(\boldsymbol{\epsilon}) \qquad (13)$$

Here, $\psi_0^+$ is the part of the elastic energy density contributing to damage growth, while $\psi_0^-$ denotes the energy portion that remains unaffected by damage progression. The
degradation function, g(d), modulates the impact of damage, inversely correlating with the damage variable d. As d evolves from 0 to 1, g(d) shifts from 1 to 0, hence the decrease in the $\psi_0^+$ component's contribution correlatively with escalating damage.

To accurately model material behavior under damage, g must satisfy:

* g(0) = 1, indicating that in absence of damage (d = 0), the material exhibits full stiffness with no degradation.

* g(1) = 0, marking the fully damaged state (d = 1) as a point of stiffness loss, eliminating the damaging energy contribution.

* g'(d) ≤ 0, ensuring that the degradation function monotonically decreases with increasing damage, reflecting the intuitive decrease in material stiffness.

* g'(1) = 0, ensuring the elastic energy contributing to further damage ceases at full degradation.

Commonly, the degradation function is represented as:

$$g(d) = (1 - d)^2 \qquad (14)$$

For this quadratic model, the elastic energy decomposition adapts as follows:

$$\psi(\boldsymbol{\epsilon}, d) = ((1 - d)^2(1 - k) + k)\psi_0^+(\boldsymbol{\epsilon}) + \psi_0^-(\boldsymbol{\epsilon}) \qquad (15)$$

where k represents a residual stiffness parameter, chosen to be as small as possible often taken between $10^{-6}$ and 0 to enhance numerical stability during the resolution process.

## The damage driving force

As in [18], it is possible to define the damage criterion based on the crack driving force, in a manner akin to Griffith's criterion, to determine the necessary energy for crack propagation:

$$\begin{cases} \bar{D} < 0 \Rightarrow \dot{d} = 0 \\ \bar{D} \ge 0 \Rightarrow \dot{d} \ge 0 \end{cases} \qquad (16)$$

$\bar{D}$ is the thermodynamic force associated to the damage variable d and given by:

$$\bar{D} = - \left[ \frac{\partial \psi(\boldsymbol{\epsilon}, d)}{\partial d} + G_c \delta_d \gamma \right] \qquad (17)$$

The evolution of damage, governed by the Kuhn Tucker conditions, is delineated as follows:

$$\left\{ \begin{array}{l} \dot{d} \ge 0 \\ -\frac{\partial \psi_e(\boldsymbol{\epsilon}, d)}{\partial d} \le G_c \delta_d \gamma \\ \dot{d} \left( -\frac{\partial \psi_e(\boldsymbol{\epsilon}, d)}{\partial d} - G_c \delta_d \gamma \right) = 0 \end{array} \right. \qquad (18)$$

The variation in the crack surface density function, $\delta_d \gamma$, is quantified as:

$$\delta_d \gamma = \frac{1}{c_0 l_c} (\alpha'(d) - 2l_c^2 \Delta d) \qquad (19)$$

This formulation allows for a generalized expression of $\bar{D}$:

$$\bar{D} = -g'(d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{G_c}{c_0 l_c} [\alpha'(d) - 2l_c^2 \Delta d] \qquad (20)$$

For a linear crack geometric function ($\alpha(d) = d$) and a degradation function of g(d) = (1 - d)$^2$, the crack driving force simplifies to:

$$\bar{D} = 2(1 - d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{3G_c}{8l_c} [1 - 2l_c^2 \Delta d] \qquad (21)$$
and the damage criteria becomes:

$$ \begin{cases} 2(1-d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{3G_c}{8l_c^2}[1-2l_c^2\Delta d] < 0 \Rightarrow \dot{d} = 0 \\ 2(1-d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{3G_c}{8l_c^2}[1-2l_c^2\Delta d] \ge 0 \Rightarrow \dot{d} \ge 0 \end{cases} \quad (22) $$

In the scenario with a quadratic crack geometric function, $\alpha(d) = d^2$, the damage driving force becomes:

$$ \tilde{D} = 2(1-d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{G_c}{2l_c}[2d - 2l_c^2\Delta d] \quad (23) $$

and, the damage criteria is described by:

$$ \begin{cases} 2(1-d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{G_c}{l_c}[d - l_c^2 \Delta d] < 0 \Rightarrow \dot{d} = 0 \\ 2(1-d)\psi_0^+(\boldsymbol{\epsilon}) - \frac{G_c}{l_c}[d - l_c^2 \Delta d] \ge 0 \Rightarrow \dot{d} \ge 0 \end{cases} \quad (24) $$

### Irreversibility criterion: discussion on the use of history variable

To ensure the irreversibility of damage growth and inhibit crack healing, [18] introduced a history-dependent mechanism. With this mechanism, a history variable, $\mathcal{H}$, captures the maximum value attained by the local contribution to the driving force, $\psi_0^+$, across the loading history:

$$ \mathcal{H} = \max_t (\psi_0^+(\boldsymbol{\epsilon})) \quad (25) $$

This history variable fulfills the Kuhn–Tucker conditions, ensuring the irreversibility of damage:

$$ \begin{cases} \mathcal{H} - \psi_0^+(\boldsymbol{\epsilon}) \ge 0 \\ \dot{\mathcal{H}} \ge 0 \\ \dot{\mathcal{H}}(\psi_0^+(\boldsymbol{\epsilon}) - \mathcal{H}) = 0 \end{cases} \quad (26) $$

Though this history function can only be used freely in the case of a quadratic crack geometric function because of its local minima at $d = 0$. At the onset of damage ($d = 0$), the critical strain energy, $\mathcal{H}_c$, necessary for crack propagation initiation, is determined by:

$$ \frac{G_c}{c_0 l_c} \alpha'(d) + \mathcal{H}_c g'(0) = 0 \Rightarrow \mathcal{H}_c = -\frac{1}{g'(0)} \frac{G_c}{c_0 l_c} \alpha'(d) \quad (27) $$

$\mathcal{H}_c$ signifies the threshold of strain energy for damage initiation. Adopting a quadratic geometric function yields

$$ \mathcal{H}_c = 0 \quad (28) $$

suggesting that damage evolves immediately upon loading. But for a linear crack geometric function ($\alpha(d) = d$, $c_0 = \frac{8}{3}$) and considering $g(d) = (1-d)^2$, we obtain:

$$ \mathcal{H}_c = \frac{3G_c}{16l_c} \quad (29) $$

When replacing $\mathcal{H}$ in the damage evolution equation for resolution of the numerical problem, in the initial stage $\mathcal{H} < \mathcal{H}_c$ will lead to negative values for the damage. Consequently, to avoid the negative damage values the history variable becomes:

$$ \mathcal{H} = \max_t (\psi_0^+(\boldsymbol{\epsilon}), \mathcal{H}_c) \quad (30) $$

The use of the history function defined in the previous section to ensure the irreversibility leads to an excessive crack spread perpendicularly to the crack propagation direction.
<figure><img src="image_6.png" /><figcaption>Fig. 2 Diffusive crack topology for linear crack geometric function with horizontal sharp crack and a numerical implementation with a history variable condition to ensure irreversibility, $l_c = 0.005$ mm</figcaption></figure>

In [27], for a linear crack geometric function ($\alpha(d) = d$), the half-bandwidth of the crack $D$ is:

$$D = 2l_c \qquad (31)$$

And for the quadratic crack geometric function ($\alpha(d) = d^2$):

$$D = +\infty \qquad (32)$$

In Fig. 2 we have a 1 × 1 mm structure with pre-existing crack defined by enforcing Dirichlet boundary condition setting $d = 1$ over the cracked zone, and with $l_c = 0.005$ mm. The structure remains unexposed to any mechanical loading. When the irreversibility is assured by the history function, it is noticeable that the damage is diffused over a zone larger than expected with $D > 10l_c = 0.05$ mm as depicted in Fig. 2.

Whereas, when using the bound constrained optimization solver on the damage variable as done in the Fig. 3, the half-bandwidth is equal to $2l_c = 0.01$ mm.

Which is why in our numerical implementation, the irreversibility of damage is ensured through the bound constrained optimization scheme on the damage variable, rather than relying on the history variable. This approach guarantees the conservation and proper sequencing of the free energy term at each time step within the coupled system.

For subsequent discussions and numerical implementations presented in this paper, we will be considering:

* a linear crack geometric function $\alpha(d) = d$,

* alongside a quadratic degradation function $g(d) = (1 - d)^2$,

* the use of bound constrained optimization scheme to ensure damage irreversibility $\dot{d} \ge 0$ and $d \in [0, 1]$, using variational inequality solver SNES provided by PETSc [3].

## The elastic energy decomposition

This section explores the methodologies employed to partition the elastic strain energy, distinguishing between components that contribute to damage evolution and those that do not influence the damage process.
<figure><img src="image_7.png" /><figcaption>Fig. 3 Diffusive crack topology for linear crack geometric function with horizontal sharp crack and a numerical implementation with solver constraints, $l_c = 0.005$ mm</figcaption></figure>

### Symmetric formulation (no energy decomposition)

In the symmetric (or isotropic) formulation, the entire energy stored in the solid bulk contributes to damage.

$$ \psi_0^+(\boldsymbol{\epsilon}) = \psi_0(\boldsymbol{\epsilon}), \quad \psi_0^-(\boldsymbol{\epsilon}) = 0 \qquad (33) $$

Thus giving:

$$ \psi(\boldsymbol{\epsilon}, d) = g(d)\psi_0(\boldsymbol{\epsilon}) \qquad (34) $$

where $\psi$ denotes the stored energy density, $g(d)$ represents the degradation function, and $\psi_0$ is the elastic strain energy density of the undamaged elastic solid

$$ \psi_0(\boldsymbol{\epsilon}) = \frac{1}{2}\kappa \mathrm{Tr}^2(\boldsymbol{\epsilon}) + \mu \mathrm{Tr}(\boldsymbol{\epsilon}^2) \qquad (35) $$

Here, $\kappa = \lambda + 2\mu/n$ signifies the bulk modulus, with $\lambda$ and $\mu$ being the Lamé coefficients, and $n$ the dimension number.

The stress tensor is thus given by:

$$ \boldsymbol{\sigma} = \frac{\partial \psi(\boldsymbol{\epsilon}, d)}{\partial \boldsymbol{\epsilon}} = g(d) \frac{\partial \psi_0(\boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}} = g(d) \boldsymbol{\sigma}_0 = g(d) (\lambda \mathrm{Trace}(\boldsymbol{\epsilon})\mathbb{I} + 2\mu \boldsymbol{\epsilon}) \qquad (36) $$

This symmetric formulation neglects the possibility of crack closure effects, which can occur under compressive loads, potentially mitigating damage progression and altering the material's response.

### The deviatoric-spherical decomposition of the strain tensor

Introduced by [2], this method builds upon the strain tensor's division into spherical and deviatoric components [14]. We write:

$$ \boldsymbol{\epsilon} = \boldsymbol{\epsilon}_d + \boldsymbol{\epsilon}_s \qquad (37) $$

with $\boldsymbol{\epsilon}_s = \frac{1}{n} \mathrm{Tr}(\boldsymbol{\epsilon})\mathbb{I}$, and $\boldsymbol{\epsilon}_d = \boldsymbol{\epsilon} - \boldsymbol{\epsilon}_s$ and $\mathbb{I}$ is the identity tensor.

According to [2], all deviatoric strains and the positive spherical strain are deemed damaging. The elastic strain energy density is then formulated as:

$$ \psi_0(\boldsymbol{\epsilon}) = \frac{1}{2}\kappa \langle \mathrm{Tr}(\boldsymbol{\epsilon}) \rangle_+^2 + \frac{1}{2}\kappa \langle \mathrm{Tr}(\boldsymbol{\epsilon}) \rangle_-^2 + \mu \boldsymbol{\epsilon}_d : \boldsymbol{\epsilon}_d \qquad (38) $$
The Macaulay brackets $\langle \cdot \rangle$ discern the positive and negative traces:

$$ \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_{+} = \max(\mathrm{Tr}(\boldsymbol{\varepsilon}), 0), \quad \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_{-} = \max(-\mathrm{Tr}(\boldsymbol{\varepsilon}), 0) \qquad (39) $$

$\psi_0$ splits into $\psi_0^+$ and $\psi_0^-$ components:

$$ \psi_0^+(\boldsymbol{\varepsilon}) = \frac{1}{2}\kappa \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_+^2 + \mu \boldsymbol{\varepsilon}_d : \boldsymbol{\varepsilon}_d, \quad \psi_0^-(\boldsymbol{\varepsilon}) = \frac{1}{2}\kappa \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_-^2 \qquad (40) $$

and the stress tensor is given by:

$$ \boldsymbol{\sigma} = g(d) \frac{\partial \psi_0^+(\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}} + \frac{\partial \psi_0^-(\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}} = g(d) \boldsymbol{\sigma}^+ + \boldsymbol{\sigma}^- \qquad (41) $$

with

$$ \boldsymbol{\sigma}^+ = \kappa \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_{+} \mathbb{I} + 2\mu \boldsymbol{\varepsilon}_d, \quad \boldsymbol{\sigma}^- = \kappa \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_{-} \mathbb{I} \qquad (42) $$

This decomposition method presumes that volumetric compression does not facilitate damage formation and propagation. Conversely, volumetric expansion indicates potential damage development, leading to the release of elastic strain energy. Moreover, damage evolution may occur irrespective of volumetric changes, as induced by shear strains contributing to damage progression.

## The spectral decomposition of the strain tensor

The spectral decomposition method, introduced by [18], effectively separates the elastic strain energy, $\psi_0$, into two components based on the strain's nature: the 'positive' component, $\psi_0^+$, which is driven by tensile strains conducive to damage evolution, and the 'negative' component, $\psi_0^-$, which results from compressive strains and is considered non-damaging, hence not influencing damage development. This differentiation is achieved through spectral decomposition of the strain tensor, $\boldsymbol{\varepsilon}$, into eigenvalues, $\varepsilon_i$, and their corresponding eigenvectors, $\mathbf{n}_i$:

$$ \boldsymbol{\varepsilon} = \sum_i \varepsilon_i \mathbf{n}_i \otimes \mathbf{n}_i \qquad (43) $$

The positive and negative components of the strain energy density are expressed as:

$$ \psi_0^+(\boldsymbol{\varepsilon}) = \frac{\lambda \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_+^2}{2} + \mu (\langle \varepsilon_1 \rangle_+^2 + \langle \varepsilon_2 \rangle_+^2 + \langle \varepsilon_3 \rangle_+^2) \qquad (44) $$

$$ \psi_0^-(\boldsymbol{\varepsilon}) = \frac{\lambda \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_-^2}{2} + \mu (\langle \varepsilon_1 \rangle_-^2 + \langle \varepsilon_2 \rangle_-^2 + \langle \varepsilon_3 \rangle_-^2) \qquad (45) $$

and the stress is computed using:

$$ \boldsymbol{\sigma} = g(d) \frac{\partial \psi_0^+(\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}} + \frac{\partial \psi_0^-(\boldsymbol{\varepsilon})}{\partial \boldsymbol{\varepsilon}} = g(d) \boldsymbol{\sigma}^+ + \boldsymbol{\sigma}^- \qquad (46) $$

with

$$ \boldsymbol{\sigma}^+ = \sum_i \left( \lambda \left\langle \sum_i \varepsilon_i \right\rangle_+ + 2\mu \langle \varepsilon_i \rangle_+ \right) \mathbf{n}_i \otimes \mathbf{n}_i \\ = \lambda \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_+ \mathbb{I} + 2\mu (\langle \varepsilon_1 \rangle_+ + \langle \varepsilon_2 \rangle_+ + \langle \varepsilon_3 \rangle_+) \qquad (47) $$

and

$$ \boldsymbol{\sigma}^- = \sum_i \left( \lambda \left\langle \sum_i \varepsilon_i \right\rangle_- + 2\mu \langle \varepsilon_i \rangle_- \right) \mathbf{n}_i \otimes \mathbf{n}_i \\ = \lambda \langle \mathrm{Tr}(\boldsymbol{\varepsilon}) \rangle_- \mathbb{I} + 2\mu (\langle \varepsilon_1 \rangle_- + \langle \varepsilon_2 \rangle_- + \langle \varepsilon_3 \rangle_-) \qquad (48) $$
The model accounts for tensile strains entirely in the damage process, attributes a partial contribution to shear strains, and disregards compressive strains. Notably, the initial term on the right side of the equations (involving the $\lambda$ coefficient) considers the sign of the strain tensor's trace—indicating volumetric change—without segmenting into individual strain components, unlike the latter term.

### The spectral decomposition of the stress tensor

The spectral decomposition of the stress tensor, introduced by [28] and building upon the foundational work of [18], applies the principles of strain tensor decomposition to the unbroken (undamaged) stress tensor, $\sigma_0$. The spectral decomposition is done along the eigenvalues $\sigma_i$ and the eigenvectors $s_i$:

$$ \sigma_0 = \sum_i \sigma_i s_i \otimes s_i \qquad (49) $$

In alignment with the strain decomposition approach, tensile stresses—indicated by positive eigenvalues—contribute to damage evolution, while compressive stresses, denoted by negative eigenvalues, are considered non-contributory to damage progression. The construction of positive ($P^+$) and negative ($P^-$) projection tensors follows the methodologies outlined by [19].

$$ P^+ = \partial_{\sigma}[\sigma_0^+], \quad P^- = I - P^+ \qquad (50) $$

The resultant stress tensors, representing the damaging ('positive') and non-damaging ('negative') components of stress, are expressed as:

$$ \sigma_0^+ = P^+ : \sigma_0 \qquad (51) $$

$$ \sigma_0^- = P^- : \sigma_0 \qquad (52) $$

Subsequently, the elastic strain energy positive and negative contributions are given by:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{1}{2} \sigma_0^+ : \boldsymbol{\epsilon}, \quad \psi_0^-(\boldsymbol{\epsilon}) = \frac{1}{2} \sigma_0^- : \boldsymbol{\epsilon} \qquad (53) $$

And the stress tensor is given as:

$$ \sigma = g(d)\sigma_0^+ + \sigma_0^- \qquad (54) $$

### The orthogonal decomposition of the strain tensor

Building upon the innovations of [12,22], this method introduces an orthogonal decomposition approach to the stiffness tensor. This technique is predicated on dividing the fourth-order stiffness tensor, $C_0$, which is fundamental to the formulation of elastic strain energy density.

The stiffness tensor, $C_0$, can be decomposed into a series of eigenvalues, $\Lambda_i$, and their corresponding second-order orthonormal eigentensors, $\omega_i$:

$$ C_0 = \sum_i \Lambda_i \omega_i \otimes \omega_i \qquad (55) $$

The square root of the stiffness tensor, $C_0^{1/2}$, is calculated to facilitate the orthogonal decomposition:

$$ C_0^{1/2} = \sum_i \Lambda_i^{1/2} \omega_i \otimes \omega_i \qquad (56) $$
This process yields the pseudo-stress tensor, $\tilde{\epsilon}$, as follows:

$$ \tilde{\epsilon} = \mathbb{C}_0^{\frac{1}{2}} : \epsilon \qquad (57) $$

Leveraging the orthogonality condition between the positive and negative parts of the strain tensor, this method distinctly separates the pseudo-stress into tensile (positive values) and compressive (negative values) components:

$$ \tilde{\epsilon} = \tilde{\epsilon}^{+} + \tilde{\epsilon}^{-} \qquad (58) $$

The corresponding positive ($\sigma^{+}$) and negative ($\sigma^{-}$) stress tensors are thus formulated as:

$$ \sigma_0^+ = \mathbb{C}_0^{\frac{1}{2}} \tilde{\epsilon}^+, \quad \sigma_0^- = \mathbb{C}_0^{\frac{1}{2}} \tilde{\epsilon}^- \qquad (59) $$

Finally, the elastic strain energy positive and negative can be computed as:

$$ \psi_0^+(\epsilon) = \frac{1}{2} \sigma_0^+ : \epsilon, \quad \psi_0^-(\epsilon) = \frac{1}{2} \sigma_0^- : \epsilon \qquad (60) $$

And the stress tensor is obtained through:

$$ \sigma = g(d)\sigma_0^+ + \sigma_0^- \qquad (61) $$

## The cleavage plane based degradation method

In this section we will be using a new method based on increasing the appropriate compliance tensor components corresponding to the weakest cleavage plane as described in [17]. This approach to be used for anisotropic materials with an embedded weakest cleavage plane, induces degradation specifically along said plane while preserving integrity in other directions. Unlike the previous methods based on decomposition reliant on the deformation response, this approach relies on the material's inherent properties to guide damage evolution.

The stress acting along the weakest cleavage plane, defined by normal vector **m**, is:

$$ \sigma_m = \mathbf{m}^T \cdot \boldsymbol{\sigma} \cdot \mathbf{m} \qquad (62) $$

We define two projection tensors $\mathbb{N}$ and $\mathbb{M}$ as:

$$ \mathcal{N}_{ijkl} = m_i m_j m_k m_l \qquad (63) $$

$$ \mathcal{M}_{ijkl} = \frac{1}{4} (m_i m_k \delta_{jl} + m_j m_l \delta_{ik} + m_k m_l \delta_{il} + m_l m_i \delta_{ik}) \qquad (64) $$

While $\mathbb{N}$ captures normal components relative to the weakest cleavage plane, $\mathbb{M}$ encompasses off-diagonal components. The equivalent forms for Young's modulus ($E_n$) and shear modulus ($G_n$) are derived from the undamaged compliance tensor $S_0 = \mathbb{C}_0^{-1}$:

$$ \frac{1}{E_n} = \mathcal{N}_{ijkl} S_{0ijkl} \qquad (65) $$

$$ \frac{1}{G_n} = 2 \left( \mathcal{M}_{ijkl} S_{0ijkl} - \frac{1}{E_n} \right) \qquad (66) $$

with the use of the Einstein summation convention.

We define the tensor $\mathbb{H}$ that denotes the component that amplifies compliance in response to tensile stress ($\sigma_m > 0$):

$$ \mathbb{H} = \left( \frac{1}{E_n} - \frac{1}{G_n} \right) \mathbb{N} + \frac{\mathbb{M}}{G_n} \qquad (67) $$
Under compression ($\sigma_m < 0$), $\mathbb{H}$ accounts solely for shear stress contributions:

$$ \mathbb{H} = \frac{1}{G_n} (M - N) \qquad (68) $$

The compliance and elasticity tensors, updated to reflect damage, are:

$$ \mathbb{S} = \mathbb{S}_0 + g_s(d)\mathbb{H} \qquad (69) $$

The degradation of the material is modeled by an increase in the compliance tensor components using the function $g_s$:

$$ g_s(d) = \frac{1 - (1 - d)^2}{(1 - d)^2} = \frac{2d - d^2}{(1 - d)^2} \qquad (70) $$

$g_s$ must satisfy:

* $g_s(d = 0) = 0$, the material retains its original compliance tensor, ensuring full stiffness without any degradation in absence of damage.

* $g_s(d = 1) = \infty$, leading to infinite values of the compliance tensor components, thereby effectively nullifying any contribution of damaging energy once the material is fully damaged.

With the updated elasticity tensor

$$ \mathbb{C} = \mathbb{S}^{-1} \qquad (71) $$

the stress-strain relationship is thus represented as:

$$ \psi(\boldsymbol{\varepsilon}, d) = \frac{1}{2} \boldsymbol{\varepsilon} : \mathbb{C} : \boldsymbol{\varepsilon} \qquad (72) $$

The stress tensor is then given by:

$$ \boldsymbol{\sigma} = \mathbb{C} : \boldsymbol{\varepsilon} \qquad (73) $$

The work of [6] also leverages the phase field method to examine fracture propagation within crystalline structures, specifically focusing on crack propagation along the weakest cleavage plane. Their methodology involves penalizing deviations from planes that are not perpendicular to the normal vector of this weakest plane. Their method differs from the new cleavage plane-based degradation method presented in this article, as the penalization technique is embedded within the crack surface density function as follows:

$$ \gamma(d, \nabla d) = \frac{3d}{8l_c} + \frac{3l_c}{8} \nabla d \cdot \nabla d + \beta \frac{3l_c}{8} (\nabla d \cdot \nabla d - (\mathbf{m} \cdot \nabla d)^2) \qquad (74) $$

where $\beta$ represents the penalization factor, and $\mathbf{m}$ denotes the normal vector to the weakest cleavage plane. When $\beta \gg 0$, fracture on planes not normal to $\mathbf{m}$ are penalized.

The phase field evolution equation, modified to accommodate this framework, is:

$$ \dot{d} = -\frac{1}{\eta}(-2(1-d)\psi^+ + \frac{3}{8l_c}G_c - \frac{3l_c}{4}G_c\Delta d - \frac{3l_c}{4}G_c\beta(\Delta d - \mathbf{m} \cdot \nabla(\mathbf{m} \cdot \nabla d))) \qquad (75) $$

The penalization technique adopted in [6] is designed for use with conventional decomposition approaches such as symmetric, volumetric-deviatoric, or strain spectral methods. Consequently, it does not address some of the primary shortcomings of these methods, which will be discussed later in this work. In contrast, the new cleavage plane-based degradation method, integrated directly into the material characteristics, emphasizes that damage evolution occurs only along the predetermined weakest cleavage plane and avoids the use of traditional energy decomposition approaches coupled with a penalization technique.
## Governing balance equations

The model seeks solutions for the displacement field, $u$, and the damage variable, $d$, while adhering to the principles of mechanical equilibrium. The mechanical equilibrium equation is formulated as:

$$ \mathrm{div}\sigma + f = 0 \qquad (76) $$

where

$$ \sigma = \frac{\partial \psi(\boldsymbol{\epsilon}, d)}{\partial \boldsymbol{\epsilon}} \qquad (77) $$

The stress considered is a composite of the degraded stress component ($\sigma_0^+$) and the residual, non-degraded part ($\sigma_0^-$):

$$ \sigma = g(d) \frac{\partial \psi_0^+(\boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}} + \frac{\partial \psi_0^-(\boldsymbol{\epsilon})}{\partial \boldsymbol{\epsilon}} = g(d)\sigma_0^+ + \sigma_0^- \qquad (78) $$

For the phase field model, the damage variable's evolution follows the Allen-Cahn equation:

$$ \dot{d} = -L \left\langle \frac{\partial \psi}{\partial d}(\boldsymbol{\epsilon}, d) - G_c \delta_d \gamma \right\rangle_+ \qquad (79) $$

where the mobility, $L$ is function of $\eta$, the viscosity of the crack:

$$ L = \frac{1}{\eta} \qquad (80) $$

A coupled monolithic algorithm is employed to solve for both the displacement and the phase field variable simultaneously implemented in Idaho National Laboratory's open source FEM code Multiphysics Object Oriented Simulation Environment (MOOSE) [15].

## Numerical implementation

This investigation embarks on a comprehensive series of numerical tests to scrutinize the effects of different elastic energy decompositions on fracture mechanics simulations. With a constant set of geometry and material properties, the study contrasts the outcomes derived from six meticulously selected models, each embodying a unique approach to elastic energy decomposition:

1. Model 1: Symmetric formulation,

2. Model 2: Deviatoric/spherical decomposition of the strain tensor,

3. Model 3: Spectral decomposition of the strain tensor,

4. Model 4: Spectral decomposition of the stress tensor,

5. Model 5: Orthogonal decomposition of the strain tensor,

6. Model 6: Cleavage plane based degradation.

The essence of this endeavor is to meticulously compare and contrast the performances of these models under various loading conditions to illuminate their respective strengths and limitations in faithfully representing the material behavior under study. Through this numerical exploration, we aim to furnish a nuanced understanding of how different decomposition methods influence simulation accuracy and the predictive capabilities of phase field models in fracture mechanics.

## Tension and compression testing on a material (integration) point

This section delves into the examination of tension and compression responses using a single material point. We employ a 3D cubic element (HEX8) with dimensions of 1 × 1 × 1
<figure><img src="image_8.png" /><figcaption>Fig. 4 Tension and compression testing geometry and boundary conditions</figcaption></figure>

<table><caption>**Table 1** Material properties</caption><thead><tr><td>Material property</td><td>Value</td></tr></thead><tbody><tr><td>E</td><td>210 × 10<sup>3</sup> MPa</td></tr><tr><td>ν</td><td>0.3</td></tr><tr><td>G<sub>c</sub></td><td>2.7 N/mm</td></tr><tr><td>l<sub>c</sub></td><td>0.04 mm</td></tr><tr><td>η</td><td>1.0 × 10<sup>-3</sup> Ns/mm<sup>2</sup></td></tr></tbody></table>

mm for this purpose. Model dimensions and boundary conditions are illustrated in Fig. 4. Material properties pertinent to the study are summarized in Table 1. The experimental setup subjects the model to incremental tensile displacement initially, followed by compressive displacement, then another tensile displacement in a cyclic displacement loading as delineated in Fig. 5, applied along the top edge in the y-direction.

With imposed $\mathbf{u}_y$ at the top boundary (Fig. 4), the strain tensor is:

$$ \boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_x & 0 & 0 \\ 0 & \epsilon_y & 0 \\ 0 & 0 & \epsilon_z \end{bmatrix} \qquad (81) $$

<figure><img src="image_9.png" /><figcaption>Fig. 5 Prescribed displacement over time</figcaption></figure>
with $\epsilon_x = \epsilon_z = -\nu\epsilon_y$, and the stress tensor is:

$$ \sigma = \begin{bmatrix} 0 & 0 & 0 \\ 0 & \sigma_y = E\epsilon_y & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (82) $$

$\nu$ being the Poisson ratio and $E$ the Young modulus for homogeneous material.

### Model 1: Symmetric formulation

Model 1 operates under the premise that the elastic strain energy, in its entirety, contributes to damage evolution and is modulated by the degradation function. This model eschews any form of elastic energy decomposition, treating all strains as equally damaging. The outcome of the finite element numerical resolution is illustrated in Fig. 6. It shows the damage progression during tensile loading and, notably during the compression phase, a marginal increase in compression stress is observed which contradicts expectations, as crack closure under compressive strains should ostensibly amplify stress—yet this anticipated response remains absent in the numerical results.

When $d = 0$, this uniaxial tensile loading yields:

$$ \psi(\boldsymbol{\epsilon}, d) = g(d)\psi_0(\boldsymbol{\epsilon}) = g(d)\frac{1}{2}\boldsymbol{\epsilon} : \mathbb{C}_0 : \boldsymbol{\epsilon} = \frac{1}{2}E\boldsymbol{\epsilon}_y^2 \qquad (83) $$

since $g(d = 0) = 1$. Upon the damage variable nearing $g(d = 1) = 0$, the formula simplifies to:

$$ \psi(\boldsymbol{\epsilon}, d) = g(d)\psi_0(\boldsymbol{\epsilon}) = 0 \qquad (84) $$

It is noteworthy that once the damage variable reaches a high enough value, there is no further contribution to the damaging energy, but in instances where the damage is still in its benign stage, potential damage evolution would have been observable alongside compression modes, provided the critical strain energy value had been attained.

### Model 2: Deviatoric/spherical decomposition of strain tensor

Model 2 employs a deviatoric/spherical decomposition of the strain tensor proposed by [2], to partition elastic strain energy into damaging and non-damaging components. As shown in Fig. 7, initial tensile loading induces an elastic response, with damage onset. Upon the application of compressive forces, the model appears to revert to its original state; however, the negligible stress values obtained contradict with the expectation of a clear recovery of stiffness and a greater stress due to the closure of cracks under compression.

The resulting strain tensor $\boldsymbol{\epsilon}$, considering linear isotropic elasticity, is segmented into deviatoric ($\boldsymbol{\epsilon}_d$) and spherical components ($\boldsymbol{\epsilon}_s$).

$$ \boldsymbol{\epsilon} = \boldsymbol{\epsilon}_d + \boldsymbol{\epsilon}_s \qquad (85) $$

$$ \boldsymbol{\epsilon} = \begin{bmatrix} -\frac{1}{3}(\nu + 1)\epsilon_y & 0 & 0 \\ 0 & \frac{2}{3}(\nu + 1)\epsilon_y & 0 \\ 0 & 0 & -\frac{1}{3}(\nu + 1)\epsilon_y \end{bmatrix} + \begin{bmatrix} \frac{1}{3}(1 - 2\nu)\epsilon_y & 0 & 0 \\ 0 & \frac{1}{3}(1 - 2\nu)\epsilon_y & 0 \\ 0 & 0 & \frac{1}{3}(1 - 2\nu)\epsilon_y \end{bmatrix} \qquad (86) $$

The decomposition leads to the definition of the damaging and non-damaging contributions to elastic strain energy as in (40). For $\epsilon_y > 0$, we obtain:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{1}{2}E\epsilon_y^2 \qquad (87) $$
<figure><img src="image_10.png" /><figcaption>Fig. 6 Tension and compression testing with symmetric formulation</figcaption></figure>

$$ \psi_0^-(\boldsymbol{\epsilon}) = 0 \qquad (88) $$

Upon damage reaching 1, only non-damaging components contribute, implying no further energy addition or damage progression:

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_0^-(\boldsymbol{\epsilon}) = 0 \qquad (89) $$

For compressive loading with $\epsilon_y < 0$ and damage at 1:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{E\epsilon_y^2(\nu + 1)}{3} \qquad (90) $$

$$ \psi_0^-(\boldsymbol{\epsilon}) = \frac{E\epsilon_y^2(1 - 2\nu)}{6} \qquad (91) $$

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_0^-(\boldsymbol{\epsilon}) = 0.067E\epsilon_y^2 \qquad (92) $$

Figure 8 presents the outcomes of a loading scenario wherein the damage remains markedly below 1 following tensile loading. It is remarkable that the stress evolution at the outset of compressive loading is characterized by a slope significantly inferior to the original (prior to damaging). This observation underlines the insufficiency of the decomposition approach to restore initial stiffness levels, thereby indicating the inadequacy of the model in effectively simulating crack closure phenomena under compressive forces. Furthermore, an intriguing development of damage under compression is observed, culminating in the stress diminishing to zero. This analysis underscores the complex interplay between damage evolution and mechanical response, particularly in the context of compressive strain effects.

In summary, the application of these first two models necessitates prudence, particularly in scenarios characterized by cyclic loading or dynamic testing conditions leading to compressive strains.
<figure><img src="image_11.png" /><figcaption>Fig. 7 Tension and compression testing with deviatoric/spherical decomposition of the strain tensor</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Tension and compression testing with deviatoric/spherical decomposition of the strain tensor—not fully damaged</figcaption></figure>
### Model 3: Spectral decomposition of the strain tensor

With this spectral decomposition of the strain tensor, from the Eqs. (44) and (45) we obtain:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{E\epsilon_y^2}{2(\nu + 1)}(1 - \nu(1 - 2\nu)) \quad (93) $$

$$ \psi_0^-(\boldsymbol{\epsilon}) = \frac{E\nu^2\epsilon_y^2}{2(\nu + 1)} \quad (94) $$

For a tensile loading. Unlike the spheric/deviatoric model, the spectral approach evaluates both positive and negative values of the strain tensor's trace, as well as individual diagonal terms, offering a nuanced understanding of strain contributions to damage.

Upon damage escalation to 1, only the non-damaging energy component remains significant:

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_0^-(\boldsymbol{\epsilon}) = 0.07E\epsilon_y^2 \quad (95) $$

Considering compression with $\epsilon_y < 0$ and $d = 1$:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{E\nu^2\epsilon_y^2}{2(\nu + 1)} \quad (96) $$

$$ \psi_0^-(\boldsymbol{\epsilon}) = \frac{E\epsilon_y^2}{2(\nu + 1)}(1 - \nu(1 - 2\nu)) \quad (97) $$

The strain energy under compression at $d = 1$ is dominated by the non-damaging component, $\psi^-$:

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_0^-(\boldsymbol{\epsilon}) = 0.43E\epsilon_y^2 \quad (98) $$

The resultant stress, therefore, reflects only the non-damaging stress, $\sigma_0^-$:

$$ \sigma = \sigma_0^- \quad (99) $$

Figure 9 illustrates the stress-strain relationship for Model 3. Due to $d \neq 0$ and the degradation of $\psi_0^+$, the slope is slightly less steep but barely noticeable as $\psi_0^+$ is already overshadowed by the value of $\psi_0^-$ in the case of compressive strains. The stiffness recovery in compressive strain is important but not complete.

### Model 4: Spectral decomposition of the stress tensor

Model 4 results are shown in Fig. 10. Upon transitioning to compression loading, the stress-strain slope mirrors the one observed during at the beginning of tension loading.

Elastic strain energy decomposition in this model, under initial tensile loading ($\sigma_y > 0$) with damage $d = 0$, is determined by the stress tensor:

$$ \sigma = \begin{bmatrix} 0 & 0 & 0 \\ 0 & \sigma_y & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad (100) $$

The model emphasizes the separation of stress tensor eigenvalues into positive and negative domains, with a uniaxial tension test predominantly engaging the positive component, and a compression test favoring the negative component. Consequently, for $\epsilon_y > 0$:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{1}{2}E\epsilon_y^2 = \frac{1}{2}\frac{\sigma_y^2}{E}, \quad \psi_0^-(\boldsymbol{\epsilon}) = 0 \quad (101) $$
<figure><img src="image_13.png" /><figcaption>**Fig. 9** Tension and compression testing with strain spectral decomposition</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Tension and compression testing with stress spectral decomposition</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11 Tension and compression testing with strain spectral decomposition—not fully damaged</figcaption></figure>

Conversely, for $\epsilon_y < 0$:

$$ \psi_0^+(\boldsymbol{\epsilon}) = 0, \quad \psi_0^-(\boldsymbol{\epsilon}) = \frac{1}{2} \frac{\sigma_y^2}{E} \qquad (102) $$

Distinctly, this model does not facilitate damage progression under compressive forces, diverging from the strain-based spectral method that allows for potential damage development after exceeding a critical strain value. Figure 11 showcases an instance where the material does not attain a fully damaged state during tensile loading but when reaching a critical value later during compression loading, resumes damage progression for the strain spectral decomposition, which is halted when tensile loading is applied again, only to continue evolving when $\psi_0^+(\boldsymbol{\epsilon})$ reaches its highest value attained in the previous tensile loading. This is not the case for the stress spectral method shown in Fig. 12 where damage evolution halts upon the application of compressive loading and remains stagnant throughout the rest of the test. This distinction underscores the model's reliance on stress tensor decomposition, which inherently disregards transverse dilatation and positive strain effects on damage enhancement.

### Model 5: Orthogonal decomposition of the strain tensor

Model 5 introduces the novel approach of orthogonal decomposition of the strain tensor, a method yet to be extensively explored in the context of tension and compression loading. As depicted in Fig. 13, the stress response under sustained compression does not lead to recovery of the initial stiffness of the undamaged material. This divergence is attributed to the spectral decomposition into positive and negative terms of the pseudo stress tensor
<figure><img src="image_16.png" /><figcaption>Fig. 12 Tension and compression testing with stress spectral decomposition—not fully damaged</figcaption></figure>

obtained through computation using the square root of the original stiffness tensor (C₀). The pseudo-stress $\tilde{\varepsilon} = C_0^{1/2} : \varepsilon$ is calculated as:

$$ \tilde{\varepsilon} = \begin{bmatrix} \frac{\sqrt{\varepsilon_{\varepsilon_y}}(1-2\nu)\sqrt{-\frac{\nu+1}{2\nu-1}-(\nu+1)}}{3\sqrt{\nu+1}} & 0 & 0 \\ 0 & \frac{\sqrt{\varepsilon_{\varepsilon_y}}(1-2\nu)\sqrt{-\frac{\nu+1}{2\nu-1}+2(\nu+1)}}{3\sqrt{\nu+1}} & 0 \\ 0 & 0 & \frac{\sqrt{\varepsilon_{\varepsilon_y}}(1-2\nu)\sqrt{-\frac{\nu+1}{2\nu-1}-(\nu+1)}}{3\sqrt{\nu+1}} \end{bmatrix} \quad (103) $$

$$ \tilde{\varepsilon} = \begin{bmatrix} -0.17\sqrt{\varepsilon_{\varepsilon_y}} & 0 & 0 \\ 0 & 0.971\sqrt{\varepsilon_{\varepsilon_y}} & 0 \\ 0 & 0 & -0.17\sqrt{\varepsilon_{\varepsilon_y}} \end{bmatrix} \quad (104) $$

The positive and negative components of the strain energy, dependent on the sign of $\tilde{\varepsilon}_i$, are mirrored in tension and compression scenarios. Consequently, the resulting elastic strain energy components in tension ($\varepsilon_y > 0$) are:

$$ \psi_0^-(\boldsymbol{\epsilon}) = \frac{E \varepsilon_y^2 \left( 0.44 \nu^2 \sqrt{-\frac{1}{2\nu-1}} + 0.22 \nu \sqrt{-\frac{1}{2\nu-1}} - 0.11 \nu \sqrt{\nu+1} - 0.22 \sqrt{-\frac{1}{2\nu-1}} + 0.22 \sqrt{\nu+1} \right)}{\sqrt{\nu+1}} \quad (105) $$

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{E \varepsilon_y^2 \left( -0.44 \nu^2 \sqrt{-\frac{1}{2\nu-1}} - 0.22 \nu \sqrt{-\frac{1}{2\nu-1}} + 0.11 \nu \sqrt{\nu+1} + 0.22 \sqrt{-\frac{1}{2\nu-1}} + 0.277 \sqrt{\nu+1} \right)}{\sqrt{\nu+1}} \quad (106) $$

$$ \psi_0^+(\boldsymbol{\epsilon}) = 0.47 \varepsilon_y^2, \quad \psi_0^-(\boldsymbol{\epsilon}) = 0.029 \varepsilon_y^2 \quad (107) $$

In compression ($\varepsilon_y < 0$), the roles of $\psi^+$ and $\psi^-$ are inverted:

$$ \psi_0^+(\boldsymbol{\epsilon}) = 0.029 E \varepsilon_y^2, \quad \psi_0^-(\boldsymbol{\epsilon}) = 0.47 E \varepsilon_y^2 \quad (108) $$

This model predicts a similar behavior to the strain spectral method, with Fig. 13 illustrating the stress–strain relationship. The slope of the stress-negative displacement is slightly less steep than the initial tensile phase. Damage evolution under compressive loading is possible due to the nonzero $\psi_0^+$ component.
<figure><img src="image_17.png" /><figcaption>**Fig. 13** Tension and compression testing with orthogonal decomposition</figcaption></figure>

**Model 6: Cleavage plane based degradation**

Model 6 explores a novel approach by selectively degrading components associated with the weakest cleavage plane, characterized by the normal vector **m** = (0, 1, 0). The response of this model, as illustrated in Fig. 14, embarks with an elastic phase during tensile loading, followed by damage evolution that escalates as stresses diminish to zero. This model shows a complete recovery of the stiffness upon the application of compressive forces, which upholds the correct representation of crack closure. This behavior underscores the critical role of **H** in the model's response under compression.

In compression, for **m** = (0, 1, 0), E reverts to *E*₀, mimicking the behavior of an undamaged material under compressive stress. This model would lead to the progression of damage when shear stresses are present which is not the case for this loading scenario since *ε*<sub>shear</sub> = 0.

**Pure shear ton a material (integration) point**

In this section, we detail a numerical implementation of a pure shear test on a single material point with the same material characteristics as previously defined in Table 1, and with the geometry and boundary conditions defined in the Fig. 15. The forced displacement at the top is set to $\Delta u_x = 10^{-5}$ mm per time increment ($\Delta t = 10^{-5}$s).

**Model 1: Symmetric formulation**

The symmetric formulation in Model 1 degrades the entirety of the elastic strain energy and this is very clearly demonstrated in the Fig. 16. As damage progresses, the stress diminishes to zero as expected.
<figure><img src="image_18.png" /><figcaption>**Fig. 14** Tension and compression testing with cleavage plane based degradation method</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 15** Geometry and boundary conditions for the shear test on a material (integration) point</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 16** Shear test with Model 1 and Model 2</figcaption></figure>
### Model 2: Deviatoric/spherical decomposition of the strain tensor

Model 2 which employs the deviatoric/spherical decomposition of the strain tensor of [2], considers all deviatoric strains to be damaging. For shear loading with $d = 0$:

$$ \boldsymbol{\epsilon} = \begin{bmatrix} 0 & \epsilon_{xy} & 0 \\ \epsilon_{xy} & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (109) $$

The decomposition defines the damaging and non-damaging contributions to elastic strain energy as:

$$ \psi_{0}^{+}(\boldsymbol{\epsilon}) = \frac{E}{\nu + 1} \epsilon_{xy}^{2}, \quad \psi_{0}^{-}(\boldsymbol{\epsilon}) = 0 \qquad (110) $$

And when the damage is around 1, only the residual energy contributes to the stress calculation:

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_{0}^{-}(\boldsymbol{\epsilon}) = 0 \qquad (111) $$

Which aligns perfectly with the results obtained in the Fig. 16 also similar to the results of the previous method.

### Model 3: Spectral decomposition of the strain tensor

Model 3 considers a spectral decomposition of the strain tensor, the division only considers the sign of the eigenvalues for the decomposition so when shear stresses are applied, the elastic strain energy is equally contributing to the damaging part of the total strain energy ($\psi_0^+$) and the residual part ($\psi_0^-$).

$$ \psi_{0}^{+}(\boldsymbol{\epsilon}) = \frac{E}{2(\nu + 1)} \epsilon_{xy}^{2} \qquad (112) $$

$$ \psi_{0}^{-}(\boldsymbol{\epsilon}) = \frac{E}{2(\nu + 1)} \epsilon_{xy}^{2} \qquad (113) $$

At the onset of the test ($d = 0$):

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_{0}^{+}(\boldsymbol{\epsilon}) + \psi_{0}^{-}(\boldsymbol{\epsilon}) = \frac{E}{1.3} \epsilon_{xy}^{2} \qquad (114) $$

And for $d = 1$:

$$ \psi(\boldsymbol{\epsilon}, d) = \psi_{0}^{-}(\boldsymbol{\epsilon}) = \frac{E}{2.6} \epsilon_{xy}^{2} \qquad (115) $$

This shows that even after full damage, there is still a significant non-degraded part of the elastic strain energy leading to a stress increase, just as depicted in the results presented in Fig. 17.

It is also noticeable that the elastic stage lasts longer than the previous two models, this is because $\psi_0^+$ in Model 3 takes half the value obtained with Model 2 or Model 1 and thus, necessitates higher stresses to reach the critical strain energy needed for the damage initiation.

### Model 4: Spectral decomposition of the stress tensor

Model 4, just like Model 3 is a spectral decomposition, separating the damaging elastic strain energy from the residual depending on the sign of the eigenvalues of the stress tensor. In a pure shear test, when expressed in the principal direction basis, the strain tensor is given by:

$$ \boldsymbol{\epsilon} = \begin{bmatrix} \epsilon_{xy} & 0 & 0 \\ 0 & -\epsilon_{xy} & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (116) $$
<figure><img src="image_21.png" /><figcaption>Fig. 17 Shear test with Model 3, Model 4, and Model 5</figcaption></figure>

The stress tensor is then:

$$ \sigma = \begin{bmatrix} \sigma_{xy} & 0 & 0 \\ 0 & -\sigma_{xy} & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (117) $$

with $\sigma_{xy} = \frac{E}{1+\nu}\epsilon_{xy}$. Thus, the elastic strain energy decomposition writes:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{1}{2}\sigma_{xy}\epsilon_{xy} = \frac{E}{2(\nu+1)}\epsilon_{xy}^2 \quad \text{and} \quad \psi_0^-(\boldsymbol{\epsilon}) = \frac{1}{2}\sigma_{xy}\epsilon_{xy} = \frac{E}{2(\nu+1)}\epsilon_{xy}^2 \qquad (118) $$

just like in the results obtained with Model 3 decomposition method (Fig. 17).

### Model 5: Orthogonal decomposition of the strain tensor

Model 5 also follows in the footsteps on the previous two spectral decomposition (Model 3 and 4), this one decomposes the pseudo-stress obtained from the square root of the elasticity tensor:

$$ \tilde{\boldsymbol{\epsilon}} = \mathbb{C}_0^{\frac{1}{2}} : \boldsymbol{\epsilon} \qquad (119) $$

Using the previously defined $\mathbb{C}_0^{\frac{1}{2}}$ in, the resulting pseudo-stress $\tilde{\boldsymbol{\epsilon}}$:

$$ \tilde{\boldsymbol{\epsilon}} = \begin{bmatrix} \sqrt{\frac{E}{2(1+\nu)}} \epsilon_{xy} & 0 & 0 \\ 0 & -\sqrt{\frac{E}{2(1+\nu)}} \epsilon_{xy} & 0 \\ 0 & 0 & 0 \end{bmatrix} \qquad (120) $$

And according to Eq. (53), we write:

$$ \psi_0^+(\boldsymbol{\epsilon}) = \frac{E}{2(\nu+1)}\epsilon_{xy}^2 \quad \text{and} \quad \psi_0^-(\boldsymbol{\epsilon}) = \frac{E}{2(\nu+1)}\epsilon_{xy}^2 \qquad (121) $$

These results show that the stress can still be increasing even after $d = 1$ due to the non-zero $\psi_0^-$ component, which matches the results obtained in Fig. 17.

### Model 6: Cleavage plane based degradation

Model 6 uses a method which degrades the components associated with the weakest cleavage plane. In this test we consider $\mathbf{m} = (0, 1, 0)$ the normal to the cleavage plane. Fig-
<figure><img src="image_22.png" /><figcaption>**Fig. 18** Shear test with cleavage plane based degradation</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 19** Geometry, boundary conditions of the 2D shear test model</figcaption></figure>

ure 18 shows the pure shear test results that match with those obtained when considering Models 1 and 2 (Fig. 16).

We can see that the stresses diminish to 0 as the damage increase, which is what is expected as a result of a pure shear test.

## Pure shear test on a notched edge model

In this section, we focus on the shear test of a square notched specimen in plane strain formulation. This test has previously been explored by [5, 18]. The geometry and boundary conditions of pre-cracked specimen are shown in Fig. 19. The mesh is made up of 14,400 QUAD elements and as for the previous numerical implementations, a coupled monolithic algorithm is used to solve for the displacement vector and phase field variable. Utilizing the same decomposition methods and material properties as previously discussed, except for the characteristic length which is chosen to be $l_c = 10^{-5}$ m, allows for a focused comparison under these specific conditions.
Figure 20 presents the force-displacement curve for all six models, where a consistent response is observed in the elastic stage across all methodologies, characterized by a steady increase in stress. This uniformity diverges as damage approaches a critical level at the crack tip, as depicted in Fig. 21.

Model 1's symmetric formulation results in a stress reduction as the crack propagates, culminating near zero as the crack bisects the structure (Fig. 21c). This model's complete degradation of elastic strain energy is mirrored in the stress reduction correlating with damage increase.

Model 2, employing deviatoric/spherical decomposition of the strain tensor, displays a precipitous stress fall as the crack fully traverses the structure. This method's distinction with the first one lies in its exclusion of negative hydrostatic stress values from damage calculation, affecting the stress-strain curve's sharper decline and influencing the crack path's direction, contrasting with Model 1's symmetric bifurcation. For Model 2, a part of the spherical strains (the negative values representing compression modes) doesn't account for damage evolution and the results represented in Fig. 21f shows that some elements appear to be under compression leading to only a downward curved crack path.

Model 3, based on [18]'s strain spectral decomposition, shows a stress decline initiation that abruptly ceases despite the crack's complete passage (Fig. 21i), contradicting expected behavior under full structural bisecting by a crack. This model's division of shear-induced strain energy into damaging and non-damaging components yields an unreliable outcome for pure shear loading, as when the damage is equal to 1, $\psi(\epsilon, d) = \bar{\psi}_0(\epsilon)$ which isn't null and will lead to an increase in the stresses after full damage. Similar results can be seen in the application of the other spectral decomposition methods, Model 4 with stress spectral decomposition (Fig. 21j-l) and Model 5 with the orthogonal decomposition (Fig. 21j-n, p). Therefore, it's advisable to exercise caution when employing these decomposition methods in scenarios involving shear loading, particularly in cases involving polycrystalline materials where normal loading may generate shear stresses within the material due to its heterogeneous microstructure, potentially triggering damage evolution, and as demonstrated, the outcomes yielded by these decomposition methods may not be suitable for such circumstances.

Conversely, Model 6's cleavage plane based degradation method results in a pronounced stiffness loss, evidenced by stress diminishing to negligible levels upon the crack extending across the entire structure, as seen in Figs. 21q-s. The model exhibits a higher initial stress peak, stemming from an extended elastic phase caused by the alteration in the critical elastic strain energy required for damage initiation. But once triggered, the damage evolution and propagation of the crack happen very quickly. Unlike other outcomes, the crack path follows a horizontal trajectory, attributed to the cleavage plane specified in this model, where its normal vector is (0, 1, 0). A crack path for a normal vector of (0.707, 0.707, 0) would have a crack path similar to the previous Models as can be seen in the Fig. 22. This method is interesting as it's only penalizing the weakest cleavage plane and not all directions in the structure, so for a direction of (0.707, 0.707, 0) there would still be a substantive force residual as seen in the Fig. 23, considering that the damage is not affecting the entirety of the stresses pertaining to the uniaxial loading direction, but only those affected to the direction of the normal to the weakest cleavage plane. This is interesting when regarding polycrystalline material where the material's inherent anisotropy could suggest direction dependent damage evolution. Experimental observation have shown
<figure><img src="image_24.png" /><figcaption>Fig. 20 Force-displacement curves for the shear tests</figcaption></figure>

<table><caption>**Table 2** Material properties</caption><thead><tr><td>Material property</td><td>Value</td></tr></thead><tbody><tr><td>E</td><td>210 × 10<sup>3</sup> MPa</td></tr><tr><td>ν</td><td>0.3</td></tr><tr><td>G<sub>c</sub></td><td>2.7 N/mm</td></tr><tr><td>l<sub>c</sub></td><td>0.02 mm</td></tr><tr><td>η</td><td>5.0 × 10<sup>-4</sup> Ns/mm<sup>2</sup></td></tr></tbody></table>

crystalline material to have a preferred plane along which cleavage fracture occurs such as the plane {111} for silicon [11], {100} for halite [16], {0001} for zinc [13].

## Compression and tension test on holed square

In this section, we examine a two-dimensional square with dimensions of 1 mm by 1 mm, containing a central circular hole with a radius of 0.2 mm. The mesh is composed of 56,576 linear triangular elements. The square is subjected to compression then tension loading as defined in Fig. 25 and the geometry and mesh can be found in Fig. 24. The material properties used for the calculations are in Table 2. The vector normal to the weakest cleavage plane in Model 6 is (0, 1, 0).

Figure 26 presents the force-displacement curve for all six models and Fig. 27 shows the crack propagation path.

The damage evolution in Models 1 and 2 under compressive loading results in significant material degradation, ultimately leading to complete structural failure. The crack propagates through the entire structure, eliminating residual stiffness and leaving no remaining load-bearing capacity. During tensile loading, this degradation persists, as the initial stiffness is degraded, and no stiffness recovery occurs under tensile modes.

Model 3, using the spectral decomposition of the strain tensor, demonstrates damage evolution under compressive loading but with significantly reduced extent compared to Models 1 and 2. This degradation manifests as a minor reduction in stiffness during tensile loading, evidenced by a slight decrease in the slope of the stress–strain curve. With continued tensile loading, the crack propagates further, leading to a stress reduction that eventually reaches zero.
Model 4 exhibits no damage evolution under compressive loading, thereby maintaining consistent stiffness throughout both compressive and tensile phases. However, upon reaching the critical elastic strain energy during tensile loading, crack propagation begins, eventually leading to a loss of load-bearing capacity.

Model 5 exhibits behavior similar to Model 3, with minor damage evolution under compressive loading. This stiffness is very slightly degraded and this degradation persists into the tensile phase, where the stress–strain slope is slightly diminished. As the crack propagates, stiffness is fully lost, causing stress to drop to negligible levels.

For Model 6, we have considered a weakest cleavage plane with a normal vector aligned with the y-axis (0, 1, 0). This energy decomposition is based on the stress tensor acting on the cleavage plane, as articulated in Eqs. 67 and 68. As illustrated in Fig. 28, the damage progression starts under compressive loads owing to the significant damaging shear stress

<figure><img src="image_25.png" /><figcaption>**Fig. 21** Comparison of crack paths in shear test for the six models</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 21 continued</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 22 Model 6 crack path for $m = (0.707, 0.707, 0)$</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Fig. 23 Force-displacement curves for the Model 6 considering two different weak cleavage planes</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 24 Compression-tension test on a holed square: geometry, boundary condition, and mesh</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 25 Prescribed displacement over time</figcaption></figure>
<figure><img src="image_31.png" /><figcaption>**Fig. 26** Force-strain curves for compression-tension tests on a holed square</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 27** Comparison of the different crack paths of the compression-tension test</figcaption></figure>

observed at the midspan between the vertical and horizontal edges. The evolution of damage under these compressive loads exhibits a distinct pattern, manifesting as a vertical sequence of short horizontal cracks. When subjected to tensile loads, the stress escalates until it attains the peak load-bearing capability of the partially damaged structure. Crack propagation then proceeds horizontally along the peripheries of the hole, with the stress subsequently diminishing to zero, as depicted in Fig. 27.
<figure><img src="image_33.png" /><figcaption>Fig. 28 Compression-tension test with Model 6: stress, strain and damage evolution</figcaption></figure>

## Conclusion

Throughout this comprehensive study, we have meticulously explored the evolution of phase field modeling in the context of brittle fracture, with a particular focus on the role of elastic strain energy decomposition methods. Our investigation sheds light on the limitations of existing models under certain loading scenarios and identifies models that demonstrate superior performance. All methods have shown to be reliable in the case of tensile loading, but some have failed in other loading scenarios. The strain volumetric-deviatoric decomposition's limitations are evident in compressive testing on pre-damaged materials, urging caution in its application, particularly under cyclic loading conditions. Spectral decomposition methods demand meticulous scrutiny in scenarios where shear stress plays a critical role, particularly in the context of polycrystalline materials with strongly anisotropic properties or preferred crystallographic orientations. Among the methods assessed, the cleavage plane based degradation method for anisotropic mate-
rial with embedded weakest cleavage plane stands out for its effectiveness, consistently delivering accurate results across various fracture modes and loading conditions. This distinction is especially relevant given the relative scarcity of phase field literature addressing these complex scenarios. The cleavage plane based degradation model, by virtue of its adaptability and precision, offers significant promise for advancing our understanding of fracture mechanics in polycrystalline materials. It provides a potent framework for capturing the nuanced behaviors inherent to brittle fracture, paving the way for more refined and accurate predictive models. As such, continued research into and application of this model represent vital steps forward, with the potential to enhance both theoretical insights and practical applications within the field of fracture mechanics.

**Author contributions**

All authors contributed to the writing and reviewing of the article.

**Funding**

This research is supported by Angers Loire Métropole.

**Data availability**

All data and material are available upon request.

**Declarations**

**Competing interests**

No competing interests to be declared.

Received: 19 September 2024 Accepted: 20 March 2025

Published online: 22 April 2025

**References**

1. Allen SM, Cahn JW. A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening. Acta Metal. 1979;27(6):1085-95.

2. Amor H, Marigo JJ, Maurini C. Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids. 2009;57(8):1209-29.

3. Balay S, Abhyankar S, Adams MF, Brown J, Brune P, Buschelman K, Dalcin L, Dener A, Eijkhout V, Gropp W, et al. Petsc users manual (rev. 3.13). Technical report, Argonne National Lab.(ANL), Argonne, IL (United States); 2020.

4. Bazant PZ, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. J Appl Mech. 1988. https://doi.org/10.1115/1.3173674.

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

6. Clayton JD, Knap J. Phase field modeling of directional fracture in anisotropic polycrystals. Comput Mater Sci. 2015;98(C):158-69.

7. Ferreira AR, Marengo A, Perego U. A phase-field gradient-based energy split for the modeling of brittle fracture under load reversal. Comput Methods Appl Mech Eng. 2024;431: 117328.

8. Francfort G, Marigo JJ. Revisiting brittle fracture as an energy minimization problem. J Mech Phys Solids. 1998;46(8):1319-42.

9. Gmati H, Mareau C, Ammar A, El Arem S. A phase-field model for brittle fracture of anisotropic materials. Int J Numer Methods Eng. 2020;121:03.

10. Griffith AA, Geoffrey T, Ingram VI. The phenomena of rupture and flow in solids. Philos Trans R Soc Lond Ser A. 1921;221(582-593):163-98.

11. Hauch JA, Holland D, Marder MP, Swinney HL. Dynamic fracture in single crystal silicon. Phys Rev Lett. 1999;82(19):3823.

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

13. Kamdar MH. Cleavage in zinc. Metal Trans. 1971;2:485-9.

14. Lancioni G, Royer-Carfagni G. The variational approach to fracture mechanics. A practical application to the French panthéon in Paris. J Elasticity. 2009;95:1-30.

15. Lindsay AD, Gaston DR, Permann CJ, Miller JM, Andrš D, Slaughter AE, Kong F, Hansel J, Carlsen RW, Icenhour C, Harbour L, Giudicelli GL, Stogner RH, German P, Badger J, Biswas S, Chapuis L, Green C, Hales J, Hu T, Jiang W, Jung YS, Matthews C, Miao Y, Novak A, Peterson JW, Prince ZM, Rovinelli A, Schunert S, Schwen D, Spencer BW, Veeraraghavan S, Recuero A, Yushu D, Wang Y, Wilkins A, Wong C. 2.0-MOOSE: enabling massively parallel multiphysics simulation. SoftwareX. 2022;20: 101202.

16. Mansouri H, Ajalloeian R. Mechanical behavior of salt rock under uniaxial compression and creep tests. Int J Rock Mech Min Sci. 2018;110:19-27.
17. Mareau C, Morel F. A continuum damage mechanics-based approach for the high cycle fatigue behavior of metallic polycrystals. Int J Damage Mech. 2019;28(6):838-56.

18. Miehe C, Hofacker M, Welschinger F. A phase field model for rate-independent crack propagation: robust algorithmic implementation based on operator splits. Comput Methods Appl Mech Eng. 2010;199(45):2765-78.

19. Miehe C, Lambrecht M. Algorithms for computation of stresses and elasticity moduli in terms of Seth-hill's family of generalized strain tensors. Commun Numer Methods Eng. 2001;17(5):337-53.

20. Miehe C, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int J Numer Methods Eng. 2010;83(10):1273-311.

21. Moës N, Belytschko T. X-fem, de nouvelles frontières pour les éléments finis. Eur J Comput Mech. 2002;11:305-18.

22. Nguyen TT, Yvonnet J, Waldmann D, He QC. Implementation of a new strain split to model unilateral contact within the phase field method. Int J Numer Methods Eng. 2020;121(21):4717-33.

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

24. Steinke C, Kaliske M. A phase-field crack model based on directional stress decomposition. Comput Mech. 2019;63:1019-46.

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

26. Vicentini F, Zolesi C, Carrara P, Maurini C, De Lorenzis L. On the energy decomposition in variational phase-field models for brittle fracture under multi-axial stress states. Int J Fract. 2024. https://doi.org/10.1007/s10704-024-00763-w.

27. Wu JY. A unified phase-field theory for the mechanics of damage and quasi-brittle failure. J Mech Phys Solids. 2017;103:72-99.

28. Zhang S, Jiang W, Tonks MR. A new phase field fracture model for brittle materials that accounts for elastic anisotropy. Comput Methods Appl Mech Eng. 2020;358: 112643.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.