# Path-following methods for phase-field simulation of quasi-static crack propagation

Flavien Loiseau, Véronique Lazarus

► To cite this version:

Flavien Loiseau, Véronique Lazarus. Path-following methods for phase-field simulation of quasi-static
crack propagation. 2026. hal-05478246

HAL Id: hal-05478246

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

Preprint submitted on 26 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Highlights

## Path-following methods for phase-field simulation of quasi-static crack propagation

Flavien Loiseau ⚙️, Véronique Lazarus ⚙️

* Identification of generic path-following methods compatible with staggered solvers.

* A new method based on the strain outside the cracked zone offers improved performance.

* Evaluation by comparison with increasingly complex linear elastic sharp crack problems.

* Toward numerically robust, self-consistent evaluation of equilibrium crack propagation paths.

* Toward rational design of microstructures for enhanced fracture toughness.
# Path-following methods for phase-field simulation of quasi-static crack propagation

Flavien Loiseau Ⓔ<sup>a,*</sup>, Véronique Lazarus Ⓔ<sup>a</sup>

<sup>a</sup>IMSIA, CNRS, EDF, ENSTA, Institut Polytechnique de Paris, 91120, Palaiseau, France

## Abstract

The variational approach to fracture, particularly through its regularization as a phase-field model, has become a widely used tool for simulating the quasi-static propagation of cracks in structures. However, classic incremental loading can induce unstable crack growth, violating the quasi-static assumption, and in some cases, leads to a loss of force balance, preventing the estimation of dissipated energy during snapback instabilities. To address this challenge, path-following methods are investigated. Their aim is to adjust the applied load so that it stays at the propagation threshold, thereby preserving the quasi-static assumption and ensuring equilibrium solutions. In this work, we apply and evaluate multiple path-following methods within the framework of variational phase-field fracture models, which are developed to regularize linear elastic variational sharp crack evolution problems. Our study pursues two objectives. First, we identify several path-following methods that are both constitutive model-independent and problem-independent, while remaining straightforward to implement. To achieve this, we focus on partitioned path-following methods based on the displacement field, which decouple the path-following control equation from the rest of the system, providing an easier integration into staggered solvers. In addition, we also introduce a new path-following method which limit the maximum strain increment outside the cracked regions. Second, through three crack propagation problems of increasing complexity, we assess the equilibrium path predictions of the variational phase-field model by comparison with the associated sharp crack model. The comparison demonstrates that the proposed path-following method offers a simple yet highly effective approach to capture the equilibrium path in phase-field fracture simulations. This method robustly maintains the quasi-static assumption, ensuring physically meaningful results. By enabling accurate estimation of the energy dissipated during snapback instabilities, it paves the way for the rational design of more resistant heterogeneous materials.

**Keywords:** Path-following methods, Linear Elastic Fracture Mechanics, Variational Griffith phase-field model, Quasi-static crack propagation, Snap-back instability.

## 1. Introduction

To ensure the safety and integrity of mechanical structures, particularly those with pre-existing defects, understanding and predicting crack propagation is crucial. To this aim, Griffith (1920) laid the foundation of Linear Elastic Fracture Mechanics (LEFM) by stating that quasi-static crack propagation is governed by the balance between the elastic potential energy $P$ and the fracture energy. Written with today's form, it states that the crack propagates when the energy release rate $G$ reaches a critical threshold $G_c$. Francfort and Marigo (1998) generalized this idea in the variational approach to fracture, which states that the structure behavior is governed by the minimization of the structure energy over any displacement field and crack extension (i.e., any sharp discontinuity in the displacement field). Based on works in image segmentation (Mumford and Shah, 1989; Ambrosio and Tortorelli, 1990, 1992), Francfort and Marigo (1998) and Bourdin et al. (2000) proposed to regularize this approach by approximating the sharp cracks with smeared cracks, leading to the definition of variational phase-field models. Thus, instead of explicitly tracking crack

*Corresponding author
Email addresses: flavien.loiseau@ensta.fr (Flavien Loiseau), veronique.lazarus@ensta.fr (Véronique Lazarus)
geometry, cracks are represented by a continuous phase field, which is smeared over a regularization length $\ell$ to prevent sharp discontinuities. Variational phase-field models have been proven to Γ-converge towards the sharp crack model: when the regularization length $\ell$ tends to zero, the regularized fracture energy tends to the non-regularized fracture energy (Bellettini and Coscia, 1994; Bourdin, 1999; Giacomini and Ponsiglione, 2003; Chambolle, 2004). In this study, we employ both the sharp crack model and a variational phase-field model. The latter is utilized with a sufficiently small regularization length $\ell$ to ensure convergence towards the sharp crack model, which serves as our reference.

This study aims to investigate the behavior of a structure with an initial crack under specified loading conditions. The primary objective is to determine the evolution of the propagation threshold as well as the crack path during the failure of the structure. To formalize this idea, we recall the notion of equilibrium path. In the context of LEFM, the *equilibrium path* is the sequence crack length and critical load, corresponding to Griffith crack propagation threshold, satisfying the structural equilibrium. It enables the identification of each critical state during crack propagation that leads to further crack propagation. Thus, the equilibrium path enables the assessment of the structural integrity under a given load at any point during crack propagation. Furthermore, this approach enables the detection of structural-level bifurcations, such as snap-through and snapback phenomena, which can lead to catastrophic failure of the structure. It can also be used to design materials and structures with increased apparent toughness (Roux et al., 2003; Xia et al., 2013; Morano et al., 2020; Lebihain et al., 2021; Triclot et al., 2024).

Determining a point on the equilibrium path using the sharp crack model involves predicting a crack extension and calculating the corresponding load amplitude to satisfy the Griffith criterion $G = G_c$. However, for variational phase-field fracture models, directly imposing the constraint $G = G_c$ is more complex due to the regularization of the crack, which introduces numerical challenges in accurately controlling the energy release rate. Consequently, a practical solution is to employ path-following methods<fnref n="1" /> to determine the equilibrium path. A comprehensive literature already exists to address the problem of determining the equilibrium path of a structure. To solve structural problems with geometric nonlinearities (finite strain, buckling, etc.), arc-length methods have been developed (Wempner, 1971; Riks, 1972, 1979) and successfully applied in finite element simulation (Crisfield, 1981; Ramm, 1981). It incorporates the load factor $\lambda$ as an additional unknown along with a control equation to adapt the loading and avoid potential instabilities that would occur under classic incremental control. Aside their implementation, what mainly differentiates the path-following methods is the choice of the control equation. In particular, the classic arc-length control equations are often limited when a snapback occurs due to localized phenomena (Crisfield, 1982; De Borst, 1987), which is the case of fracture and damage as strain localization occurs in a small region of the domain. Other load control equations have been proposed to overcome this issue and are reviewed by Rastiello et al. (2022). De Borst (1987) introduced an indirect displacement control method based on local displacement measures (e.g., a crack mouth opening displacement), but its success relies on the selection of arbitrary control points. Rastiello et al. (2022) discussed its limitations in complex scenarios, such as multi-cracking, where finding adequate control points becomes highly challenging. To tackle this problem, Geers (1999) formalized this type of approach using a subset of the degrees of freedom in the control equation. Chen and Schreyer (1990) introduced a maximum strain increment control method, independent of the constitutive model and the problem geometry. Alternative approaches include the local arc-length method by May and Duan (1997) and another arc-length method applied only in the process zone by Pohl et al. (2014). The idea behind these methods is to limit the crack growth between two load steps.

While these methods can efficiently track the equilibrium path, they do not ensure that the condition $G = G_c$ is met. Consequently, the solution path may include points that do not correspond to the threshold of crack propagation. In the context of quasi-static crack propagation, such solutions are elastic, in contrast to dissipative solutions where the crack actually propagates. To enforce following a dissipative equilibrium path, various control equations based on dissipative quantities have been suggested, such as controlling the load factor using the maximum elastic predictor (Lorentz and Badel, 2004), or employing energy release methods (Gutiérrez, 2004; Verhoosel et al., 2009). Other approaches, some specific to phase-field models, include using the fracture dissipation (Singh et al., 2016), domain integrals of history variables (Fayezioghani et al., 2019), and internal variable increments (Rastiello et al., 2019). While these methods are often more practical and robust, their implementation often require significant changes to classic phase-field fracture

<footnote n="1">¹The terminology for these methods varies across disciplines. In mathematics and physics, they are commonly designated as numerical continuation methods. Within mechanical engineering, they are predominantly referred to as arc-length methods. In the context of damage/fracture mechanics literature, they are frequently termed path-following methods. For consistency with the latter, we adopt the term path-following methods throughout this work.</footnote>
solvers. The systematic literature review, presented in Appendix A, on the application of path-following methods in phase-field fracture shows that the number of phase-field fracture studies employing path-following is limited (May et al., 2016; Singh et al., 2016; Wu et al., 2022; Zambrano et al., 2023, among others), despite the numerous works on phase-field fracture models. We suppose that this limited adoption is partially due to the (relative) complexity of their numerical implementation, compared to the classic staggered scheme (alternate minimization) employed for phase-field fracture models.

In this study, we investigate path-following techniques for quasi-static crack propagation problems within the LEFM framework. To this aim, we employ both a sharp crack model based on Griffith theory and a variational phase-field fracture model. Our primary objective is to implement and assess a variational phase-field fracture model combined with various path-following control equations. We specifically aim to find a method that enables us to trace the equilibrium path of structures undergoing fracture. Ideally, this method should fulfill several criteria. It should (i) be accurate and match the results of the LEFM sharp crack model, (ii) be problem independent to accommodate complex geometries and loading conditions commonly encountered in engineering applications, (iii) be constitutive model independent so that it is easily applied to other behavior (e.g., multi-physics coupling), and (iv) be minimally intrusive to preserve the simple algorithmic structure of the classic staggered solvers (Bourdin et al., 2000; Miehe et al., 2010; Farrell and Maurini, 2017). To find a method that fullfills those criteria, we compare different existing path-following methods. Furthermore, building upon the work of Chen and Schreyer (1990), we introduce a novel path-following method for the phase-field fracture model that limits the maximum strain increment in the structure, excluding already cracked regions. A path-following method is also developed, for the LEFM sharp crack model, serving as our reference solution. The methods developed for this work are implemented in two open-access codes: one for the LEFM approach (Loiseau and Lazarus, 2026) and one for the phase-field approach (Loiseau and Lazarus, 2025a). The codes are used to assess the convergence and accuracy of phase-field approaches through three numerical benchmarks of increasing complexity. Based on the results, we provide a critical discussion on the use of path-following methods in phase-field fracture, including recommendations for practitioners and potential extension to this work.

This study is organized as follows. In Section 2, we introduce the concept of the equilibrium path in the context of quasi-static crack propagation. First, we discuss this concept within the framework of Linear Elastic Fracture Mechanics (LEFM), corresponding to the sharp crack model. We then present the transition to the variational approach to fracture, leading to the introduction of variational phase-field fracture models (smeared crack model). Section 3 focuses on the numerical solution of the path-following problem. We begin by introducing the numerical solution to the sharp crack model. Subsequently, we present the incorporation of the solution to the path-following control equation in phase-field fracture simulations, focusing on a staggered solver. Then, we present the path-following methods examined in this study. In particular, we propose a new path-following method based on the maximum strain increment on a domain excluding the cracked regions. In Section 4, we present numerical examples to compare the different methods. These examples aim to assess the solution of the smeared crack models by comparing it to the solution of the reference sharp crack model and to evaluate the different path-following equations to determine the most suitable one. Section 5 focuses on a broader examination of the results and potential extensions of this work. Finally, Section 6 summarizes the main conclusions of the study.

## 2. Equilibrium path problem for quasi-static crack propagation

In this section, we introduce path-following in the context of quasi-static crack propagation, as well as its treatment within the frameworks of Griffith theory and the variational phase-field models for fracture.

### 2.1. Crack propagation problem in Griffith theory

We consider an elastic body $\Omega$ containing a crack $\Gamma$, as depicted in Figure 1. The body $\Omega$ is clamped on the boundary $\partial\Omega_0$ and subjected to an external surface force $T$ on the boundary $\partial\Omega_T$. The external force is assumed to be proportional to the load factor $\lambda(t)$; thus, it is recast as $T(t, x) = \lambda(t)\bar{T}(x)$, where $\bar{T}(x)$ carries the spatial variation of the force. We assume that no additional displacements or body forces are imposed to clarify the presentation. The extension to other loads is discussed in Section 5.2. Crack surfaces are also assumed to be free of any external load. The external work associated with the load is given by $W_{\text{ext}} = \int_{\partial\Omega_T} T \cdot u \, dS$. The elastic energy of the body $\Omega$ is denoted $\mathcal{E}$, and the potential energy $\mathcal{P} = \mathcal{E} - W_{\text{ext}}$. In a 2D setting, these quantities are defined by unit thickness. The
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of a crack propagation problem with prescribed force load.</figcaption></figure>

static elastic problem (without crack propagation) consists of finding the displacement field that minimizes the potential energy $P$ of the body $\Omega$. Griffith (1920) proposed to extend this minimization principle to study the rupture of a solid with a crack $\Gamma$ of length $a$. To this aim, he introduced the surface energy $G_c$, also called the fracture energy, as the energy (per surface unit) required to generate new free surfaces<fnref n="2" />. Thus, he states that the equilibrium position of the crack must minimize the energy $P + G_c a$. Deriving this energy with respect to the crack length $a$ leads to the definition of the energy release rate

$$G = -\frac{\partial P}{\partial a} \quad \text{at constant load.} \qquad (1)$$

To ensure the stability of a given crack equilibrium state, the inequality $G \le G_c$ must be imposed, as $G > G_c$ implies the existence of an excess of energy, hence dynamic crack propagation. Moreover, introducing the irreversibility of crack propagation and deriving the total energy with respect to the crack length $a$ gives the following quasi-static crack propagation conditions:

$$(G - G_c)\dot{a} = 0, \quad G \le G_c, \quad \dot{a} \ge 0. \qquad (2)$$

In other words, two distinct cases can occur in this static framework:

$$\begin{cases} G < G_c \Rightarrow \dot{a} = 0, & \text{which means that the crack does not propagate,} \\ G = G_c \Rightarrow \dot{a} \ge 0, & \text{indicating that the crack propagation becomes possible.} \end{cases} \qquad (3)$$

## 2.2. Presentation of the equilibrium path problem

In the case $G > G_c$, there is an excess of energy in the system. This situation arises when $\frac{\partial G}{\partial a} > 0$ at constant load. It implies that, once the threshold is reached, the critical load has to be decreased with increasing crack length to maintain $G = G_c$. In practice, it results in unstable crack growth and potentially catastrophic structural failure. Within the presented analytical framework, this condition violates the quasi-static assumption. Our objective is to determine the sequence of loads that maintains the crack propagation threshold $G = G_c$, thereby identifying overloads that may induce such structural instability. In particular, we aim to characterize the equilibrium path of the structure. In the context of quasi-static crack propagation, the *equilibrium path* denotes the sequence of critical equilibrium states that satisfy the Griffith criterion $G = G_c$ during crack growth. With the assumption introduced in the previous section, an equilibrium state thus comprises three key components: the critical load factor $\lambda_c$ at which propagation occurs, the crack path $\Gamma$, and the corresponding displacement field $u(x)$ of the structure. The solution to this path-following problem requires determining the equilibrium states ($u(t, x)$, $\lambda_c(t)$, $\Gamma(t)$) such that, for each time $t$, the displacement field $u(t, x)$ solves the elastic problem for the fixed crack $\Gamma(t)$, while satisfying the Griffith criterion $G = G_c$ exactly.

To illustrate this concept, we consider a compact tension (CT) specimen containing a short crack along the expected path of the left crack, as depicted in Figure 2a. This case represents a simplified version of the problem that will be discussed in Section 4.2. We consider two loading approaches to apply the point load: either by prescribing the displacement $u = \lambda u_c$ or the force $F = \lambda F_c$, where $u_c$ and $F_c$ denote the critical displacements and forces,

<footnote n="2">²While considered as a reversible energy in Griffith (1920), we consider here that $G_c$ includes all the energy dissipated during propagation.</footnote>
<figure><img src="image_3.png" /><figcaption>(a) CT with crack under punctual loads.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Force-displacement curve.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c) Load factor-crack length curve.</figcaption></figure>

<caption>Figure 2: Illustration of the equilibrium path for a crack bridging problem.</caption>

corresponding respectively, to $G(a_0, u_c) = G_c$ under displacement control and $G(a_0, F_c) = G_c$ under force control. The force-displacement equilibrium path of the structure is illustrated in Figure 2b. This curve represents the force-displacement pairs verifying $G = G_c$ during crack propagation. The load begins with an elastic phase until the first propagation threshold is reached at point A. Then, the left crack starts to propagate, gradually reducing the rigidity of the structure. This first propagation lasts until the left crack merges with the right crack, resulting in a sudden drop in structural rigidity. Then, another elastic phase occurs between the origin and point C, where the structure is reloaded. Thereafter, the crack continues to propagate from the right tip of the newly merged crack. This force-displacement curve shows that the critical force, which verifies $G = G_c$, decreases as soon as the left crack starts to grow at point A, indicating structural softening. Moreover, it also shows that the structural response exhibits a snapback instability starting at point B, as both force and displacement must be decreased to maintain equilibrium until the cracks merge.

This aspect is better illustrated in Figure 2c, which depicts the evolution of the load factor as $G = G_c$ is maintained during crack propagation. Under increasing displacement control, the load factor initially increases between points A and B as the left crack length grows, indicating stable crack propagation. A peak load factor is reached at point B, after which the load factor decreases, signaling the onset of unstable crack growth and corresponding snapback instability. This instability corresponds to the bridging of the two cracks, followed by an elastic equilibrium phase until the second phase of crack propagation begins. This second phase remains stable for most of the crack growth until another unstable propagation occurs at point D, when the crack interacts with the specimen boundary, ultimately leading to complete structural failure. For increasing force control, the interpretation is more straightforward: the first critical load factor (at point A) is the maximum value across all crack lengths. Once the load exceeds this threshold, unstable crack propagation occurs immediately, leading to the total failure of the structure.

These observations have significant implications for numerical simulations of quasi-static crack propagation. Incremental loading schemes are commonly employed, where the load factor is updated iteratively as $\lambda_n = \lambda_{n-1} + \Delta\lambda$, with $\Delta\lambda > 0$. In the presented example, displacement control allows the problem to be solved until the first snapback instability occurs at point B. Following this instability, a crack jump is observed until the next equilibrium state is reached (e.g., when cracks merge at a crack length of 0.5). However, there is no guarantee that this equilibrium point will be attained in numerical simulations. Under force control and in the presence of structural softening, the situation is even more critical. Once the crack propagation threshold is exceeded, the force balance can not be satisfied by any future state of the structure beyond this threshold. In numerical simulations, this leads to solver divergence.

Building on these insights, this paper aims to determine the equilibrium path of the structure as a central focus, from which the different scenarios can be discussed. The following subsections will formalize the path-following problem within the frameworks of LEFM and the variational approach to fracture, providing a rigorous foundation for analyzing the equilibrium path as well as its numerical and physical implications.

## 2.3. Sharp crack model via Griffith theory

Within the quasistatic framework of Griffith theory, an equilibrium state at the instant $t$, where the crack state is $\Gamma(t)$, is achieved by solving the elastic problem under the constraint that the energy release rate equals the critical energy release rate ($G = G_c$). To obtain the complete quasistatic equilibrium path, one must (i) find the crack propagation
threshold, (ii) determine the propagation angle, (iii) increment the crack, and (iv) repeat the procedure for the next equilibrium state. In this study, the Maximum Energy Release Rate (MERR) criterion is used to determine the direction of crack growth (Erdogan and Sih, 1963; Palaniswamy, 1972; Hussain et al., 1973; Nuismer, 1975). We introduce the energy release rate $G^*(s, \varphi)$ accounting for an infinitesimal straight crack increment of length $s$ without a bifurcation angle $\varphi$. To define this quantity, we recall that Amestoy and Leblond (1992) expressed the stress intensity factors $K_p^*(\varphi)$ across a kink of length $s$ as

$$K_p^*(s, \varphi) = F_{pq}(\varphi)K_p + \sqrt{s}TG_p(s) + O(s), \quad (4)$$

where $K_p$, respectively $T$, is the mode $p$ stress intensity factor, respectively $T$-stress, of the crack with the kinked increment. The functions $F_{pq}$ and $G_p$ are universal functions determined by Amestoy and Leblond (1992). The energy release rate $G^*(\varphi)$ across the kink is then given by Irwin's formula (Irwin, 1957)

$$G^*(s, \varphi) = \frac{1}{E'} (K_I^*(s, \varphi) + K_{II}^*(s, \varphi)). \quad (5)$$

In the following, we will take the kink length $s$ sufficiently small to neglect the effect of $T$-stress. According to the MERR criterion, the crack propagates when $G^*(\varphi) = G_c$, with the bifurcation angle $\varphi$ minimizing the potential energy of the body, thereby maximizing the energy release rate. The crack bifurcation angle is thus given by

$$\varphi(t) = \arg\max_{\varphi'} G^*(t, \varphi'). \quad (6)$$

In this framework, the equilibrium path is governed by the crack propagation problem with a controlled crack increment $\Delta a$ between two load steps. For a given crack $\Gamma(t)$, the problem governing the equilibrium state $(u(t, x), \varphi(t), \lambda(t))$ is described by the following system of coupled equations:

$$u(t, x) = \arg\min_{u'} \mathcal{P}(u', \Gamma(t)), \quad \varphi(t) = \arg\max_{\varphi'} G^*(t, \varphi'), \quad G^*(\varphi(t)) = G_c. \quad (7)$$

By incrementally growing the crack by a length $\Delta a$ in the direction $\varphi$ and repeating the procedure, the full equilibrium path can be reconstructed.

## 2.4. Variational approach to fracture

Francfort and Marigo (1998) extended Griffith's theory by formalizing it within a variational framework. Let $\mathcal{U}_\Gamma$ denote the space of admissible displacements, defined as

$$\mathcal{U}_\Gamma = \{u \in H^1(\Omega \setminus \Gamma) \mid \forall x \in \partial\Omega_0, u(x) = 0\}, \quad (8)$$

where $H^1(\Omega \setminus \Gamma)$ represents a classic Sobolev space on the crack body. In this context, the elastic energy and the external work are expressed as

$$\mathcal{E}(u, \Gamma) = \int_{\Omega \setminus \Gamma} \frac{1}{2} \varepsilon(u) : \mathbf{E} : \varepsilon(u) \,dx, \quad \mathcal{W}_{\text{ext}}(u) = \int_{\partial\Omega_{\Gamma}} \mathbf{T} \cdot \mathbf{u} \,dx, \quad (9)$$

where $\varepsilon(u) = \frac{1}{2}(\nabla u + \nabla^T u)$ is the strain tensor and $\mathbf{E}$ is the elasticity tensor. Griffith's theory, as presented in equation (2), can be reformulated as the minimization problem

$$(u, \Gamma) = \arg \min_{\substack{u' \in \mathcal{U}_{\Gamma} \\ \Gamma' \supseteq \Gamma_0}} \mathcal{E}(u', \Gamma') - \mathcal{W}_{\text{ext}}(u') + G_c \mathcal{H}(\Gamma'), \quad (10)$$

where $\mathcal{H}$ denotes the Hausdorff measure of the crack set $\Gamma$. The constraint $\Gamma' \supseteq \Gamma_0 = \Gamma(t_0)$, with $t_0 < t$, ensures crack growth, reflecting the irreversibility of cracking and preventing crack healing. This framework represents a crack as a stress-free surface that can propagate within the domain $\Omega$, inducing a sharp discontinuity in the displacement field and allowing for a displacement jump across the crack. For a given crack length, the next equilibrium state can be determined by solving the problem

$$(u, \Gamma) = \arg \min_{\substack{u' \in \mathcal{U}_{\Gamma} \\ \Gamma' \supseteq \Gamma_0}} \mathcal{E}(u', \Gamma') - \mathcal{W}_{\text{ext}}(u') + G_c \mathcal{H}(\Gamma'), \quad (11)$$

$$\text{subject to } \mathcal{H}(\Gamma) - \mathcal{H}(\Gamma_0) = \Delta a,$$

where the constraint fixes the crack growth rate.
## 2.5. Variational phase-field model for fracture

As the crack corresponds to a discontinuity in the elastic body, the minimization problem in equation (10) is a free-discontinuity problem, which is a well-studied problem (Ambrosio et al., 2000). In the context of image segmentation, free-discontinuity problems have been regularized by replacing sharp discontinuities with a smeared band of non-vanishing width (Mumford and Shah, 1989; Ambrosio and Tortorelli, 1992). Francfort and Marigo (1998) proposed applying the same regularization to the minimization problem in equation (10), and later implemented in the work of Bourdin et al. (2000). Applying this regularization replaces the explicit description of the discontinuity surface $\Gamma$ with a smeared field, called the crack phase field and denoted by $\alpha(\mathbf{x})$. By convention, we consider that the material at coordinate $\mathbf{x}$ is uncracked when $\alpha(\mathbf{x}) = 0$ and fully cracked when $\alpha(\mathbf{x}) = 1$. The Hausdorff measure of the crack set $\Gamma$ is then approximated with

$$ \mathcal{H}_{\ell}(\alpha) = \frac{1}{c_w} \int_{\Omega} \frac{w(\alpha)}{\ell} + \ell \nabla \alpha \cdot \nabla \alpha \, d\mathbf{x} \xrightarrow{\ell \to 0} \mathcal{H}(\Gamma), \quad (12) $$

where $\ell$ is a regularization parameter to which the thickness of the crack phase profile is proportional. The minimization problem in equation (10) governing the state of the domain becomes

$$ (\mathbf{u}, \alpha) = \arg \min_{\substack{\mathbf{u} \in \mathcal{U} \\ \alpha \in \mathcal{A}}} \mathcal{E}(\mathbf{u}, \alpha) - \mathcal{W}_{\text{ext}}(\mathbf{u}) + G_c \mathcal{H}_{\ell}(\alpha) \quad (13) $$

where $\mathcal{U}$ and $\mathcal{A}$ are respectively the admissibility space of the displacement field $\mathbf{u}$ and of the crack phase field $\alpha$. The irreversibility of crack growth is imposed through an irreversibility constraint on the phase field. Given an initial phase field $\alpha_0(\mathbf{x})$, the phase field can stay constant or only grow, leading to the admissibility space

$$ \mathcal{A} = \{ \alpha \in H^1(\Omega) \mid \forall \mathbf{x} \in \Omega, \alpha(\mathbf{x}) \ge \alpha_0(\mathbf{x}) \}. \quad (14) $$

Moreover, the elastic energy and external work functionals become

$$ \begin{aligned} \mathcal{E}(\mathbf{u}, \alpha) &= \int_{\Omega} \frac{1}{2} a(\alpha) \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \, d\mathbf{x}, \\ \mathcal{D}(\alpha) &= \frac{G_c}{c_w} \int_{\Omega} \frac{w(\alpha)}{\ell} + \ell \nabla \alpha \cdot \nabla \alpha \, d\mathbf{x}. \end{aligned} \quad (15) $$

The function $a(\alpha)$ describes how the crack phase affects the initial elastic properties. The function $w(\alpha)$ describes the local dissipation due to the crack phase growth. The scaling constant $c_w$ is derived from the choice of $w$. For the choice of these functions, we retain the AT1 model, proposed by Pham et al. (2011) and shown to be pertinent against experiments by Tanné et al. (2018). In the AT1 model, the functions are

$$ a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha, \quad c_w = 8/3. \quad (16) $$

The energy functional involved in the regularized minimization, in equation (13), i.e. the smeared crack model, has been shown to $\Gamma$-converge towards its sharp crack counterpart, in equation (10), when $\ell \to 0$ (Braides, 1998). In practice, this means that the solution of equation (13) converges to the solution of Griffith's theory, as given in equation (2), when $\ell \to 0$. Therefore, the convergence should also be recovered on the numerical solution in the absence of numerical biases. This key property will be used to evaluate the different path-following methods introduce for phase-field simulations.

For a given crack represented by the phase field $\alpha_0$, the next equilibrium state can be determined by solving the problem

$$ \begin{aligned} (\mathbf{u}, \alpha) = & \arg \min_{\substack{\mathbf{u} \in \mathcal{U} \\ \alpha \in \mathcal{A}}} \mathcal{E}(\mathbf{u}, \alpha) - \mathcal{W}_{\text{ext}}(\mathbf{u}) + G_c \mathcal{H}_{\ell}(\alpha) \\ & \text{subject to } \mathcal{H}_{\ell}(\alpha) - \mathcal{H}_{\ell}(\alpha_0) = \Delta a. \end{aligned} \quad (17) $$

Thanks to the convergence equation (12), the crack length $\mathcal{H}_{\ell}(\alpha)$ tends to its unregularized counterpart $\mathcal{H}(\alpha)$. Thereby, the constraint in the smeared crack path-following problem equation (17) converges to the constraint in the sharp crack path-following problem in equation (11) when the regularization length $\ell$ vanishes.
### 3. Numerical solutions to the path-following problem

This section focuses on the numerical solution of the path-following problem using two distinct approaches: (i) the sharp crack model based on Griffith theory, and (ii) the variational phase-field fracture model. In other words, it presents numerical solution for the problems in equation (11) (sharp crack model) and in equation (17) (phase-field model). For the variational phase-field fracture models, the constraint on crack increment $\Delta a$ will be modified to facilitate its numerical treatment, leading to the introduction of different path-following constraints.

Throughout this section, we will leverage the linearity of the elastic solution with respect to the load factor $\lambda$, as in the work of Lazarus (2003). Consequently, quantities determined for a unit load factor $\lambda = 1$ will be overlined $\bar{\cdot}$. Among others, it is important to note that the following relationships hold

$$ \mathbf{u} = \lambda \bar{\mathbf{u}}, \quad \mathcal{P} = \lambda^2 \bar{\mathcal{P}}, \quad G^*(\varphi) = \lambda^2 \bar{G}^*(\varphi). \qquad (18) $$

This property will be employed throughout this section to adapt the various path-following problems and develop associated numerical methods.

#### 3.1. Sharp crack model solver

For the sharp crack model, following the equilibrium path involves adapting the load factor $\lambda$ to prevent instabilities during crack propagation. We recall that the goal is to ensure that the crack propagates steadily, consistently meeting the Griffith criterion $G^*(\varphi) = G_c$ throughout its progression. To this aim, we exploit the linearity property presented before and assume that the (critical) load factor $\lambda_c$ at the propagation threshold is unknown. In this case, the problem of equation (7) is recast as

$$ \bar{\mathbf{u}}(\mathbf{x}) = \arg \min_{\bar{\mathbf{u}}'} \bar{\mathcal{P}}(\bar{\mathbf{u}}', a), \qquad (19) $$

$$ \varphi = \arg \max_{\varphi'} \bar{G}^*(\varphi'), \qquad (20) $$

$$ \lambda_c^2 \bar{G}^*(\varphi) = G_c. \qquad (21) $$

For a given existing crack, these equations are solved sequentially in the order in which they are presented. First, the finite element method is employed to solve the elastic problem in equation (19). From this solution, the Stress Intensity Factors are determined, here using the I-integral method (Yu and Kuna, 2021) with domain energy integrals (Destuynder et al., 1981). Then, they are injected into the Amestoy and Leblond (1992) formula to determine the energy release rate $\bar{G}^*(\varphi)$ after a kink. The equation (20) is then solved using a gradient descent (a Newton method can also be employed). Finally, the critical load factor $\lambda_c$ directly derives from equation (21). Once the unknowns ($\bar{\mathbf{u}}, \varphi, \lambda_c$) are determined, the crack is incremented with a segment of fixed length $\Delta a$ in the direction $\varphi$. The whole procedure is detailed in Algorithm 1 and its implementation is available in an open repository (Loiseau and Lazarus, 2026) using the FEiCSx library (Baratta et al., 2023).

**Algorithm 1** Numerical solution to the path-following problem for the sharp crack model.

**Require:** Elastic body geometry, initial crack with tip position $\mathbf{x}_0^c$, imposed unit load $\bar{T}$.

1: **for** each load step $n$ **do**

2: Mesh the body with the explicit geometric crack

3: Solve the linear elastic problem (using FEM) to get the displacement for a unit load factor ($\lambda = 1) \rightarrow \bar{\mathbf{u}}_n(\mathbf{x})$

4: Determine the stress intensity factors using the I-integral method (Yu and Kuna, 2021)

5: Define the energy release rate function $G^*(\varphi)$ using Amestoy-Leblond formula (Amestoy and Leblond, 1992)

6: Solve the minimization in equation (20) using a Newton method to get the crack propagation angle $\rightarrow \varphi_n$

7: Compute the critical load factor using the equation (21) $\rightarrow \lambda_{c,n}$

8: (Optional) Compute the scaled displacement field using $\mathbf{u}_n(\mathbf{x}) = \lambda_{c,n} \bar{\mathbf{u}}_n(\mathbf{x}) \rightarrow \mathbf{u}_n(\mathbf{x})$

9: Grow the crack: $\mathbf{x}_n^c = \mathbf{x}_{n-1}^c + \Delta a [\cos(\varphi_n) \mathbf{e}_1 + \sin(\varphi_n) \mathbf{e}_2]$

10: **end for**
### 3.2. Variational phase-field model staggered solver

Let us now focus on the numerical solution of variational phase-field models. The energy functional $\mathcal{E}$ is non-convex, which can cause convergence issues when using a classic monolithic approach (i.e., gradient descent or Newton's method). To avoid this issue, Bourdin et al. (2000) noted that the energy function $\mathcal{E}_\ell$ is polyconvex, i.e., it is : (1) convex in displacement $\mathbf{u}$ at fixed crack phase $\alpha$, and (2) convex in crack-phase $\alpha$ at fixed displacement $\mathbf{u}$. Thus, the authors proposed a staggered solution strategy, often called the alternate minimization algorithm, which iterates between two minimization steps : (1) the minimization with respect to displacement $\mathbf{u}$ at fixed crack phase $\alpha$, and (2) the minimization with respect to crack phase $\alpha$ at fixed displacement $\mathbf{u}$.

The displacement field that minimizes the energy functional (at fixed crack-phase $\alpha$) is a stationary point of the functional $\mathcal{E}_u(\mathbf{u}) = \mathcal{E}(\mathbf{u}, \alpha_{\text{fixed}})$. It means that the (Gâteaux) derivative of the functional with respect to the displacement field $\mathbf{u}$ is null, leading to the variational problem

$$ \text{Find } \mathbf{u} \in \mathcal{U} \text{ such that } \int_{\Omega} a(\alpha)\varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}^*) \mathrm{d}\mathbf{x} + \int_{\partial\Omega_T} \mathbf{T} \cdot \mathbf{u}^* \mathrm{d}\mathbf{x} = 0, \quad \forall \mathbf{u}^* \in \mathcal{U}, \qquad (22) $$

which is a standard linear elasticity problem.

The crack phase field that minimizes the energy functional (at fixed displacement $\mathbf{u}$) is a stationary point of the functional $\mathcal{E}_\alpha(\alpha) = \mathcal{E}(\mathbf{u}_{\text{fixed}}, \alpha)$. It means that the (Gâteaux) derivative of this functional with respect to the solution crack field $\alpha$ must be equal to zero. This leads to the variational problem,

$$ \text{Find } \alpha \in \mathcal{A} \text{ such that } \int_{\Omega} \frac{1}{2} a'(\alpha) \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \alpha^* \mathrm{d}\mathbf{x} + \frac{G_c}{c_w} \int_{\Omega} \frac{w'(\alpha)\alpha^*}{\ell} + 2\ell \nabla\alpha \cdot \nabla\alpha^* \mathrm{d}\mathbf{x} = 0, \quad \forall \alpha^* \in \mathcal{A}, \qquad (23) $$

subject to the irreversibility constraint

$$ \forall \mathbf{x} \in \Omega, \quad \alpha(\mathbf{x}) \ge \alpha_0(\mathbf{x}), \qquad (24) $$

where $\alpha_0$ denotes the initial phase field (before the current load). Thus, this sub-problem is a constrained nonlinear problem.

The alternate minimization consists in alternatively solving those two subproblems until convergence. The procedure for one load step is described in Algorithm 2.

Following classic solvers, the load control strategy typically involves incrementing the load by a fixed increment between successive steps. However, this method may induce dynamic crack propagation if the energy release rate increases with crack length, as it might exceed the critical threshold ($G > G_c$) at constant load. As this control strategy cannot capture the equilibrium path, we will refer to it as the no path-following strategy. The following section introduces path-following methods to avoid dynamic propagation by adjusting the load to maintain $G = G_c$.
**Algorithm 2** Computation of the load step *n* with the classic staggered solution algorithm of variational phase-field model (alternate minimization).

**Require:** Domain mesh, previous state ($u_0(\mathbf{x})$, $\alpha_0(\mathbf{x})$), imposed load $\mathbf{T}$.

1: Initialize the iteration counter $i = 0$

2: **while** not converged **do**

3: Increment the iteration counter $i + 1 \to i$

4: Solve the elastic sub-problem in equation (22) to get the displacement $\to \mathbf{u}^{(i)}(\mathbf{x})$

5: Solve the phase field sub-problem in equation (23) to get the phase field $\to \alpha^{(i)}(\mathbf{x})$

6: Check convergence: $\|\mathbf{u}^{(i)} - \mathbf{u}^{(i-1)}\|_{\infty} < \text{tol}_u$ and $\|\alpha^{(i)} - \alpha^{(i-1)}\|_{\infty} < \text{tol}_{\alpha}$

7: **end while**

**Algorithm 3** Computation of the load step *n* with the adapted staggered solution algorithm of the variational phase-field model (alternate minimization), incorporating the solution of the control equation.

**Require:** Domain mesh, previous state ($u_0(\mathbf{x})$, $\alpha_0(\mathbf{x})$), imposed unit load $\bar{\mathbf{T}}$.

1: Initialize the iteration counter $i = 0$

2: **while** not converged **do**

3: Increment the iteration counter $i + 1 \to i$

4: Solve the elastic sub-problem in equation (26) to get the displacement for a unit load ($\lambda = 1 \to \bar{\mathbf{u}}^{(i)}(\mathbf{x})$)

5: Solve the control equation to get the load factor $\lambda^{(i)}$ satisfying $G = G_c$ (see sections 3.3.1 to 3.3.6)

6: Compute the displacement field $\mathbf{u}^{(i)}(\mathbf{x}) = \lambda^{(i)}\bar{\mathbf{u}}^{(i)}(\mathbf{x})$

7: Solve the phase field sub-problem in equation (23) to get the phase field $\to \alpha^{(i)}(\mathbf{x})$

8: Check convergence: $\|\mathbf{u}^{(i)} - \mathbf{u}^{(i-1)}\|_{\infty} < \text{tol}_u$ and $\|\alpha^{(i)} - \alpha^{(i-1)}\|_{\infty} < \text{tol}_{\alpha}$

9: **end while**

### 3.3. Variational phase-field model path-following solver

Within the phase-field framework, directly prescribing the crack length, as proposed in equation (17), is feasible but requires significant modifications to the solution method (Singh et al., 2016). Instead of directly constraining the crack growth, an indirect control method can be employed using a quantity $\Delta\tau$ related to crack growth, which we refer to as the control quantity in this study. The control equation can be expressed as

$$ \Delta\tau(\lambda, \mathbf{u}_0, \alpha_0, \mathbf{u}, \alpha, \dots) = \Delta\tau_{\text{imp}}, \quad (25) $$

where $\Delta\tau_{\text{imp}}$ is the prescribed step size. The control equation (25) adjusts the load factor $\lambda$ such that the control quantity $\Delta\tau$ increases by the step size $\Delta\tau_{\text{imp}}$ between two consecutive load steps. Then, the goal is to find a control quantity $\Delta\tau$ that relates to the crack growth, allowing for indirect control of the crack growth rate.

An overview of existing work on path-following methods applied to phase-field fracture is provided in Appendix A. It evidences that among the numerous studies on phase-field fracture, only a few employ path-following methods. We suppose that this limited adoption is partially due to the complexity of their numerical implementation. To maintain the simplicity of the alternate minimization, we rely on control equations based on displacement, as they enable factoring the load factor from the elastic sub-problem using the linearity property. Hence, it decouples the control equation from the linear displacement sub-problem. Thanks to the staggered approach, the nonlinear crack phase sub-problem is already decoupled from the load factor $\lambda$. Those decoupling enables to treat the control equation and each sub-problem separately.

We can now focus on the integration of the solution of the control equation in the alternate minimization presented in Algorithm 2. Replacing $\mathbf{u}(\mathbf{x}) = \lambda\bar{\mathbf{u}}(\mathbf{x})$ and cancelling the load factor, the displacement sub-problem of equation (22) becomes

$$ \text{Find } \bar{\mathbf{u}} \in \mathcal{U} \text{ such that } \int_{\Omega} a(\alpha) \varepsilon(\bar{\mathbf{u}}) : \mathbf{E} : \varepsilon(\bar{\mathbf{u}}^*) d\mathbf{x} + \int_{\partial\Omega_{\tau}} \bar{\mathbf{T}} \cdot \bar{\mathbf{u}}^* d\mathbf{x} = 0, \quad \forall \bar{\mathbf{u}}^* \in \bar{\mathcal{U}}. \quad (26) $$

Note that this formulation is now independent of the load factor $\lambda$, as expected. Based on this reformulation, the staggered scheme presented in Algorithm 2 can be adapted to incorporate the path-following control equation. The
<table><caption>Table 1: Synthetic summary of the path-following methods. For the detailed notations, refer to the subsection associated with the method.</caption><thead><tr><th>Name</th><th>Control Equation</th><th>Unit</th><th>Interpretation of control quantity</th></tr></thead><tbody><tr><td>CNDI</td><td>√<span style="text-decoration: overline;">∑<sub>x<sub>c</sub>∈C</sub> (u(x<sub>c</sub>) - u<sub>0</sub>(x<sub>c</sub>))<sup>2</sup></span> - Δτ<sub>imp</sub> = 0</td><td>m</td><td>Nodal displacement increment at x<sub>c</sub></td></tr><tr><td>CLAL</td><td>max<sub>x∈Ω</sub> (<span style="text-decoration: overline;">(u(x) - u<sub>0</sub>(x))</span> · <span style="text-decoration: overline;">u<sub>0</sub>(x)</span>)/||u<sub>0</sub>(x)||) - Δτ<sub>imp</sub> = 0</td><td>m</td><td>Max displacement increment in Ω</td></tr><tr><td>CMSI</td><td>max<sub>x∈Ω</sub> (<span style="text-decoration: overline;">Δε(x) : <span style="text-decoration: overline;">ε(u<sub>0</sub>(x))</span></span></span>)/||ε(u<sub>0</sub>(x))||) - Δτ<sub>imp</sub> = 0</td><td>-</td><td>Max strain increment in Ω</td></tr><tr><td>CCDV</td><td>Y(u(x), α(x)) - Y(u<sub>0</sub>(x), α<sub>0</sub>(x)) - Δτ<sub>imp</sub> = 0</td><td>J</td><td>Dissipation per load step</td></tr><tr><td>CMTCDV</td><td>Y<sub>max</sub>(u(x), α(x)) - Y<sub>max,0</sub> - Δτ<sub>imp</sub> = 0</td><td>J</td><td>Dissipation per load step</td></tr><tr><td>CMSIOC</td><td>max<sub>x∈Ω</sub> (<span style="text-decoration: overline;">(1 - α<sub>0</sub>(x))Δε(x) : <span style="text-decoration: overline;">ε(u<sub>0</sub>(x))</span></span></span>)/||ε(u<sub>0</sub>(x))||) - Δτ<sub>imp</sub> = 0</td><td>-</td><td>Max strain increment in Ω \ Γ</td></tr></tbody></table>

approach involves splitting the solution of the elastic sub-problem into two distinct steps: first, solving the elastic sub-problem for a unit load factor $\lambda = 1$ to determine $\bar{u}$; and second, solving the control equation to determine the critical load factor $\lambda_c$. The adapted procedure is summarized in Algorithm 3. Since the phase-field sub-problem remains unchanged, minimal adaptations are required to implement the path-following control equation in an existing staggered solver.

The remainder of this section focuses on presenting various control equations and their respective solutions, intending to identify a control quantity related to the cracking process that can be used to control it. In particular, the novel Control by Maximum Strain Increment Outside the Crack is introduced. The solutions are presented in a way that seamlessly integrates into Algorithm 3. Each subsection is organized as follows: (1) presentation of the control equation and its physical interpretation, and (2) determination of its solution. For conciseness, the different control equations along with a brief interpretation are gathered in Table 1. The methods presented here have been implemented using the FEniCSx library (Baratta et al., 2023), and is available in an open repository (Loiseau and Lazarus, 2025a).

In the following subsections, to solve the control equations, we assume that the solution from the previous load step is known, including the displacement field and the crack phase. Additionally, we consider being at iteration $i$ of the staggered algorithm. Some control equations will also use the crack phase field $\alpha^{(i-1)}$ from the current load step at the previous iteration of the staggered scheme. When the iteration number is not specified for quantities at the current load step, it is implied that we are referring to the quantity at the current iteration $i$.

### 3.3.1. Control by Nodal Displacement Increment (CNDI)

The CNDI, proposed by De Borst (1987), corresponds to a classic control strategy in experimental studies. The core idea is to indirectly control the load through a combination of displacements at specific control points, denoted as $C = \{x_c^1, x_c^2, \dots\}$. These control points are arbitrarily selected prior to running the simulation. For instance, this approach can replicate the control by Crack Mouth Opening Displacement (CMOD) in a three-point bending test. The associated control equation can be formalized as

$$ \sqrt{\sum_{x_c \in C} (u(x_c) - u_0(x_c))^2} - \Delta\tau_{\text{imp}} = 0. \qquad (\text{CNDI}) $$
The linearity property $u = \lambda \bar{u}$ enables the solution of the control equation. The equation is recast as a second-order polynomial in the load factor $\lambda$

$$a\lambda^2 + b\lambda + c = 0, \quad \text{where} \quad \begin{cases} a = \sum_{x_c \in C} \bar{u}(x_c) \cdot \bar{u}(x_c), \\ b = 2 \sum_{x_c \in C} \bar{u}(x_c) \cdot u_0(x_c), \\ c = -\Delta \tau_{\text{imp}}^2 + \sum_{x_c \in C} u_0(x_c) \cdot u_0(x_c). \end{cases} \qquad (27)$$

Note that there are two real solutions to this problem. We choose the solution corresponding to the largest (and positive) load factor

$$\lambda = \frac{-b + \sqrt{\Delta}}{2a}, \quad \Delta = b^2 - 4ac. \qquad (28)$$

The solving strategy involves computing the polynomial coefficients $a, b$, and $c$ as well as the discriminant $\Delta$. Then, the formula in equation (28) directly provides the solution for the load factor $\lambda$.

### 3.3.2. Control by Local Arc Length (CLAL)

Selecting the control node can be complex and, in some cases, even impossible. This challenge has inspired other authors to develop methods that automatically choose the appropriate degree(s) of freedom (May and Duan, 1997; Geers, 1999; Pohl et al., 2014, among others). These techniques are generally referred to as local arc-length methods. The core idea is to define a constraint that applies to the degrees of freedom most affected by the dissipation process. This can be done either by defining a local control quantity or by applying the classic arc-length constraint to degrees of freedom that are automatically selected in the process zone. Here, we focus on control by maximum displacement increment, which aims to limit the largest nodal displacement increment between two load steps within the structure. The control equation for this method is given by

$$\max_{x \in \Omega} \left( (u(x) - u_0(x)) \cdot \frac{u_0(x)}{\|u_0(x)\|} \right) - \Delta \tau_{\text{imp}} = 0. \qquad (\text{CLAL})$$

Note that the step size $\Delta \tau_{\text{imp}}$ represents a displacement, making it straightforward to interpret and select.

The control equation can be solved by substituting the displacement field with its expression $u = \lambda \bar{u}$ using the linearity property. This substitution yields

$$\max_{x \in \Omega} (\lambda c_1(x) + c_0(x)) - \Delta \tau_{\text{imp}} = 0, \quad \text{with} \quad \begin{cases} c_0(x) = -u_0(x) \cdot \frac{u_0(x)}{\|u_0(x)\|}, \\ c_1(x) = \bar{u}(x) \cdot \frac{u_0(x)}{\|u_0(x)\|}. \end{cases} \qquad (29)$$

To solve this type of equation, we employ the nested interval algorithm proposed by Lorentz and Badel (2004). For conciseness, the algorithm is detailed in Appendix B.

### 3.3.3. Control by Maximum Strain Increment (CMSI)

Chen and Schreyer (1990) proposed the Control by Maximum Strain Increment (CMSI), recently revisited by Rastiello et al. (2022). This method can also be considered as a local arc-length method, as it aims to limit the maximum strain increment within the domain. The control equation is given by

$$\max_{x \in \Omega} \left( \Delta \varepsilon(x) : \frac{\varepsilon(u_0(x))}{\|\varepsilon(u_0(x))\|} \right) - \Delta \tau_{\text{imp}} = 0, \qquad (\text{CMSI})$$

where the strain increment is $\Delta \varepsilon(x) = \varepsilon(\lambda \bar{u}(x) - u_0(x))$. This method is similar to the CLAL method presented before, with the key difference being that the critical node becomes a critical element. Note that, since the numerical strain
measure depends on the element size, the step size parameter $\Delta\tau_{imp}$ in the CMSI must be adjusted according to the mesh size $h$.

By substituting the strain increment with its expression, we obtain

$$ \max_{\mathbf{x} \in \Omega} (\lambda c_1(\mathbf{x}) + c_0(\mathbf{x})) - \Delta\tau_{\text{imp}} = 0, \text{ with} \quad \begin{cases} c_0(\mathbf{x}) = -\varepsilon(\mathbf{u}_0(\mathbf{x})) : \frac{\varepsilon(\mathbf{u}_0(\mathbf{x}))}{\|\varepsilon(\mathbf{u}_0(\mathbf{x}))\|}, \\ c_1(\mathbf{x}) = \varepsilon(\bar{\mathbf{u}}(\mathbf{x})) : \frac{\varepsilon(\mathbf{u}_0(\mathbf{x}))}{\|\varepsilon(\mathbf{u}_0(\mathbf{x}))\|} \end{cases} \quad (30) $$

To solve this equation, we employ the nested interval algorithm proposed by Lorentz and Badel (2004) and detailed in Appendix B.

### 3.3.4. Control by Crack Driving Variable (CCDV)

Zambrano et al. (2023) proposed the Control by Crack Driving Variable (CCDV) method for phase-field fracture simulations. As we do not use the history variable, we propose an alternative formulation of the control equation of Zambrano et al. (2023). The crack driving variable is obtained by taking the Gâteaux derivative of the energy functional with respect to the crack phase-field. By evaluating this derivative in the direction $\alpha^* = \alpha$, we obtain

$$ \frac{G_c}{c_w} \int_{\Omega} \frac{w'(\alpha)\alpha}{\ell} + 2\ell \nabla\alpha \cdot \nabla\alpha \, d\mathbf{x} = Y(\mathbf{u}, \alpha), \quad \text{where } Y(\mathbf{u}, \alpha) = - \int_{\Omega} \frac{1}{2} a'(\alpha)\alpha \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \, d\mathbf{x}. \quad (31) $$

The right-hand side represents the quantity driving the crack phase-field, analogous to the energy release rate. Thus, controlling this quantity allows direct control over the crack evolution and the dissipation process. Accordingly, we introduce the control equation

$$ Y(\mathbf{u}(\mathbf{x}), \alpha(\mathbf{x})) - Y(\mathbf{u}_0(\mathbf{x}), \alpha_0(\mathbf{x})) - \Delta\tau_{\text{imp}} = 0. \quad (\text{CCDV}) $$

In practice, the explicit expression of the crack phase-field $\alpha$ as a function of the load factor $\lambda$ is unknown. However, we adhere to the staggered solution approach and use the crack phase-field $\alpha^{(i-1)}$ from the previous alternate minimization iteration.

The solution of the control equation of the CCDV is obtained by substituting the displacement with its expression. It gives

$$ \lambda^2 Y(\bar{\mathbf{u}}, \alpha) - Y(\mathbf{u}_0, \alpha_0) - \Delta\tau_{\text{imp}} = 0 \quad \Rightarrow \quad \lambda = \pm \sqrt{\frac{Y(\mathbf{u}_0, \alpha_0) + \Delta\tau_{\text{imp}}}{Y(\bar{\mathbf{u}}, \alpha^{(i-1)})}}. \quad (32) $$

We select the solution corresponding to the positive load factor.

### 3.3.5. Control by Maximum over Time of Crack Driving Variable (CMTCDV)

The Control by Crack Driving Variable (CCDV) method may fail in scenarios where the energy release rate decreases during crack propagation. To align more closely with the core method proposed by Zambrano et al. (2023), we modify the CCDV to control the maximum over time $Y_{\text{max}}$ of the crack driving variable $Y$, defined as

$$ Y_{\text{max}} = - \int_{\Omega} \max_{t} \left( \frac{1}{2} a'(\alpha) \alpha \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \right) d\mathbf{x}. \quad (33) $$

This approach corresponds to using the history variable introduced by Miehe et al. (2010) in the control equation. However, note that the irreversibility constraint is still imposed through a variational inequality.

The control equation for the Control by Maximum of Crack Driving Variable (CMTCDV) is given by

$$ Y_{\text{max}}(\mathbf{u}(\mathbf{x}), \alpha(\mathbf{x})) - Y_{\text{max},0} - \Delta\tau_{\text{imp}} = 0. \quad (\text{CMTCDV}) $$

Due to the maximum over time, this equation is nonlinear in $\lambda$. To solve this nonlinear equation, a scalar nonlinear solver can be employed. In practice, we use the `scipy.optimize. root_scalar` function from the SciPy library (Virtanen et al., 2020) to solve this equation.
### 3.3.6. New proposition: Control by Maximum Strain Increment Outside of the Crack (CMSIOC)

We introduce the Control by Maximum Strain Increment Outside of the Crack (CMSIOC), a novel control method inspired by the CMSI approach proposed by Chen and Schreyer (1990), presented in Section 3.3.3. In this method, we replace the maximization over the whole domain $\Omega$ with a maximization over the uncracked domain $\Omega \setminus \Gamma$. This modification excludes the strain increments that may occur in the cracked region ($\alpha = 1$), which are proportional to crack opening. By doing so, it ensures that the controlled maximum strain increment is localized in the cracking zone ahead of the crack tip, where strain localization occurs due to the growth of the crack. This control equation thus enforces a physically meaningful limit on crack propagation. Rather than geometrically restricting the maximization domain, we employ a degradation function to weight the maximum strain increment, leading to

$$ \max_{x \in \Omega} \left( (1 - \alpha_0(\mathbf{x})) \Delta\varepsilon(\mathbf{x}) : \frac{\varepsilon(\mathbf{u}_0(\mathbf{x}))}{\|\varepsilon(\mathbf{u}_0(\mathbf{x}))\|} \right) - \Delta\tau_{\text{imp}} = 0. \qquad (\text{CMSIOC}) $$

The control equation employs the phase field $\alpha_0$ from the previous load step as a reference. This ensures that the region of active crack propagation (occurring between the current and previous load steps) is included in the maximization domain. A similar reasoning is used in the $\kappa - \omega$ control equation proposed by Fayezioghani et al. (2019) with the current damage. As for the CMSI, the numerical strain measure depends on the size of an element, the step size parameter $\Delta\tau_{\text{imp}}$ in the CMSIOC must be adapted according to the mesh size $h$.

Similarly to the CMSI, the strain increment can be substituted by its expression, leading to

$$ \max_{x \in \Omega} (\lambda c_1(\mathbf{x}) + c_0(\mathbf{x})) - \Delta\tau_{\text{imp}} = 0, \quad \text{with} \quad \left\{ \begin{array}{l} c_0(\mathbf{x}) = -(1 - \alpha_0(\mathbf{x}))\varepsilon(\mathbf{u}_0(\mathbf{x})) : \frac{\varepsilon(\mathbf{u}_0(\mathbf{x}))}{\|\varepsilon(\mathbf{u}_0(\mathbf{x}))\|}, \\ c_1(\mathbf{x}) = (1 - \alpha_0(\mathbf{x}))\varepsilon(\bar{u}(\mathbf{x})) : \frac{\varepsilon(\mathbf{u}_0(\mathbf{x}))}{\|\varepsilon(\mathbf{u}_0(\mathbf{x}))\|} \end{array} \right. \qquad (34) $$

To solve this equation, we employ the nested interval algorithm proposed by Lorentz and Badel (2004) and detailed in Appendix B.

## 3.4. Numerical considerations for phase-field fracture simulations

This subsection outlines key numerical adjustments required for accurate phase-field fracture simulations, including corrections for the energy release rate, initialization of initial cracks, and measurement of crack length. That is:

* **Correction of $G_c$.** To ensure the correct effective energy release rate, we use the relation proposed by Bourdin et al. (2008),

$$ G_c = G_c^{\text{PF}} \left( 1 + \frac{h}{c_w \ell} \right), \qquad (35) $$

where $h$ denotes the mesh size along the crack path. This relation accounts for the one-element-wide plateau in the crack phase profile induced by the (continuous) finite element discretization.

* **Representation of initial cracks.** According to a previous study (Loiseau and Lazarus, 2025b), we introduce initial cracks by removing a one-element-wide band of elements along the initial crack path and imposing a zero Dirichlet boundary condition, $\alpha(\mathbf{x}) = 0$, on the crack faces. By comparing it to the sharp crack model, we have demonstrated that this technique circumvents the mesh bias introduced by other crack initialization techniques, enabling the accurate estimation of the crack propagation threshold.

* **Measurement of crack length.** We evaluate the crack length in phase-field simulations using the regularized Hausdorff measure defined in equation (12). However, since the initial crack region is partially excluded from the mesh, a correction is necessary. Specifically, the missing region $\Omega_{\text{missing}} = [0, a_0] \times [0, h]$ must be accounted for in the crack length measurement $H_c$. The corrected crack length $a_{\text{PF}}$ is thus

$$ a_{\text{PF}} = \underbrace{H_c(\alpha)}_{\text{in mesh equation (12)}} + \underbrace{H_{\text{missing}}}_{\text{outside mesh}} = H_c(\alpha) + \frac{1}{c_w} \int_{\Omega_{\text{missing}}} \underbrace{\frac{w(\alpha)}{\ell}}_{w(\alpha(\mathbf{x}))=1} + \ell \nabla\alpha \cdot \nabla\alpha \, d\mathbf{x} = H_c(\alpha) + \frac{ha_0}{c_w \ell}. \qquad (36) $$

This expression applies to both the AT1 and AT2 phase-field models.
<table><caption>Table 2: Step size employed for the first numerical example (SENT test).</caption><thead><tr><th></th><th>No PF</th><th>CNDI</th><th>CLAL</th><th>CMSI</th><th>CCDV</th><th>CMTCDV</th><th>CMSIOC</th></tr></thead><tbody><tr><td>Step size Δτ<sub>imp</sub></td><td>6 × 10<sup>-8</sup> m</td><td>2 × 10<sup>-8</sup> m</td><td>2 × 10<sup>-7</sup> m</td><td>0.05</td><td>6.75 × 10<sup>-3</sup> J</td><td>6.75 × 10<sup>-3</sup> J</td><td>0.01</td></tr></tbody></table>

## 4. Selected numerical examples

This section aims to demonstrate the ability of phase-field fracture models to replicate the equilibrium path predicted by the sharp crack model in crack propagation problems. In particular, the capabilities and performances of the various path-following equations are evaluated. To achieve this, the path-following methods are applied to three numerical examples of increasing complexity, illustrated in Figure 3. For each example, the equilibrium path obtained using the phase-field model is compared to that of the sharp crack model through three means:

1. the force vs. displacement (F-u) curve, which characterizes the overall structural response,

2. the load factor vs. crack length (λ-a) curve, illustrating the evolution of the critical load during crack propagation,

3. snapshots of the deformed domain for different crack lengths, providing a qualitative assessment of the applied load and crack opening.

For the three examples, the meshes for phase-field fracture simulations are partitioned into two regions (illustrated in Figure 4): (1) a region with coarse elements of size $h_{max}$, and (2) a region where the crack is expected with refined elements of size $h$. The choice of mesh size $h$ is based on the regularization length $\ell$ (i.e., $h = \ell/N$, where $N$ is specified in each example). For the last two meshes, a transition region, denoted $h_{max} \to h$, where the element size evolves linearly between both values, is also introduced.

### 4.1. Example 1: Single Edge Notched Tensile (SENT) test

The first numerical example involves a square elastic body with an initial crack subjected to a tensile load, as illustrated in Figure 3a. This configuration is known as the Single Edge Notched Tension (SENT) test. It follows the setup described by Singh et al. (2016). This simple academic example is commonly studied in the phase-field fracture community (Bourdin et al., 2008; Miehe et al., 2010; Vignollet et al., 2014; Singh et al., 2016, among many others). It exhibits a snapback instability leading to dynamic crack propagation above the threshold under both increasing displacement and force control. Here, it serves two purposes: (1) to evaluate the capacity of the path-following constraints to capture the equilibrium path, and (2) to investigate the mesh convergence in $h$ of the phase-field fracture simulations toward LEFM-based sharp crack predictions for a sufficiently small regularization length $\ell$.

#### 4.1.1. Implementation details

In this study, the regularization length $\ell = 3/200 \times L = 1.5 \times 10^{-5}$ m is kept constant. According to the LEFM prediction, the crack path is expected to remain straight, allowing us to employ a structured mesh aligned with the crack path to minimize numerical artifacts in phase field simulations. The bulk mesh size is $h_{max} = 2\ell$, and mesh size in the crack region is varied $h = \ell/N$ with $N \in \{2, 4, 6, 8\}$. The step size parameter is reported in Table 2 for each path-following method. For the CNDI, the two control nodes are placed $3 \times 10^{-4}$ m above and below the initial crack tip.

#### 4.1.2. Comparison of the control equations

Figure 5 shows the F-u and λ-a equilibrium paths for the SENT specimen, using the LEFM sharp crack model as the reference and comparing it with the phase-field fracture model using the different path-following methods. In both plots, the phase-field results obtained with each path-following constraint perfectly overlap until the constraint fails to track the equilibrium path. The point of mismatch for each path-following constraint is marked, and the subsequent portion of the curve is depicted as a dotted line. First, the equilibrium path is not captured at all when no path-following constraint is applied, as anticipated. It results in a significant crack jump during a single load step, violating the quasi-static assumption, which might lead to physically unrealistic results. In contrast, all the other methods capture the equilibrium path predicted by the LEFM sharp crack model up to a certain point. Hence, as far as the equilibrium
<figure><img src="image_6.png" /></figure>

Plane strain

$$E = 230.77 \text{ MPa} \quad \nu = 0.43$$

$$G_c = 2700 \text{ J} \cdot \text{m}^{-2} \quad \ell = 15 \times 10^{-6} \text{ m}$$

(a) Single Edge Notched Tensile (SENT) specimen under tension.

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

Plane strain

$$E = 230 \text{ MPa} \quad \nu = 0.3$$

$$G_c = 2700 \text{ J} \cdot \text{m}^{-2} \quad \ell \text{ varied}$$

(b) Compact Tension (CT) specimen with a secondary crack (force load).

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

Plane stress

$$E = 3 \text{ GPa} \quad \nu = 0.36$$

$$G_c = 700 \text{ J} \cdot \text{m}^{-2} \quad \ell \text{ varied}$$

(c) Compact Tension (CT) specimen with an eccentric hole inducing mixed-mode I+II sollicitation (force load).

<caption>Figure 3: Illustration of the numerical examples with increasing complexity.</caption>
<figure><img src="image_9.png" /><figcaption>Figure 4: Element size distribution in the meshes.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 5: Equilibrium path of the SENT specimen determined with the sharp crack model (LEFM) and the phase-field model with different path-following constraint. The dashed parts of the curves correspond to a crack jump and the extremities of the jump are highlighted with markers (for CMTCDV, CMSI, and CMSIOC, the marker denotes the end of the simulation).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 6: Equilibrium path of the SENT specimen determined with the sharp crack model (LEFM) and the phase-field model (CMSI) with different mesh sizes.</figcaption></figure>

path is captured, the choice of the path-following method does not matter. Three methods, the CMSI, the CMSIOC, and the CMTCDV, successfully capture the full equilibrium path. The other methods capture a more or less significant portion of the equilibrium path but fail when the crack becomes large and begins to interact with the boundary. An abrupt drop in the critical load factor occurs when only a small ligament remains, corresponding to a strong unloading of the structure.

### 4.1.3. Convergence with mesh size

The previous results were obtained using solutions converged in mesh size *h*. Here, we analyze the convergence behavior of the phase-field solution as the mesh size *h* is varied. For this analysis, we focus on the CMSI, as it fully captures the snapback instability and is representative of a working path-following method. The bulk mesh size is fixed at $h_{max} = 2\ell$, while the mesh size *h* along the expected crack path is varied. Figure 6 illustrates the outcomes of this mesh convergence study. Due to the use of a structured mesh that perfectly aligns with the expected crack path, convergence is achieved rapidly. Furthermore, the correction formula provided in equation (35) removes the influence of the mesh size on the effective energy release rate. It ensures that the post-peak behavior in the *F*–*u* curve aligns accurately with the reference solution, regardless of the mesh size. Here, we attribute the minor remaining errors to the incomplete convergence in bulk mesh size $h_{max}$, and in the regularization length $\ell$.

Moreover, Figure 7 compares the deformed configurations obtained from the converged phase-field model (colored by the phase field) with those from the sharp crack model. The snapshots of both models are synchronized based on the crack length. Although the critical imposed displacement is slightly higher in the LEFM solution (see the bottom and top deformed boundaries), the crack opening profiles align remarkably well, except near the crack tip, where differences are expected due to the regularization effect. This result illustrates the close correspondence between the LEFM sharp crack model and the variational phase-field fracture model, when the latter is converged.

### 4.2. Example 2: *Compact Tension (CT) specimen under force load with a secondary crack*

The second numerical example, described in Figure 3b, is a pre-cracked CT specimen with a secondary crack along the expected crack path. This secondary crack is introduced to induce a snapback, which occurs when the left crack bridges with the other one. Moreover, this example features force boundary conditions on the top hole. As this example displays a structural softening, it can not be tackled without path-following methods. Indeed, after reaching the propagation threshold, no successive equilibrium state can balance the first critical force.

While this problem might seem to be similar to the *en passant* problem (Melin, 1983), only one crack tip is loaded up to its propagation threshold. Thus, the crack tips do not repel each other, and the crack follows a straight path.
<figure><img src="image_12.png" /><figcaption>Figure 7: Two snapshots of the phase-field solution in the SENT specimen for various crack lengths (extracted from the CMSI simulations). The mesh is deformed with 5 times the displacement field and the color represents the crack phase field (dark is sane, and light is cracked). The red outline represents the deformed boundary of the LEFM solution. Note that the two superposed configurations have been translated vertically to center the crack.</figcaption></figure>

<table><caption>Table 3: Step size employed for the second numerical example (CT with crack).</caption><thead><tr><th></th><th>No PF</th><th>CNDI</th><th>CLAL</th><th>CMSI</th><th>CCDV</th><th>CMTCDV</th><th>CMSIOC</th></tr></thead><tbody><tr><td>Step size Δτ<sub>imp</sub></td><td>10<sup>7</sup> m</td><td>10<sup>-6</sup> m</td><td>5 × 10<sup>-7</sup> m</td><td>0.25</td><td>10<sup>-2</sup> J</td><td>10<sup>-2</sup> J</td><td>10<sup>-2</sup></td></tr></tbody></table>

Note that we also observe a small interpenetration of the crack lips on the right side of the right crack. While we could account for this contact, it would significantly complicate the problem. Since the objective here is to assess the path-following methods, we neglect the contact between the crack faces. This assumption is applied to both sharp crack and phase-field models, thereby ensuring that it does not influence the comparison of the methods.

For this example, we proceed as follows. In Section 4.2.1, we present some implementation details specific to this problem. In Section 4.2.2, we compare the different path-following methods to determine which methods can capture the equilibrium path. This analysis is conducted on a coarse problem, i.e., with a coarse mesh and a large regularization length, to reduce computational cost. Additionally, this first step is done without any comparison to the reference sharp crack model. As the equilibrium path is unaffected by the path-following method, the goal is to select a single method to use for the remainder of this example. In Section 4.2.3, we perform a mesh convergence study on both the LEFM sharp crack model and the phase-field model (with the selected path-following method). In this study, the bulk mesh size $h_{max}$ is varied. The other numerical parameters, i.e. $ℓ$ and $h$, are kept constant. Finally, in Section 4.2.4, we carry out a numerical $\Gamma$-convergence study. To this end, we perform phase-field simulations with different regularization lengths $ℓ$ (at a constant $ℓ/h$) and compare them to the reference sharp crack model.

### 4.2.1. Implementation details

The step size parameter is reported in Table 3 for each path-following method. In phase-field simulations, an unstructured mesh is employed to avoid any mesh-induced bias (Negri, 2003, 2007) that could influence the predicted crack path. The mesh sizes will be specified in each study. Moreover, to facilitate comparison, the crack length $a$ represents the sum of the lengths of the existing cracks. Hence, the initial crack length is $a_0 = L/2 + L/8$.
<figure><img src="image_13.png" /><figcaption>Figure 8: Equilibrium path of the CT specimen with a secondary crack determined using a coarse phase-field model with different path-following constraint. The dashed parts of the curves correspond to a crack jump and the extremities of the jump are highlighted with markers (for CMTCDV, CMSI, and CMSIOC, the marker denotes the end of the simulation). For readability, the x-axis of the force-displacement curve has been limited to the range $u \in [0, 4 \times 10^{-5}]$.</figcaption></figure>

### 4.2.2. Comparison of the path-following methods

Figure 8 illustrates the equilibrium path obtained from the phase field model using various path-following equations. Those simulations were performed with a large regularization length ($\ell = L/64$) and a coarse mesh ($h_{max} = L/64$, $h = \ell/10$) to assess, with a minimal computational cost, which methods correctly capture the equilibrium path. Figure 8b shows that the load factor reaches a peak as soon as the crack starts to propagate. It means that for any longer crack length, the structure will be unable to support the load, and the crack will propagate, leading to unstable crack propagation. As there is a force load, it also means that no cohesive forces will be able to balance this load, leading to the failure of the classic solver with no path-following. Moreover, this example also confirms the previous observation of Section 4.1.2 that while the equilibrium is successfully captured, the different methods yield results that perfectly superimpose. Hence, the predicted equilibrium path is independent of the path-following method. In this example, the three methods that fully capture the snapback are the CMTCDV, the CMSI, and the CMSIOC. Note that the methods that failed this example also failed the previous (SENT) example.

For the following studies, we want to focus on using a single path-following method. A criterion is needed to select a single method among the CMTCDV, the CMSI, and the CMSIOC. As a selection criterion, we propose to aim for a constant crack growth rate. It ensures that each load step yields the same crack increment, while skipping the elastic phases to limit unnecessary calculations.

To evaluate which methods best fit this criterion, Figure 9 represents the evolution of crack length as a function of the (dimensionless) load step for the CMTCDV, the CMSI, and the CMSIOC. In this figure, horizontal segments indicate load increases without crack growth (elastic phase). Vertical jumps signify unstable crack propagation between two successive load steps. Thus, horizontal and vertical lines are to be avoided. The ideal method, according to our selection criterion, provides a linear increase in crack length with load steps.

For the CMSI, we observe that the slope is first linear, and then it decreases when the crack length increases, indicating a strong decrease in crack growth rate. Indeed, when the crack length is small, the maximum strain increment is in the cracking region ahead of the crack tip. Its control leads to a nearly linear crack growth with load steps. However, when the crack length is sufficiently large, the maximum strain increment occurs in the cracked region ($\alpha = 1$), and the increment is proportional to a crack opening displacement (COD). Hence, the CMSI acts as a limiter of COD. However, as the crack grows, a constant increases of COD leads to a increasingly small crack increments due to the growing lever arm. It leads to the strong decrease in crack growth rate, thus requiring a large number of load steps to reach full failure.
<figure><img src="image_14.png" /><figcaption>Figure 9: Evolution of the crack length with the load step during a simulation of the CT specimen with a secondary crack. The dashed line represents a constant crack growth rate.</figcaption></figure>

In contrast, the CMTCDV and CMSIOC methods closely approach the proposed ideal, demonstrating good performance in controlling crack growth. In particular, the restriction of the CMSIOC to the maximum strain increment outside of the crack explains its better performance than the CMSI. Indeed, during the whole simulation, the maximum strain increment is always located in the cracking region ahead of the crack tip. Hence, the CMSIOC always acts as a crack growth limiter. Although we do not detail these results, the CMTCDV demonstrated limited robustness in our numerical experiments, as small values of the step $\Delta\tau_{imp}$ occasionally lead to stagnation in the load step. Thus, we choose to use the CMSIOC path-following methods in the following studies.

### 4.2.3. Convergence with bulk mesh size

Having identified path-following methods that successfully capture the equilibrium path, we now focus on a convergence study with respect to the bulk mesh size $h_{max}$. The regularization length is fixed at $\ell = L/64$, and the mesh size along the crack path is set to $h = \ell/10$. The bulk mesh size is defined as $h_{max} = L/N_{max}$, where $N_{max} \in \{64, 128, 256\}$. With these parameters, the mesh resolution along the crack path is sufficiently refined to ensure convergence. However, $\Gamma$-convergence may not yet be fully achieved, implying that the phase-field model may not perfectly align with the sharp crack model. Note that the numerical $\Gamma$-convergence is tackled in Section 4.2.4.

Figure 10 presents both the force-displacement and load factor-crack length curves for the different bulk mesh sizes $h_{max}$. The force-displacement curves in Figure 10a show that both the phase-field and sharp crack models exhibit similar convergence behavior with respect to the bulk mesh size $h_{max}$. For both models, the bulk mesh size affects the overall structural stiffness resolution. The remaining difference between the LEFM and the phase-field models arises from the incomplete $\Gamma$-convergence affecting the accuracy of the phase-field model.

Surprisingly, Figure 10b shows that the load factor-crack length curve remains unaffected by variations in the bulk mesh size $h_{max}$. This behavior might be explained by recalling that the load factor at the propagation threshold is

$$ \lambda = \sqrt{\frac{G_c}{G}}, \quad \text{where} \quad G = -\frac{\partial P}{\partial a} = -\frac{1}{2} \frac{\partial R}{\partial a} u^2, \qquad (37) $$

where R represents the overall structural rigidity of the specimen. This indicates that the load factor $\lambda$ depends on the derivative of the structural rigidity R with respect to the crack length a. From those results, it seems that even a coarse bulk mesh is sufficient to accurately resolve this derivative, thereby leaving the load factor-crack length curve (nearly) unchanged.
<figure><img src="image_15.png" /><figcaption>Figure 10: Equilibrium paths of the CT specimen with a secondary crack obtained with the phase-field fracture model for varying bulk mesh sizes $h_{\max}$ (CMSIOC path-following method).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 11: Equilibrium path of the CT specimen with a secondary crack obtained with the phase-field fracture model for different regularization length $\ell$ (CMSIOC path-following method).</figcaption></figure>

#### 4.2.4. Γ-convergence

To conclude this study, we numerically assess the Γ-convergence of the phase-field fracture model toward the sharp crack model. For this analysis, the bulk mesh size is fixed at $h_{\max} = L/128$, while the mesh size along the expected crack path is set to $h = \ell/10$. The regularization length $\ell$ is varied as $\ell = L/N$, with $N \in \{64, 128, 256\}$. Since the CMSIOC method is employed, the step size $\Delta\tau_{\text{imp}}$ is adapted with the mesh size to maintain an approximately consistent crack growth per load step. In this example, the step size is adjusted according to $\Delta\tau_{\text{imp}} = \frac{N}{64} \times 0.01$.

Figure 11 shows the equilibrium paths obtained for the different values of the regularization length $\ell$. As expected, the results demonstrate convergence of the phase-field fracture model toward the sharp crack model in both the force-displacement and load factor-crack length curves. A basic, but important, remark is that we are practically limited in our capacity to achieve Γ-convergence by our computational resources, as the crack profile must be sufficiently refined.
<figure><img src="image_17.png" /><figcaption>Figure 12: Equilibrium path of the CT specimen with an eccentric hole obtained with the phase-field fracture model for different regularization length $ℓ$ (CMSIOC path-following method).</figcaption></figure>

### 4.3. Example 3: CT problem under force load with an eccentric hole (mixed mode I+II)

The third numerical example, described in Figure 3c, analyzes a pre-cracked CT specimen containing a hole, reproducing an experimental setup detailed by Cavuoto et al. (2022). The presence of the hole induces mode mixity (I+II) as the crack propagates, resulting in a curved crack path. Following the work of Triclot et al. (2023), this example has to be modeled using force boundary conditions to obtain reliable results. As it presents a structural softening, the solution of this problem requires using a path-following method to maintain equilibrium during the whole crack propagation, and avoid the loss of force balance. The primary objective of this example is to evaluate whether the variational phase-field fracture model can accurately reproduce the results of the LEFM sharp crack model in a mixed-mode (I+II) crack propagation scenario.

#### 4.3.1. Implementation details

The phase-field simulations are again conducted on an unstructured mesh to minimize mesh-induced bias (Negri, 2003, 2007), which could otherwise influence the predicted crack path. Both the phase-field and sharp crack models employ a uniform bulk mesh size of $h_{max} = W/128$, while the mesh size in the expected crack region is refined to $h = ℓ/8$. Given its demonstrated efficiency in prior examples, the CMSIOC path-following method is exclusively employed. Furthermore, identical crack propagation paths were obtained using the LEFM sharp crack model with two distinct crack increments: $\Delta a = W/256$ and $\Delta a = W/512$. This convergence demonstrates that the smaller increment ($\Delta a = W/512$) yields sufficiently refined results, ensuring the accuracy of the LEFM reference simulations.

#### 4.3.2. Γ-convergence

Figure 12 illustrates the $F-u$ and $λ-a$ curves for varying regularization lengths $ℓ$. In contrast to previous examples, convergence with respect to $ℓ$ is more demanding. More precisely, it is achieved at much smaller values of $ℓ$, as the crack path itself must be accurately predicted. Misprediction of the crack path directly influences the force-displacement relationship, leading to a more complex convergence behavior. Despite this increased complexity, the phase-field model demonstrates convergence toward the LEFM sharp crack model, confirming its validity in capturing the essential features of crack propagation under mixed-mode loading conditions.

Let us now evaluate the convergence of the crack path predicted by the phase-field model. Figure 13 compares the phase-field crack path (region where $α > 0.8$) against the crack path derived from the LEFM sharp crack model. The phase-field crack path exhibits a strong dependence on the regularization length $ℓ$ when $ℓ$ is too large, and the model fails to resolve the sharp crack trajectory accurately. As $ℓ$ decreases, the phase-field model gradually converges
<figure><img src="image_18.png" /><figcaption>(a) Crack path (full specimen)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Crack path (zoom)</figcaption></figure>

Figure 13: Comparison of the crack path, in the CT specimen with an eccentric hole, predicted by the sharp crack model (black) and using the phase-field model (region where $\alpha \ge 0.8$) for different regularization lengths. The phase-field simulations are conducted using varying regularization lengths $l$ and the CMSIOC path-following method.
<table><caption>Table 4: Summary of the comparison of different path-following methods (*: Depends on control points choice, —: Not applicable).</caption><thead><tr><th></th><th>No PF</th><th>CNDI</th><th>CLAL</th><th>CCDV</th><th>CMTCDV</th><th>CMSI</th><th>CMSIOC</th></tr></thead><tbody><tr><td>Equilibrium path</td><td>✗</td><td>✓*</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Problem indep.</td><td>—</td><td>✗</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Model indep.</td><td>—</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td></tr></tbody></table>

toward the LEFM solution. However, even the most converged phase-field crack path remains slightly offset from the LEFM crack path. This discrepancy might be attributed to numerical errors introduced during the discretization of the variational phase-field fracture model, including mesh-induced bias (Negri, 2003, 2007).

# 5. Discussion

In this section, we critically examine the results of the numerical examples from multiple perspectives. In Section 5.1, we discuss the various path-following methods presented and assessed in this work, intending to provide practical recommendations to practitioners. In Section 5.2, we discuss extensions of these methods to more complex loading scenarios. Finally, in Section 5.3, we examine the usefulness of path-following methods for phase-field fracture simulations.

## 5.1. Assessment of the different path-following methods

### 5.1.1. Overview of the methods

In this study, we evaluated six distinct path-following methods for phase-field fracture simulations. Among these, we introduced the Control by Maximum Strain Increment Outside the Crack (CMSIOC), which, to the best of our knowledge, has not been previously documented in the literature. We selected these methods as they require no modifications to the elastic or phase-field sub-problems or their respective solvers. The only addition to the classic alternate minimization scheme is an extra step to solve the path-following control equation, which remains decoupled from the other sub-problems. This property makes their implementation relatively straightforward compared to the overall complexity of the phase-field solver. The different methods were implemented using FEniCSx, in the open-source code pragma (Loiseau and Lazarus, 2025a), and applied to simulate crack propagation in three numerical examples of increasing complexity. This approach provides a comprehensive assessment of their respective capabilities. A summary of our assessment is presented in Table 4 and is discussed in more detail below. Specifically, it shows that only the CMSI and the CMSIOC could determine the equilibrium path while being independent of both the problem and the constitutive model.

### 5.1.2. Performance and limitations of each method

All methods captured the equilibrium path, but the CLAL and CCDV methods failed to capture snapback behavior in our numerical examples fully. This limitation renders them less suitable for general engineering applications. The CNDI method also struggled in our tests, as its success heavily depends on the choice of control points. While this choice may be intuitive for simple academic problems or standard tests, it becomes increasingly challenging for more complex scenarios, such as multi-cracking. We now focus on the methods that successfully captured the equilibrium path in all examples: CMSI, CMTCDV, and CMSIOC. The CMTCDV method depends on the specific constitutive model used, making it less versatile than CMSI or CMSIOC. Additionally, its implementation is more complex, requiring a dedicated (scalar) nonlinear solver. We also observed, though not rigorously documented here, that CMTCDV can lead to stagnation of the load factor if the step size is poorly chosen, halting any crack growth in the simulation.

Both CMSI and CMSIOC performed well in our numerical tests. They are independent of both the problem and the constitutive model (note that CMSIOC reduces to CMSI for models without cracks or damage). However, they share two drawbacks: their step size must be adapted to the mesh size, and they require the implementation of the nested interval algorithm. In theory, another drawback is that the method could follow an unloading elastic path, which never occurred in the numerical examples. We consider these drawbacks minor compared to their numerous advantages. The
<figure><img src="image_20.png" /><figcaption>Figure 14: Comparison of the normalized control quantities for the CMSI and the CMSIOC for a long crack length.</figcaption></figure>

key feature distinguishing CMSIOC from CMSI is its ability to provide a nearly constant crack growth rate per load step, which is highly advantageous from a practical standpoint.

### 5.1.3. Why does the CMSIOC work well?

The main idea behind the CMSIOC is to control the crack increment by limiting the maximum strain increment outside the cracked regions ($\alpha = 1$), as these regions significantly impact the performance of the CMSI. For illustration, we consider the example from Section 4.2. The Figure 14 compares the normalized control quantity for both the CMSI and the CMSIOC. For this load step, the maximum strain increment is located inside the cracked region close to the initial crack tip. This strain increment is proportional to the Crack Opening Displacement (COD), and the CMSI becomes a control in COD. This is what occured with the CMSI in the Section 4.2.2, which lead to a strong decrease in the crack growth rate and required a lot of load step to reach full failure. However, we want direct control over the crack growth, rather than over its opening. As shown in Figure 14, excluding the cracked regions ($\alpha = 1$) from the control quantity in the CMSIOC moves the maximum strain increment to the region near the (regularized) stress singularity, ahead of the current crack tip. This maximum strain increment outside the crack reflects the localized high-strain field induced by crack growth ahead of the current crack tip. By limiting this increment, the CMSIOC method effectively controls crack growth.

### 5.1.4. Recommendations

To conclude on this assessment of the different path-following methods, we recommend using either the CNDI or the CMSIOC. The CNDI is advantageous for well-known problems (e.g., to simulate standard experiments), as it requires some *a priori* knowledge, or at least some strong intuition, to choose the control points. Despite this choice of control points being problem-dependent, their simplicity of implementation makes the CNDI an interesting approach. The CMSIOC, while requiring the implementation of the nested interval algorithm, is very appealing thanks to its generality (both problem and constitutive model independent), robustness, and is very practical as it yields a nearly constant crack growth per load step.

## 5.2. Extension to more complex loading scenarios

In this study, we focused on boundary conditions with proportional displacements or tractions. However, other load cases can be treated with no or minor changes. For instance, the extension to controlled body forces is straightforward. This section discusses the applicability of the proposed methods to two more complex loading scenarios: cases with combined constant and proportional loads, and cases with non-proportional loads.
### 5.2.1. Combining constant and proportional loads

An extension of interest involves scenarios where a constant load is combined with a proportional load (e.g., gravity combined with an applied force). In such cases, the displacement field can be decomposed using the principle of superposition as $u = \lambda \bar{u} + u_{\text{const}}$, thanks to the linearity of the elastic problem. This decomposition requires solving two elastic sub-problems: one for the proportional component $\bar{u}$ and another for the constant component $u_{\text{const}}$. While the control equations for the path-following methods require minor adjustments, no fundamental changes are needed for their solution.

### 5.2.2. Non-proportional loads

A more challenging extension involves non-proportional loads, i.e. loads where $u_{\text{imp}}(x, \lambda) \neq \lambda \bar{u}_{\text{imp}}(x)$. For instance, one could consider a moving load translating along a structure, described by $u_{\text{imp}}(x, \lambda) = u'_{\text{imp}}(x - \lambda L)$. In this case, the linearity between the displacement field $u$ and the load factor $\lambda$ no longer holds, i.e., $u(x) \neq \lambda \bar{u}(x)$. However, the load can be linearized using a Taylor series expansion:

$$ u_{\text{imp}}(\lambda_0 + \Delta\lambda) = u_{\text{imp}}(\lambda_0) + \Delta\lambda \frac{\partial u_{\text{imp}}(\lambda)}{\partial \lambda} \Rightarrow u(x) = u_0(x) + \Delta\lambda \bar{u}(x), \quad (38) $$

which provides a linear relationship between the displacement field $u$ and the increment in load factor $\Delta\lambda$. The methods proposed in this study can be adapted by initializing the simulation with an estimate of the first critical load factor $\lambda_0$ and then solving the control equations for $\Delta\lambda$ during the simulation. Then, the displacement field decomposes as $u(x) = u_{\lambda_0} + \Delta\lambda \bar{u}$, so the modifications to the control equations are analogous to those required for a constant load. When using this linearization technique, it is essential to ensure that the increments in load factor $\Delta\lambda$ remain sufficiently small to maintain accuracy.

## 5.3. When and why use path-following methods?

Path-following methods are essential tools for accurately characterizing the equilibrium path of failing structures (here with initial cracks). They enable determining the sequence of critical loads leading to the failure of the structure. Several aspects of their usefulness are discussed below.

### 5.3.1. To ensure robustness and stability of the numerical resolution

From a purely numerical point of view, enforcing quasi-static crack propagation enhances the robustness of numerical solvers. Path-following methods enable the simulation of unstable phenomena such as snapbacks and snap-throughs. They also improve the stability of the staggered solver and distribute the computational effort more evenly across load steps (see Figure 9). To further enhance numerical stability, adaptive step-size schemes can also be employed to maintain a nearly constant number of iterations per load step or ensure uniform crack increment sizes (see the review of Rastiello et al. (2022)). However, it is worth noting that path-following methods typically result in longer simulation times compared to unstable counterparts. As will be discussed below, this trade-off is justified by the significantly greater wealth of information they provide about the structural response.

### 5.3.2. To ensure self-consistency by forcing the quasi-static assumption

Fracture models, as well as any model exhibiting strain-softening behavior, are prone to structural instabilities such as snapbacks. During these phases, unstable crack propagation occurs if the energy release rate $G$ exceeds the critical value $G_c$, violating the Griffith criterion (Maggiorelli and Negri, 2025). In quasi-static phase-field simulations, this violation leads to significant crack growth during a single load increment. The energy minimization is then performed over the entire crack jump rather than incremental crack growth, which undermines the physical validity of the quasi-static framework. There is no guarantee that a single local minimization step over a large crack increment yields the same state as a sequence of local minimizations over small increments. This might be particularly critical for simulating crack propagation in heterogeneous media where local pinning on tougher zones may occur (Roux et al., 2003; Vasoya et al., 2016; Henry, 2024).
<figure><img src="image_21.png" /><figcaption>Figure 15: Illustration of the energy balance during a snapback instability.</figcaption></figure>

### 5.3.3. To reveal potential catastrophic dynamical phases

The sequence of critical loads reveal potential dynamic crack propagation phases for a given loading scenario together with their associated critical loads. For instance, coming back to the example of the introduction (see Figure 1), applying an increasing displacement would result in (i) a non-propagation phase until A, (ii) quasistatic propagation between A and C, (iii) dynamic merging up to the second crack, (iv) quasistatic propagation from C to D, (v) catastrophic failure of the remaining ligament. In contrary applying an increasing force would result in a non-propagation phase followed by a catastrophic failure above A. In other words, this sequence provides a complete picture of the structural response under specific loading conditions. This capability is particularly valuable for assessing structural instability (i.e., potential catastrophic failure), which is important in the design and analysis of structures.

### 5.3.4. To handle force boundary conditions

Path-following methods are particularly advantageous in scenarios where applied loads cannot be equilibrated. For instance, incremental force boundary conditions, often required to replicate experimental setups (Triclot et al., 2023), can lead to numerical challenges if a propagating crack isolates the loaded region. In such cases, no cohesive forces remain to equilibrate the applied load, resulting in a loss of problem well-posedness and solver failure. As shown in the numerical examples of Section 4.2 and Section 4.3, path-following methods address this issue by adjusting the load factor to recover equilibrium, enabling the simulation to capture the entire crack propagation process up to structural failure.

### 5.3.5. To design optimal tougher materials

The path-following method introduced up to now, is strictly speaking limited to low crack speed, for which the inertial effect can be neglected (Freund, 1998). However, introducing material heterogeneities to induce dynamical snapback phases is known to be an efficient way to increase the effective fracture toughness by increasing the energy needed to break the material (Roux et al., 2003; Morano et al., 2020; Lebihain et al., 2021; Triclot et al., 2024). An accurate determination of the energy dissipated during these phases is a prerequisite for optimal patterning of these heterogeneities. This energy can be estimated using the force-displacement equilibrium curves. Taking as illustration the typical example of Figure 15, it corresponds to the blue area. For the propagation of a single crack it can be done using classical LEFM sharp crack models. As a perspective, the phase-field path following methods introduced here, will permit to study situations where more complex patterns, as multi-cracking for instance, are involved.

## 6. Conclusion

This study focuses on determining the equilibrium path of structures in the context of quasi-static crack propagation problems. We employed two fracture models: a LEFM sharp crack model based on Griffith's theory, and a variational
phase-field fracture model. For the sharp crack model, we proposed and implemented a path-following algorithm, which served as our reference for evaluating the performance of the phase-field model. With the phase-field model, we implemented several path-following methods, including a novel approach: the CMSIOC. Inspired by the work of Chen and Schreyer (1990), it limits the maximum strain increment outside the crack. It corresponds to controlling the strain increment in the localized high-strain ahead of the crack tip, which is induced by crack growth. All the path-following methods studied in this work were selected for their simplicity of implementation. Indeed, all can be integrated into phase-field fracture simulations with only minor modifications to the classic staggered solver based on alternate minimization. Through three numerical examples, we demonstrated that the CMSIOC is independent of both constitutive model and problem specifics, and successfully captures the equilibrium path in all numerical examples. Additionally, it provides a nearly constant crack increment per load step, which is often considered ideal. Our findings also confirmed that, once mesh convergence and Γ-convergence are achieved, the variational phase-field model closely approximates the predictions of the LEFM sharp crack model.

Beyond comparing the performance of path-following methods, this work also highlights their advantages from both numerical and physical considerations. Numerically, these methods enhance solver stability, especially in scenarios where structural instability may arise. As a result, path-following methods enable the determination of the sequence of critical states and instability phases within a structure, providing a comprehensive understanding of its response. From a physical standpoint, these methods enforce the quasi-static assumption by ensuring incremental crack growth, which guarantees physically meaningful incremental crack propagation. They also enable the analysis of problems involving structural softening and the application of force boundary conditions. Such problems cannot be addressed using classic incremental load control, as the force balance cannot be maintained once crack propagation initiates without reducing the load. Last but not least, they enable precise quantification of the energy dissipated during a snapback instability paving the way for the rational optimal design of material heterogeneities in order to maximize effective macroscopic fracture resistance.

## Acknowledgements

The authors would like to thank Breno Ribeiro Nogueira and Alexandre Daby-Seesaram for their valuable discussions and insights regarding path-following methods.

For the purpose of Open Access, a CC-BY public copyright licence (https://creativecommons.org/licenses/by/4.0/) has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## Fundings

The work was supported by Agence de l'Innovation de Défense - AID - via Centre Interdisciplinaire d'Etudes pour la Défense et la Sécurité - CIEDS - (projects 2022 - FracAddi), and by Agence Nationale de la Recherche (ANR-23-CE51-0054 3FAM).

## Appendix A. Systematic literature review on path-following methods applied to phase-field fracture

To assess the application of path-following methods in phase-field fracture, a systematic literature review was conducted. The search was performed on December 12, 2025, using the Scopus database with the following query:

TITLE-ABS-KEY ( "load control" OR "path following" OR "path-following" OR "arc length" OR "arc-length" OR "continuation") AND
<table><caption>Table A.5: Summary of study on phase-field using path-following methods. The references are organized based on the research group that conducted the study. The references with a star does not appear in the systematic review and have been added by the authors.</caption><thead><tr><th>Reference</th><th>Control quantity</th><th>Irreversibility</th><th>Scheme</th></tr></thead><tbody><tr><td>Vignollet et al. (2014)*</td><td>Energy release (Gutiérrez, 2004)</td><td>History variable</td><td>Both</td></tr><tr><td>May et al. (2016)</td><td>Energy release (Gutiérrez, 2004)</td><td>History variable</td><td>Staggered</td></tr><tr><td>Singh et al. (2016)</td><td>Fracture dissipation (Verhoosel et al., 2009)</td><td>History variable</td><td>Both</td></tr><tr><td>Duda et al. (2015)*</td><td>Nodal displacements (De Borst, 1987)</td><td>History variable</td><td>Staggered</td></tr><tr><td>Zambrano et al. (2023)</td><td>Crack driving variable</td><td>History variable</td><td>Staggered</td></tr><tr><td>Wu (2018)</td><td>Nodal displacements (De Borst, 1987)</td><td>Variational inequality</td><td>Staggered</td></tr><tr><td></td><td>Energy release (May et al., 2016)</td><td>Variational inequality</td><td>Staggered</td></tr><tr><td>Raghu et al. (2019)</td><td>Classic arc-length</td><td>History variable</td><td>Staggered</td></tr><tr><td>Wambacq et al. (2021)</td><td>Energy release (Gutiérrez, 2004)</td><td>Variational inequality</td><td>Monolithic</td></tr><tr><td>León Baldelli and Maurini (2021)</td><td>Second-order continuation algorithm</td><td>Variational inequality</td><td>Staggered</td></tr><tr><td>Liu et al. (2022)</td><td>Empirical load adaptation (based on Δα<sub>max</sub>)</td><td>History variable</td><td>Staggered</td></tr><tr><td>Wu et al. (2022)<sup>3</sup></td><td>De Borst (1987)</td><td>History variable</td><td>Staggered</td></tr><tr><td>Bharali et al. (2022)</td><td>Fracture dissipation Singh et al. (2016)</td><td>History variable</td><td>Monolithic</td></tr><tr><td>Bharali et al. (2025)</td><td>Fracture dissipation Singh et al. (2016)</td><td>History variable</td><td>Both</td></tr><tr><td>Liu and Moran (2023)</td><td>Classic arc-length</td><td>History variable</td><td>Monolithic</td></tr><tr><td>Börjesson et al. (2022)*</td><td>Fracture dissipation (Verhoosel et al., 2009)</td><td>History variable</td><td>Monolithic</td></tr><tr><td>Yang and Shen (2023)</td><td>Increment of internal energy</td><td>History variable</td><td>Staggered</td></tr><tr><td>Ma et al. (2024)</td><td>Singh et al. (2016)</td><td>-</td><td>-</td></tr><tr><td>Aranda and Segurado (2025a)*</td><td>Crack length (close to Singh et al. (2016))</td><td>History variable</td><td>Monolithic</td></tr><tr><td>Aranda and Segurado (2025b)</td><td>(Aranda and Segurado, 2025a)</td><td>History variable</td><td>Monolithic</td></tr><tr><td>Levy et al. (2025)</td><td>(León Baldelli and Maurini, 2021)</td><td>Variational inequality</td><td>Staggered</td></tr></tbody></table>
"phase-field"
AND
"fracture")

The results of this search are summarized in Table A.5. Additionally, a few relevant references marked with an asterisk (*) were included, as they were not captured by the initial query but were deemed pertinent to this review. Several studies were also excluded for specific reasons:

* The work of Almi and Negri (2020) was excluded because it employs arc-length methods in a purely theoretical context without practical implementation.

* Russ and Waisman (2020) was deemed outside the scope of this review.

* Wu (2021) was excluded due to our inability to read the full text.

Overall, this review reveals that the use of path-following methods in the phase-field fracture domain is low compared to the available literature in this domain. Still, the existing studies cover a diverse range of control equations for both staggered and monolithic numerical schemes. However, it is notable that most of these methods either lack generality (problem or model dependency), or necessitate substantial modifications to classic phase-field solvers.

## Appendix B. Nested Interval Algorithm

Consider a control equation of the form

$$ \max_{x \in \Omega} (\lambda c_1(x) - c_0(x)) - \Delta \tau_{\text{imp}} = 0. \qquad \text{B.1} $$

This form of equation often appears when solving constraint equations based on a local control quantity. Note that this equation may have zero or multiple solutions. To find an interval $[\lambda_{\min}, \lambda_{\max}]$ that contains all possible solutions, we employ the nested interval algorithm proposed by Lorentz and Badel (2004). The procedure consists of the following steps:

1. Initialize the interval with $\lambda_{\min} = -\infty$ and $\lambda_{\max} = +\infty$.

2. Loop over the points $x$ (i.e., the nodes).

(a) Compute the load factor associated with this element $\lambda(x) = \frac{c_0(x)+\Delta\tau_{\text{imp}}}{c_1(x)}$,

(b) If $c_1(x) > 0$, then update the upper bound: $\lambda_{\max} = \max(\lambda_{\max}, \lambda(x))$,

(c) If $c_1(x) < 0$, then update the lower bound: $\lambda_{\min} = \min(\lambda_{\min}, \lambda(x))$.

For simplicity, in this work, we always select the load factor as $\lambda = \lambda_{\max}$. If the interval is empty, the procedure should be restarted with a smaller step size. Note that this case never occurred in our numerical simulations.

## References

Almi, S., Negri, M., 2020. Analysis of Staggered Evolutions for Nonlinear Energies in Phase Field Fracture. Archive for Rational Mechanics and Analysis 236, 189–252. URL: https://doi.org/10.1007/s00205-019-01468-4, doi:10.1007/s00205-019-01468-4.

Ambrosio, L., Fusco, N., Pallara, D., 2000. Functions of Bounded Variation and Free Discontinuity Problems. Oxford Mathematical Monographs, Oxford University Press, Oxford, New York.

Ambrosio, L., Tortorelli, V.M., 1990. Approximation of functional depending on jumps by elliptic functional via t-convergence. Communications on Pure and Applied Mathematics 43, 999–1036. URL: http://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160430805, doi:10.1002/cpa.3160430805.

Ambrosio, L., Tortorelli, V.M., 1992. On the Approximation of Free Discontinuity Problems. Bollettino Della Unione Matematica Italiana 6B, 105–123. WOS:A1992HW06100006.

Amestoy, M., Leblond, J.B., 1992. Crack paths in plane situations—II. Detailed form of the expansion of the stress intensity factors. International Journal of Solids and Structures 29, 465–501. URL: https://www.sciencedirect.com/science/article/pii/002076839290210K, doi:10.1016/0020-7683(92)90210-K.

Aranda, P., Segurado, J., 2025a. A Crack-Length Control Technique for Phase-Field Fracture in FFT Homogenization. International Journal for Numerical Methods in Engineering 126, e7664. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.7664, doi:10.1002/nme.7664. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.7664.
Aranda, P., Segurado, J., 2025b. Effective toughness estimation by FFT based phase field fracture: Application to composites and poly-crystals. Mechanics Research Communications 149, 104498. URL: https://www.sciencedirect.com/science/article/pii/S0093641325001314, doi:10.1016/j.mechrescom.2025.104498.

Baratta, I.A., Dean, J.P., Dokken, J.S., Habera, M., Hale, J.S., Richardson, C.N., Rognes, M.E., Scroggs, M.W., Sime, N., Wells, G.N., 2023. DOLFINx: The next generation FEniCS problem solving environment. URL: https://zenodo.org/records/10447666, doi:10.5281/zenodo.10447666.

Bellettini, G., Coscia, A., 1994. Discrete approximation of a free discontinuity problem. Numerical Functional Analysis and Optimization 15, 201–224. URL: https://doi.org/10.1080/01630569408816562, doi:10.1080/01630569408816562. publisher: Taylor & Francis _eprint: https://doi.org/10.1080/01630569408816562.

Bharali, R., Goswami, S., Anitescu, C., Rabczuk, T., 2022. A robust monolithic solver for phase-field fracture integrated with fracture energy based arc-length method and under-relaxation. Computer Methods in Applied Mechanics and Engineering 394, 114927. URL: https://www.sciencedirect.com/science/article/pii/S0045782522001992, doi:10.1016/j.cma.2022.114927.

Bharali, R., van der Meer, F.P., Larsson, F., Jänicke, R., 2025. A time step-size computing arc-length method for the phase-field hydraulic fracture model. Computer Methods in Applied Mechanics and Engineering 436, 117687. URL: https://www.sciencedirect.com/science/article/pii/S0045782524009411, doi:10.1016/j.cma.2024.117687.

Börjesson, E., Remmers, J.J.C., Fagerström, M., 2022. A generalised path-following solver for robust analysis of material failure. Computational Mechanics 70, 437–450. URL: https://doi.org/10.1007/s00466-022-02175-w, doi:10.1007/s00466-022-02175-w.

Bourdin, B., 1999. Image segmentation with a finite element method. ESAIM: Mathematical Modelling and Numerical Analysis 33, 229–244. URL: https://www.esaim-m2an.org/articles/m2an/abs/1999/02/m2an811/m2an811.html, doi:10.1051/m2an:1999114. number: 2 Publisher: EDP Sciences.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2000. Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids 48, 797–826. URL: http://www.sciencedirect.com/science/article/pii/S0022509699000289, doi:10.1016/S0022-5096(99)00028-9.

Bourdin, B., Francfort, G.A., Marigo, J.J., 2008. The Variational Approach to Fracture. Journal of Elasticity 91, 5–148. URL: https://doi.org/10.1007/s10659-007-9107-3, doi:10.1007/s10659-007-9107-3.

Braides, A., 1998. Approximation of Free-Discontinuity Problems. volume 1694 of Lecture Notes in Mathematics. Springer, Berlin, Heidelberg. URL: http://link.springer.com/10.1007/BFb0097344, doi:10.1007/BFb0097344.

Cavuoto, R., Lenarda, P., Misseroni, D., Paggi, M., Bigoni, D., 2022. Failure through crack propagation in components with holes and notches: An experimental assessment of the phase field model. International Journal of Solids and Structures 257, 111798. URL: https://www.sciencedirect.com/science/article/pii/S0020768322002852, doi:10.1016/j.ijsolstr.2022.111798.

Chambolle, A., 2004. An approximation result for special functions with bounded deformation. Journal de Mathématiques Pures et Appliquées 83, 929–954. URL: https://www.sciencedirect.com/science/article/pii/S0021782404000285, doi:10.1016/j.matpur.2004.02.004.

Chen, Z., Schreyer, H.L., 1990. A numerical solution scheme for softening problems involving total strain control. Computers & Structures 37, 1043–1050. URL: https://www.sciencedirect.com/science/article/pii/004579499090016U, doi:10.1016/0045-7949(90)90016-U.

Crisfield, M., 1982. Local instabilities in the non-linear analysis of reinforced concrete beams and slabs. Proceedings of the Institution of Civil Engineers 73, 135–145. URL: https://www.icevirtuallibrary.com/doi/10.1680/iicep.1982.1876, doi:10.1680/iicep.1982.1876. publisher: ICE Publishing.

Crisfield, M.A., 1981. A fast incremental/iterative solution procedure that handles “snap-through”. Computers & Structures 13, 55–62. URL: https://www.sciencedirect.com/science/article/pii/0045794981901085, doi:10.1016/0045-7949(81)90108-5.

De Borst, R., 1987. Computation of post-bifurcation and post-failure behavior of strain-softening solids. Computers & Structures 25, 211–224. URL: https://linkinghub.elsevier.com/retrieve/pii/0045794987901441, doi:10.1016/0045-7949(87)90144-1.

Destuynder, P., Djaoua, P.E.M., Chesnay, L., 1981. Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile. Mathématical Methods in the Applied Sciences 3, 70–87. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/mma.1670030106, doi:10.1002/mma.1670030106. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/mma.1670030106.

Duda, F.P., Ciarbonetti, A., Sánchez, P.J., Huespe, A.E., 2015. A phase-field/gradient damage model for brittle fracture in elastic–plastic solids. International Journal of Plasticity 65, 269–296. URL: https://www.sciencedirect.com/science/article/pii/S0749641914001892, doi:10.1016/j.ijplas.2014.09.005.

Erdogan, F., Sih, G.C., 1963. On the Crack Extension in Plates Under Plane Loading and Transverse Shear. Journal of Basic Engineering 85, 519–525. URL: https://doi.org/10.1115/1.3656897, doi:10.1115/1.3656897.

Farrell, P., Maurini, C., 2017. Linear and nonlinear solvers for variational phase-field models of brittle fracture. International Journal for Numerical Methods in Engineering 109, 648–667. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.5300, doi:10.1002/nme.5300. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.5300.

Fayezioghani, A., Vandoren, B., Sluys, L.J., 2019. A posteriori performance-based comparison of three new path-following constraints for damage analysis of quasi-brittle materials. Computer Methods in Applied Mechanics and Engineering 346, 746–768. URL: https://www.sciencedirect.com/science/article/pii/S0045782518304584, doi:10.1016/j.cma.2018.09.014.

Francfort, G.A., Marigo, J.J., 1998. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids 46, 1319–1342. URL: http://www.sciencedirect.com/science/article/pii/S0022509698000349, doi:10.1016/S0022-5096(98)00034-9.

Freund, L.B., 1998. Dynamic Fracture Mechanics. Cambridge University Press. Google-Books-ID: Cij0jwEACAAJ.

Geers, M.G.D., 1999. Enhanced solution control for physically and geometrically non-linear problems. Part I—the subplane control approach. International Journal for Numerical Methods in Engineering 46, 177–204. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/\%28SICI\%291097-0207\%2819990920\%2946\%3A2\%3C177\%3A\%3AAID-NME668\%3E3.0.CO\%3B2-L, doi:10.1002/(SICI)1097-0207(19990920)46:2<177::AID-NME668>3.0.CO;2-L. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/\%28SICI\%291097-0207\%2819990920\%2946\%3A2\%3C177\%3A\%3AAID-NME668\%3E3.0.CO\%3B2-L.
Giacomini, A., Ponsiglione, M., 2003. A Discontinuous Finite Element Approximation of Quasi-static Growth of Brittle Fractures. Numerical Functional Analysis and Optimization 24, 813–850. URL: https://doi.org/10.1081/NFA-120026378, doi:10.1081/NFA-120026378. publisher: Taylor & Francis _eprint: https://doi.org/10.1081/NFA-120026378.

Griffith, A.A., 1920. The phenomena of rupture and flow in solids. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character 221, 163–198. URL: http://royalsocietypublishing.org/doi/abs/10.1098/rsta.1921.0006, doi:10.1098/rsta.1921.0006.

Gutiérrez, M.A., 2004. Energy release control for numerical simulations of failure in quasi-brittle solids. Communications in Numerical Methods in Engineering 20, 19–29. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cnm.649, doi:10.1002/cnm.649. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cnm.649.

Henry, H., 2024. Pinning of crack fronts by hard and soft inclusions: A phase field study. Physical Review E 109, 025002. URL: https://link.aps.org/doi/10.1103/PhysRevE.109.025002, doi:10.1103/PhysRevE.109.025002. publisher: American Physical Society.

Hussain, M., Pu, S., Underwood, J., 1973. Strain Energy Release Rate for a Crack Under Combined Mode I and Mode II, in: 1973 National Symposium on Fracture Mechanics, Part II, ASTM International, Washington, United States. pp. 1–27. URL: https://hal.science/hal-04166069, doi:10.1520/STP33130S. issue: 2.

Irwin, G.R., 1957. Analysis of Stresses and Strains Near the End of a Crack Traversing a Plate. Journal of Applied Mechanics 24, 361–364. URL: https://doi.org/10.1115/1.4011547, doi:10.1115/1.4011547.

Lazarus, V., 2003. Brittle fracture and fatigue propagation paths of 3D plane cracks under uniform remote tensile loading. International Journal of Fracture 122, 23–46. URL: https://doi.org/10.1023/B:FRAC.0000005373.73286.5d, doi:10.1023/B:FRAC.0000005373.73286.5d.

Lebihain, M., Ponson, L., Kondo, D., Leblond, J.B., 2021. Effective toughness of disordered brittle solids: A homogenization framework. Journal of the Mechanics and Physics of Solids 153, 104463. URL: https://www.sciencedirect.com/science/article/pii/S0022509621001393, doi:10.1016/j.jmps.2021.104463.

León Baldelli, A.A., Maurini, C., 2021. Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. Journal of the Mechanics and Physics of Solids 152, 104424. URL: https://www.sciencedirect.com/science/article/pii/S0022509621001010, doi:10.1016/j.jmps.2021.104424.

Levy, M., Vicentini, F., Yosibash, Z., 2025. PFM for crack nucleation and propagation in isotropic heterogeneous materials. Computational Mechanics URL: https://doi.org/10.1007/s00466-025-02711-4, doi:10.1007/s00466-025-02711-4.

Liu, X., Réthoré, J., Lubrecht, A.A., 2022. An efficient matrix-free preconditioned conjugate gradient based multigrid method for phase field modeling of fracture in heterogeneous materials from 3D images. Computer Methods in Applied Mechanics and Engineering 388, 114266. URL: https://www.sciencedirect.com/science/article/pii/S0045782521005764, doi:10.1016/j.cma.2021.114266.

Liu, Y., Moran, B., 2023. Tearing a neo-Hookean sheet. Part I: insights into the crack tip fields by a phase field damage model. International Journal of Fracture 240, 71–91. URL: https://doi.org/10.1007/s10704-022-00676-6, doi:10.1007/s10704-022-00676-6.

Loiseau, F., Lazarus, V., 2025a. `fragma`: A finite element solver for simulating crack propagation using the phase-field approach to fracture. URL: https://doi.org/10.5281/zenodo.15654994, doi:10.5281/zenodo.15654995.

Loiseau, F., Lazarus, V., 2025b. How to introduce an initial crack in phase field simulations to accurately predict the linear elastic fracture propagation threshold? Journal of Theoretical, Computational and Applied Mechanics URL: https://jtcam.episciences.org/15987, doi:10.46298/jtcam.15198. publisher: Episciences.org.

Loiseau, F., Lazarus, V., 2026. `gcrack` - A finite element solver for simulating crack propagation based on linear elastic fracture mechanics. URL: https://doi.org/10.5281/zenodo.18247511, doi:10.5281/zenodo.18247512.

Lorentz, E., Badel, P., 2004. A new path-following constraint for strain-softening finite element simulations. International Journal for Numerical Methods in Engineering 60, 499–526. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.971, doi:10.1002/nme.971. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.971.

Ma, R., Sun, W., Guo, T., 2024. Phase-field model for ductile fracture in the stress resultant geometrically exact shell. International Journal for Numerical Methods in Engineering URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-85187494595&doi=10.1002%2fonme.7462&partnerID=40&md5=24d01a25490cb3eb092899b024912d1d, doi:10.1002/nme.7462.

Maggiorelli, E., Negri, M., 2025. Energy release and Griffith's criterion for phase-field fracture. URL: http://arxiv.org/abs/2502.10351, doi:10.48550/arXiv.2502.10351. arXiv:2502.10351 [math].

May, I.M., Duan, Y., 1997. A local arc-length procedure for strain softening. Computers & Structures 64, 297–303. URL: https://www.sciencedirect.com/science/article/pii/S0045794996001721, doi:10.1016/S0045-7949(96)00172-1.

May, S., Vignollet, J., de Borst, R., 2016. A new arc-length control method based on the rates of the internal and the dissipated energy. Engineering Computations 33, 100–115. URL: https://doi.org/10.1108/EC-02-2015-0044, doi:10.1108/EC-02-2015-0044. publisher: Emerald Group Publishing Limited.

Melin, S., 1983. Why do cracks avoid each other? International Journal of Fracture 23, 37–45. URL: https://doi.org/10.1007/BF00020156, doi:10.1007/BF00020156.

Miehe, C., Hofacker, M., Welschinger, F., 2010. A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering 199, 2765–2778. URL: https://www.sciencedirect.com/science/article/pii/S0045782510001283, doi:10.1016/j.cma.2010.04.011.

Morano, C., Zavattieri, P., Alfano, M., 2020. Tuning energy dissipation in damage tolerant bio-inspired interfaces. Journal of the Mechanics and Physics of Solids 141, 103965. URL: https://www.sciencedirect.com/science/article/pii/S0022509620302003, doi:10.1016/j.jmps.2020.103965.

Mumford, D., Shah, J., 1989. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics 42, 577–685. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/cpa.3160420503, doi:10.1002/cpa.3160420503. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/cpa.3160420503.

Negri, M., 2003. A finite element approximation of the Griffith's model in fracture mechanics. Numerische Mathematik 95, 653–687. URL: https://doi.org/10.1007/s00211-003-0456-y, doi:10.1007/s00211-003-0456-y.
Negri, M., 2007. Convergence analysis for a smeared crack approach in brittle fracture. Interfaces and Free Boundaries 9, 307-330. URL: https://ems.press/journals/ifb/articles/1478, doi:10.4171/ifb/166.

Nuismer, R.J., 1975. An energy release rate criterion for mixed mode fracture. International Journal of Fracture 11, 245-250. URL: https://doi.org/10.1007/BF00038891, doi:10.1007/BF00038891.

Palaniswamy, K., 1972. Crack Propagation under General In-Plane Loading. phd. California Institute of Technology. URL: https://doi.org/10.7907/M7WW-G964, doi:10.7907/M7WW-G964.

Pham, K.H., Amor, H., Marigo, J.J., Maurini, C., 2011. Gradient Damage Models and Their Use to Approximate Brittle Fracture. International Journal of Damage Mechanics 20, 618-652. URL: http://journals.sagepub.com/doi/10.1177/1056789510386852, doi:10.1177/1056789510386852.

Pohl, T., Ramm, E., Bischoff, M., 2014. Adaptive path following schemes for problems with softening. Finite Elements in Analysis and Design 86, 12-22. URL: https://www.sciencedirect.com/science/article/pii/S0168874X14000225, doi:10.1016/j.finel.2014.02.005.

Raghu, P., Rajagopal, A., Reddy, J.N., 2019. Thermodynamically Consistent Variational Approach for Modeling Brittle Fracture in Thick Plates by a Hybrid Phase Field Model. Journal of Applied Mechanics 87. URL: https://doi.org/10.1115/1.4045236, doi:10.1115/1.4045236.

Ramm, E., 1981. Strategies for Tracing the Nonlinear Response Near Limit Points, in: Wunderlich, W., Stein, E., Bathe, K.J. (Eds.), Nonlinear Finite Element Analysis in Structural Mechanics. Springer, Berlin, Heidelberg, pp. 63-89. doi:10.1007/978-3-642-81589-8_5.

Rastiello, G., Oliveira, H.L., Millard, A., 2022. Path-following methods for unstable structural responses induced by strain softening: a critical review. Comptes Rendus. Mécanique 350, 205-236. URL: https://comptes-rendus.academie-sciences.fr/mecanique/articles/10.5802/crmeca.112/, doi:10.5802/crmeca.112.

Rastiello, G., Riccardi, F., Richard, B., 2019. Discontinuity-scale path-following methods for the embedded discontinuity finite element modeling of failure in solids. Computer Methods in Applied Mechanics and Engineering 349, 431-457. URL: https://www.sciencedirect.com/science/article/pii/S004578251930101X, doi:10.1016/j.cma.2019.02.030.

Riks, E., 1972. The Application of Newton's Method to the Problem of Elastic Stability. Journal of Applied Mechanics 39, 1060-1065. URL: https://doi.org/10.1115/1.3422829, doi:10.1115/1.3422829.

Riks, E., 1979. An incremental approach to the solution of snapping and buckling problems. International Journal of Solids and Structures 15, 529-551. URL: https://www.sciencedirect.com/science/article/pii/0020768379900817, doi:10.1016/0020-7683(79)90081-7.

Roux, S., Vandembroucq, D., Hild, F., 2003. Effective toughness of heterogeneous brittle materials. European Journal of Mechanics - A/Solids 22, 743-749. URL: https://www.sciencedirect.com/science/article/pii/S0997753803000780, doi:10.1016/S0997-7538(03)00078-0.

Russ, J.B., Waisman, H., 2020. A novel topology optimization formulation for enhancing fracture resistance with a single quasi-brittle material. International Journal for Numerical Methods in Engineering 121, 2827-2856. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6334, doi:10.1002/nme.6334. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.6334.

Singh, N., Verhoosel, C.V., de Borst, R., van Brummelen, E.H., 2016. A fracture-controlled path-following technique for phase-field modeling of brittle fracture. Finite Elements in Analysis and Design 113, 14-29. URL: https://www.sciencedirect.com/science/article/pii/S0168874X15001870, doi:10.1016/j.finel.2015.12.005.

Tanné, E., Li, T., Bourdin, B., Marigo, J.J., Maurini, C., 2018. Crack nucleation in variational phase-field models of brittle fracture. Journal of the Mechanics and Physics of Solids 110, 80-99. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509617306543, doi:10.1016/j.jmps.2017.09.006.

Triclot, J., Corre, T., Gravouil, A., Lazarus, V., 2023. Key role of boundary conditions for the 2D modeling of crack propagation in linear elastic Compact Tension tests. Engineering Fracture Mechanics 277, 109012. URL: https://www.sciencedirect.com/science/article/pii/S0013794422007354, doi:10.1016/j.engfracmech.2022.109012.

Triclot, J., Corre, T., Gravouil, A., Lazarus, V., 2024. Toughening effects of out-of-crack-path architected zones. International Journal of Fracture URL: https://doi.org/10.1007/s10704-024-00811-5, doi:10.1007/s10704-024-00811-5.

Vasoya, M., Lazarus, V., Ponson, L., 2016. Bridging micro to macroscale fracture properties in highly heterogeneous brittle solids: weak pinning versus fingering. Journal of the Mechanics and Physics of Solids 95, 755-773. URL: https://www.sciencedirect.com/science/article/pii/S0022509615303604, doi:10.1016/j.jmps.2016.04.022.

Verhoosel, C.V., Remmers, J.J.C., Gutiérrez, M.A., 2009. A dissipation-based arc-length method for robust simulation of brittle and ductile failure. International Journal for Numerical Methods in Engineering 77, 1290-1321. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.2447, doi:10.1002/nme.2447. _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.2447.

Vignollet, J., May, S., de Borst, R., Verhoosel, C.V., 2014. Phase-field models for brittle and cohesive fracture. Meccanica 49, 2587-2601. URL: https://doi.org/10.1007/s11012-013-9862-0, doi:10.1007/s11012-013-9862-0.

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., 2020. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 17, 261-272. URL: https://www.nature.com/articles/s41592-019-0686-2, doi:10.1038/s41592-019-0686-2. publisher: Nature Publishing Group.

Wambacq, J., Ulloa, J., Lombaert, G., François, S., 2021. A dissipation-based path-following technique for the phase-field approach to brittle and ductile fracture. International Journal for Numerical Methods in Engineering 122, 3919-3940. URL: https://onlinelibrary.wiley.com/doi/10.1002/nme.6687, doi:10.1002/nme.6687.

Wempner, G.A., 1971. Discrete approximations related to nonlinear theories of solids. International Journal of Solids and Structures 7, 1581-1599. URL: https://www.sciencedirect.com/science/article/pii/0020768371900382, doi:10.1016/0020-7683(71)90038-2.

Wu, J., 2021. On the theoretical and numerical aspects of the unified phase-field theory for damage and failure in solids and structures. Lixue Xuebao/Chinese Journal of Theoretical and Applied Mechanics 53, 301-329. URL: https://www.scopus.com/inward/record.uri?eid=2-s2.0-85102399545&doi=10.6052/2f0459-1879-20-295&partnerID=40&md5=33a959120eba25972908c13418c9f53d, doi:10.6052/0459-1879-20-295.

Wu, J.Y., 2018. Robust numerical implementation of non-standard phase-field damage models for failure in solids. Computer Methods in Applied
Mechanics and Engineering 340, 767-797. URL: https://www.sciencedirect.com/science/article/pii/S0045782518303050, doi:10.1016/j.cma.2018.06.007.

Wu, Z., Guo, L., Hong, J., 2022. Improved Staggered Algorithm for Phase-Field Brittle Fracture with the Local Arc-Length Method. Computer Modeling in Engineering and Sciences 134, 611-636. URL: https://www.techscience.com/CMES/v134n1/49446, doi:10.32604/cmes.2022.020694. publisher: Tech Science Press.

Xia, S.M., Ponson, L., Ravichandran, G., Bhattacharya, K., 2013. Adhesion of heterogeneous thin films—I: Elastic heterogeneity. Journal of the Mechanics and Physics of Solids 61, 838-851. URL: https://www.sciencedirect.com/science/article/pii/S0022509612002359, doi:10.1016/j.jmps.2012.10.014.

Yang, S., Shen, Y., 2023. An acceleration scheme for the phase field fatigue fracture simulation with a concurrent temporal homogenization method. Computer Methods in Applied Mechanics and Engineering 416, 116294. URL: https://www.sciencedirect.com/science/article/pii/S0045782523004188, doi:10.1016/j.cma.2023.116294.

Yu, H., Kuna, M., 2021. Interaction integral method for computation of crack parameters K-T – A review. Engineering Fracture Mechanics 249, 107722. URL: https://www.sciencedirect.com/science/article/pii/S0013794421001739, doi:10.1016/j.engfracmech.2021.107722.

Zambrano, J., Toro, S., Sánchez, P.J., Duda, F.P., Méndez, C.G., Huespe, A.E., 2023. An arc-length control technique for solving quasi-static fracture problems with phase field models and a staggered scheme. Computational Mechanics URL: https://doi.org/10.1007/s00466-023-02388-7, doi:10.1007/s00466-023-02388-7.