# Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture

Andrés A León Baldelli, Corrado Maurini

► To cite this version:

Andrés A León Baldelli, Corrado Maurini. Numerical bifurcation and stability analysis of variational gradient-damage models for phase-field fracture. Journal of the Mechanics and Physics of Solids, 2021, Volume 152, 152 (104424), 10.1016/j.jmps.2021.104424. hal-03106368v2

HAL Id: hal-03106368

https://hal.science/hal-03106368v2

Submitted on 2 Apr 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Numerical bifurcation and stability analysis of variational
gradient-damage models for phase-field fracture

Andrés A. León Baldelli¹, Corrado Maurini²

¹IMSIA, CNRS, EDF, CEA, ENSTA Paris, Institut Polytechnique de Paris, UMR
9219, Palaiseau 91732, France

²CNRS, Institut Jean Le Rond d'Alembert, Sorbonne University, UMR 7190, 75005, Paris, France, corrado.maurini@sorbonne-universite.fr

April 2, 2021

Abstract

Gradient damage models used in phase-field approaches to brittle fracture are charac-
terised by material softening and instabilities. We present novel numerical techniques for
the bifurcation and stability analysis along quasi-static evolution paths as well as practical
tools to select stable evolutions. Our approach stems from the variational approach to
fracture and the theory of rate-independent irreversible processes whereby a quasi-static
evolution is formulated in terms of incremental energy minimisation under unilateral con-
straints. Focussing on the discrete setting obtained with finite elements techniques, we
discuss the links between bifurcation criteria for an evolution and stability of equilibrium
states. Key concepts are presented through the analytical solution of a two-degrees-
of-freedom model featuring a continuum family of bifurcation branches. We introduce
numerical methods to i) assess (second-order) stability conditions for time-discrete evo-
lutions subject to damage irreversibility, and ii) to select possible stable evolutions based
on an energetic criterion. Our approach is based on the solution of a coupled eigen-
value problem which accounts for the time-discrete irreversibility constraint on damage.
Several numerical examples illustrate that this approach allows us to filter out unstable
solutions provided by standard (first-order) minimisation algorithms as well as to effec-
tively compute stable evolution paths. We demonstrate our purpose on a multifissuration
problem featuring complex fracture patterns, to show how the eigenvalue analysis enables
to compute and retrieve morphological properties of emerging cracks.

# 1 Introduction

Fracture and damage are dissipative irreversible phenomena. Their models describe the cre-
ation of cracks in the form of sharp interfaces or the evolution damage as a diffuse reduction
of the material stiffness, associating these processes to a specific energy dissipation per unit
surface or per unit volume, respectively. Modern variational approaches formulate damage
and crack evolution problems as energy minimisation problems in the framework of the the-
ory of quasi-static rate-independent processes (Mielke and Roubíček, 2015; Pham and Marigo,
2010). They led to the establishment of a precise energetic link between a special class of
softening gradient-damage models and brittle fracture theories (Bourdin et al., 2000; Franc-
fort and Marigo, 1998; Marigo et al., 2016; Pham et al., 2011b). This special class of gradient
damage models, commonly referred to as phase-field fracture models, are nowadays among
the most flexible and robust tools for the numerical simulation of brittle fracture phenomena.

Quasi-static evolution problems arising in fracture and the associated softening damage
models are strongly nonlinear. They can admit multiple solutions, or none. In this frame-
work, a stability and bifurcation analysis is important to select among the several possible
solutions, or evolution paths, the physically relevant ones (Bigoni, 2012). For irreversible
and non-smooth processes like damage, fracture, or plasticity, this analysis is subtler than in
the conservative elastic case. Similar problems of nonlinear continuum mechanics have been
largely studied in the past, leading to the formulation of specific criteria for non-bifurcation
and stability of solutions. The classical example of plastic buckling of a two-degrees-of-
freedom plastic structure, alias the ‘Shanley column’ (Shanley, 1947), illustrates the main
counter-intuitive phenomena that are possible in this setting. One can observe continuous
families of bifurcation points, the presence of special singularities, and the distinction be-
tween the loading threshold for the first primary bifurcation and that leading to the loss of
stability (Bazant and Tabbara, 1992; Benallal et al., 1993; Fedelich and Ehrlacher, 1997; Hill,
1957; Hutchinson, 1974; Nguyen, 1987; Nguyen, 2000; Petryk, 1993; Potier-Ferry, 1985). In
these problems, the non-smoothness renders the incremental problem intrinsically nonlinear.
A common solution for establishing sufficient conditions for uniqueness is to resort to the
analysis of the so-called linear comparison solid (Hill, 1957). Classical analyses based on local
continuum theories consider pointwise conditions for uniqueness and stability, like the loss of
ellipticity of the linear comparison solid or the loss of the complementing condition on the
boundary and material interfaces (Benallal et al., 1993). In gradient theories, like gradient-
damage or gradient-plasticity, the non-local nature of the incremental problem calls for global
(in space) uniqueness and stability conditions. The energetic theory of rate-independent pro-
cesses (see Mielke and Roubíček, 2015) provides a natural framework for their formulation.
Several recent works have performed analytical stability and bifurcation analyses of gradient
damage models in the one-dimensional (Benallal and Marigo, 2007; Pham et al., 2011a) or
the three-dimensional context (Pham and Marigo, 2012), providing theoretical tools for the
understanding of the nucleation of cracks (Tanné et al., 2018) or the morphogenesis of com-
plex patterns (Bourdin et al., 2014; Sicsic et al., 2014) in phase-field fracture. Numerical
techniques for stability and bifurcations analysis have been proposed for plasticity (Petryk
and Thermann, 1992). Yet, at the best of our knowledge, similar numerical tools are not
available for the numerical solution of gradient damage and phase-field fracture models.

The goal of this paper is to present a numerical method for the bifurcation and stability
analysis of gradient damage models used in phase-field fracture. We illustrate on several test
cases how it can improve the prediction of standard algorithms for the numerical solution
of the evolution problem. Differently from previous attempts (Beaurain, 2011; Beaurain et
al., 2011), we focus on the bifurcation and stability analysis in the time-discrete framework,
where the evolution problem is formulated as an incremental energy minimisation under an
irreversibility constraint on damage, bounded to be pointwise larger than its value at the
previous time-step. Hence, we apply standard methods of the theory of bound-constrained
optimisation (see e.g. Nocedal and Wright, 2006) to determine non-bifurcation and stability
criteria. In addition to classical algorithms used to solve phase-field fracture problems, we
explicitly test the second-order optimality conditions for the time-discrete energy minimisation
problem by solving a suitable eigenvalue problem on a reduced Hessian matrix. Hence, we
exploit the available information to improve the time-continuation of the numerical solution
of the evolution problem. We use standard finite element techniques for spatial discretisation,
but the approach is applicable to other methods.
This paper is organized as follows. In Section 2, we formulate the evolution problem at hand and give our definitions of stability in the time-continuous and time-discrete setting. In Section 3, we treat in detail a simple model problem of a two-degree-of-freedom system composed by two softening springs connected in series. This example allows us to highlight the richness of possible phenomena that can be encountered even in such a simple setting where analytical solutions are available. Section 4 presents the numerical algorithms used for the bifurcation and stability analysis. Section 5 is devoted to the illustration of a series of numerical experiments highlighting the main properties of stability and bifurcation in two representative systems, namely the one dimensional bar in traction and a thin film bonded onto a stiff substrate. Conclusions and perspectives are drawn in Section 6.

We briefly introduce the notation conventions adopted in this paper. We will indifferently use a lower case letter for scalar-, vector-, or tensor-valued fields depending on the space variable $x$. For example, $u: x \in \mathbb{R}^d \to u(x) \in \mathbb{R}^d$ and $\alpha: x \in \mathbb{R}^d \to \alpha(x) \in \mathbb{R}$ will denote the displacement and damage fields, $d$ being the dimension of the physical space. A superposed dot will denote the time derivative, e.g. $\dot{u}_t(x) = \frac{d}{dt}u_t(x, t)$, whilst a prime sign will denote the derivative of a function with respect to its argument, e.g. $w'(\alpha) = \frac{dw(\alpha)}{d\alpha}$. We use boldface fonts for $n$-dimensional vectors of finite element coefficients for spatial fields, e.g. $\mathbf{u}$ and $\alpha$ will denote the finite element coefficients of $u$ and $\alpha$. Matrices will be denoted by capital letters, using a boldface font for finite element stiffness matrices. For example, a linear system obtained after a finite element discretisation will read as $\mathbf{Au} = \mathbf{b}$. As usual, $H^1(\Omega, \mathbb{R}^d)$ denotes the Sobolev space of functions with values in $\mathbb{R}^d$ which are square integrable with square integrable first derivatives.

## 2 Formulation of the time-evolution problem and stability

We present here the variational formulation of the time evolution problem for gradient damage models in the framework of quasi-static rate-independent processes and generalised standard materials, following (Halphen and Nguyen, 1975; Mielke and Roubíček, 2015; Pham and Marigo, 2010). We focus mainly on the time-discrete case, which is the setting adopted for the numerical analysis presented in this paper. To avoid the technicalities related to the infinite-dimensional nature of the continuum damage model, we quickly simplify the presentation introducing the finite-dimensional space-discretization of the continuum problem. We define the criteria for the stability and uniqueness of solutions, and illustrate them through the analytical solution of a two-degree-of-freedom model problem.

### 2.1 The total energy functional of gradient damage models

Let be $u$ the displacement field (possibly vector-valued), and $\alpha$ a scalar damage field defined on the domain $\Omega$. We will denote by $e(u) = (\nabla u + \nabla^T u)/2$ the linearised strain field. We consider the classical model of isotropic gradient damage in the small strain limit used in the phase-field Ambrosio-Tortorelli (Ambrosio and Tortorelli, 1992) approximation of variational brittle fracture (Francfort and Marigo, 1998). The bulk energy stored in the volume element is of the form

$$W(e, \alpha, \nabla \alpha) = \frac{1}{2} a(\alpha) A e(u) \cdot e(u) + w(\alpha) + w_1 \ell^2 \nabla \alpha \cdot \nabla \alpha, \quad (1)$$

where $a(\alpha)$ and $w(\alpha)$ are two scalar functions representing the stiffness modulation and the energy dissipation along an homogeneous damage process, respectively. We assume that $a$
is a monotonically decreasing function from $a(0) = 1$ to $a(1) = \eta_\ell \ll \ell/L \ll 1$ and $w$ is a monotonically increasing function from $w(0) = 0$ to $0 < w(1) = w_1 < +\infty$, where we denote by $L$ the typical size of the domain $\Omega$. We consider the case of isotropic elasticity where the fourth-order elasticity tensor is such that $A e = \lambda \operatorname{tr}(e) I + 2\mu e$, $\lambda$ and $\mu$ being the Lamé parameters and $I$ the identity tensor. The constant $w_1$ is the specific fracture energy and $\ell$ is an internal length of the material. As shown in (Pham et al., 2011b; Tanné et al., 2018), for brittle materials with stress softening and mode-I cracks, these parameters can be determined from the knowledge of the toughness $G_c$ and the critical stress $\sigma_c$. The analysis of the solutions in the form of damage localisations, regarded as smeared versions of sharp cracks, gives the following relation between the fracture toughness $G_c$ and the specific fracture energy $w_1$

$$G_c = 4c_w w_1 \ell, \quad \text{where} \quad c_w = \int_0^1 \sqrt{w(\alpha)} d\alpha.$$

We consider a quasi-static rate-independent evolution parametrised by a “time” variable $t$ and introduce the total energy functional at time $t$

$$\mathcal{E}_t(u, \alpha) := \int_{\Omega} W(e(u(x)), \alpha(x), \nabla\alpha(x)) \, dx - \mathcal{F}_t(u(x)), \quad (2)$$

where $\mathcal{F}_t$ is a linear functional of the displacement representing the work of conservative external body forces and surface tractions applied at time $t$. We denote by

$$\mathcal{V}_t := \{v \in H^1(\Omega, \mathbb{R}^d), v = \bar{u}_t \text{ on } \partial_u \Omega\}, \\ \mathcal{D}(\alpha_-) := \{\alpha \in H^1(\Omega) : \alpha = 0 \text{ on } \partial_\alpha \Omega, \alpha \ge \alpha_-\},$$

the space of admissible displacements (compatible with non-homogeneous Dirichlet boundary conditions $\bar{u}_t$) and admissible damage fields (from a damaged state $\alpha_-$). Boundary conditions are applied on the subsets $\partial_u \Omega$ and $\partial_\alpha \Omega$ of the boundary $\partial \Omega$. To simplify the presentation, we assume that the Dirichlet boundary conditions on the damage field, if present, are homogeneous. Non-homogeneous boundary conditions, often used in practice to prescribe, e.g., the presence of an initial crack, can be accounted for without any complication. Assuming an existing state of damage $\alpha_-(x) > 0$, taking $(u, \alpha)$ in $H^1(\Omega, \mathbb{R}^d) \times \mathcal{D}(\alpha_-)$ guarantees a finite value of the energy (2) and the irreversibility of damage with respect to the previous damage level $\alpha_-$. Also, the irreversibility implies that $w(\alpha)$ represents a dissipative contribution to the internal energy (1). The formulation and the basic ideas presented in the rest of this paper can be extended without any major changes to more complex material models simply by replacing the definition of the strain energy $W$. One can include, for example, traction-compression asymmetry (Amor et al., 2009; Freddi and Royer-Carfagni, 2010), anisotropy (Bleyer and Alessi, 2018; Li et al., 2015), or large deformations (Del Piero et al., 2007).

## 2.2 Finite-element discretisation in space

In the rest of this paper, we focus on the analysis of the space-discrete problem with $n_u$ displacement degrees-of-freedom $\mathbf{u} = \{u_1, \dots, u_{n_u}\}^T$ and $n_\alpha$ damage degrees-of-freedom $\alpha = \{\alpha_1, \dots, \alpha_{n_\alpha}\}^T$. A similar discrete problem is obtained from (2) after a finite element dis-
cretisation of the displacement and damage field $(u, \alpha) \in \mathcal{V}_t \times \mathcal{D}(\alpha_-)$ in the form

$$u(x) \approx u_t^{(\text{FE})}(x) := \sum_{i=1}^{n_u} \phi_h^{(u)}(x) u_h + \sum_{i=1}^{\bar{n}_u} \phi_{n_u+h}^{(u)}(x) \bar{u}_{t,h}, \quad (3)$$

$$\alpha(x) \approx \alpha_t^{(\text{FE})}(x) := \sum_{i=1}^{n_\alpha} \phi_h^{(\alpha)}(x) u_h + \sum_{i=1}^{\bar{n}_\alpha} \phi_{n_\alpha+h}^{(\alpha)}(x) \bar{\alpha}_{t,h}, \quad (4)$$

where $\phi_h^{(u)} \in H_0^1(\Omega, \mathbb{R}^d)$ and $\phi_h^{(\alpha)} \in H_0^1(\Omega)$ are the finite element basis functions for the displacement and damage fields respectively; the known coefficients $\{\bar{u}_{t,h}\}_{h=1}^{\bar{n}_u}, \{\bar{\alpha}_{t,h}\}_{h=1}^{\bar{n}_\alpha}$ are calculated to satisfy the possible $\bar{n}_u + \bar{n}_\alpha$ Dirichlet boundary conditions in the finite element approximation. We denote by $\mathbf{y} := (\mathbf{u}, \mathbf{\alpha}) \in \mathbb{R}^n$, with $n = n_u+n_\alpha$, the global vector describing the unknowns of the discrete system.

We consider only the case of finite element discretisation adopting $P_1$ finite elements with linear Lagrange polynomials as basis functions for the damage variable, for which the point-wise unilateral constraint $\alpha \ge \alpha_-$ translates in the component-wise unilateral constraint on the finite element coefficients $\alpha \ge \alpha_-$. Hence, we denote by

$$\mathcal{C}(\alpha_-) = \{\mathbf{y} := (\mathbf{u}, \mathbf{\alpha}) \in \mathbb{R}^n \equiv \mathbb{R}^{n_u} \times \mathbb{R}^{n_\alpha} : \mathbf{\alpha} \ge \alpha_-\}$$

the convex cone of admissible states from the state with damage $\alpha_-$ and by

$$\mathcal{C}_+ := \mathcal{C}(\mathbf{0}) = \{\mathbf{z} := (\mathbf{v}, \beta) \in \mathbb{R}^n \equiv \mathbb{R}^{n_u} \times \mathbb{R}^{n_\alpha} : \beta \ge 0\}$$

the cone of admissible variations such that, for any $\mathbf{y} \in \mathcal{C}(\alpha_-)$, $\mathbf{y} + \mathbf{z} \in \mathcal{C}(\alpha_-)$ for all $\mathbf{z} \in \mathcal{C}_+$.

The energy of the system in terms of the finite element coefficients $\mathbf{u} \in \mathbb{R}^{n_u}, \mathbf{\alpha} \in \mathbb{R}^{n_\alpha}$ is defined by

$$E_t(\mathbf{y}) = E_t(\mathbf{u}, \mathbf{\alpha}) := \mathcal{E}_t(u_t^{(\text{FE})}, \alpha_t^{(\text{FE})}). \quad (5)$$

We assume the energy to be sufficiently smooth to be expanded with the following second-order Taylor series

$$E_t(\mathbf{y}_t + \epsilon \mathbf{z}) - E_t(\mathbf{y}_t) = \epsilon E'_t(\mathbf{y}_t)(\mathbf{z}) + \frac{\epsilon^2}{2} E''_t(\mathbf{y}_t)(\mathbf{z}) + o(\epsilon^2), \quad (6)$$

where

$$E'_t(\mathbf{y}_t)(\mathbf{z}) := \left. \frac{d}{d\epsilon} E_t(\mathbf{y}_t + \epsilon \mathbf{z}) \right|_{\epsilon=0} = \mathbf{g}_t(\mathbf{y}_t) \cdot \mathbf{z}, \quad \mathbf{g}_t(\mathbf{y}_t) := \left[ \frac{\partial E_t(\mathbf{y}_t)}{\partial \mathbf{u}} \quad \frac{\partial E_t(\mathbf{y}_t)}{\partial \mathbf{\alpha}} \right],$$

$$E''_t(\mathbf{y}_t)(\mathbf{z}) := \left. \frac{d^2}{d\epsilon^2} E_t(\mathbf{y}_t + \epsilon \mathbf{z}) \right|_{\epsilon=0} = \mathbf{H}_t(\mathbf{y}_t) \mathbf{z} \cdot \mathbf{z}, \quad \mathbf{H}_t(\mathbf{y}_t) := \left[ \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u}^2} \quad \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u} \partial \mathbf{\alpha}} \quad \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{\alpha} \partial \mathbf{u}} \right].$$

## 2.3 Time-continuous evolution problem

Following the modern energetic approach to quasi-static rate-independent evolution processes (Mielke and Roubíček, 2015), a solution of the quasi-static process in the time interval $t \in (0, T)$ is defined as a one-parameter family of states $\{(\mathbf{u}_t, \mathbf{\alpha}_t)\}_{t \in [0, T]}$ satisfying three principles of *irreversibility*, *stability*, and *energy balance*. We present them under the strong hypothesis that the solution is smooth in time so that the following time-derivatives

$$\dot{\mathbf{y}}_t = (\dot{\mathbf{u}}_t, \dot{\mathbf{\alpha}}_t) := \lim_{\Delta t \to 0^+} \frac{\mathbf{y}_{t+\Delta t} - \mathbf{y}_t}{\Delta t}, \quad \dot{E}_t(\mathbf{y}_t) := \lim_{\Delta t \to 0^+} \frac{E_{t+\Delta t}(\mathbf{y}_t) - E_t(\mathbf{y}_t)}{\Delta t}$$

are well-defined. For each $t \in (0, T)$, the three requirements read as follows
1. **Time-continuous irreversibility.** The damage field $t \mapsto \alpha_t$ is non-decreasing in time

$$ \dot{\alpha}_t \geq 0. \qquad (7) $$

2. **Directional state-stability.** The solution is a local minimiser of the energy

$$ y_t = (u_t, \alpha_t) \in \arg\min_{y \in \mathcal{C}(\alpha_t)} E_t(y), \qquad (8) $$

which means that $y_t \in \mathcal{C}(\alpha_t)$ and

$$ \forall \mathbf{z} \in \mathcal{C}_+, \quad \exists \bar{\epsilon} > 0 : \quad \forall \epsilon \in [0, \bar{\epsilon}], \quad E_t(\mathbf{y}_t + \epsilon \mathbf{z}) - E_t(\mathbf{y}_t) \geq 0. \qquad (9) $$

3. **Energy balance.** The power exerted by the internal forces should be equal to the power exerted by the external forces. Because of the mechanical equilibrium, this condition is tantamount to requiring that the power expended by the internal generalised forces associated to the damage variable vanishes (Mielke and Roubíček, 2015; Pham and Marigo, 2010):

$$ \frac{\partial E_t}{\partial \alpha}(y_t) \cdot \dot{\alpha}_t = 0. \qquad (10) $$

**Equilibrium paths.** Retaining only the first-order term in the series expansion (6), one can show that the stability condition (8) requires the following first-order local minimality (or optimality) conditions

$$ E'_t(\mathbf{y}_t)(\mathbf{z}) = \frac{\partial E_t(\mathbf{y}_t)}{\partial \mathbf{u}} \cdot \mathbf{v} + \frac{\partial E_t(\mathbf{y}_t)}{\partial \boldsymbol{\alpha}} \cdot \boldsymbol{\beta} \geq 0, \quad \forall \mathbf{z} \in \mathcal{C}_+ \Leftrightarrow \frac{\partial E_t}{\partial \mathbf{u}}(\mathbf{y}_t) = \mathbf{0}, \quad \frac{\partial E_t}{\partial \boldsymbol{\alpha}}(\mathbf{y}_t) \geq \mathbf{0}, \qquad (11) $$

which, together with the irreversibility (7) and the energy balance (10), form the following system of variational inequalities that must be verified at each instant

$$ \frac{\partial E_t}{\partial \mathbf{u}}(\mathbf{y}_t) = \mathbf{0}, \quad \frac{\partial E_t}{\partial \boldsymbol{\alpha}}(\mathbf{y}_t) \geq \mathbf{0}, \quad \dot{\boldsymbol{\alpha}}_t \geq \mathbf{0}, \quad \frac{\partial E_t}{\partial \boldsymbol{\alpha}}(\mathbf{y}_t) \cdot \dot{\boldsymbol{\alpha}}_t = 0. \qquad (12) $$

The last complementary condition states that, for each degree-of-freedom, damage can evolve only if the corresponding derivative of the energy vanishes, and that, vice-versa, the derivative can be positive only if the damage variable hits the corresponding irreversibility constraint. In the following, we will denote an *equilibrium path* a family of states $\{\mathbf{y}_t\}_{t \in (t_0, t_1)}$ respecting the system of inequalities above in some time interval $(t_0, t_1)$ and define an *equilibrium* each single $\mathbf{y}_t$.

**Active and inactive constraints.** For a given equilibrium $\mathbf{y}_t$, let us define the subspace of the state space $\mathbb{R}^n$ where the unilateral constraints due to irreversibility are *inactive* at time $t$ as

$$ \mathcal{N}_t(\mathbf{y}_t) := \left\{ \mathbf{z} = (\mathbf{v}, \boldsymbol{\beta}) \in \mathbb{R}^n : E'_t(\mathbf{y}_t)(\mathbf{z}) = 0 \Leftrightarrow \frac{\partial E_t}{\partial \boldsymbol{\alpha}}(\mathbf{y}_t) \cdot \boldsymbol{\beta} = 0 \right\}. \qquad (13) $$

This vector space includes all the displacement degrees-of-freedom and only the damage degrees-of-freedom for which the derivative of the energy vanishes. It determines the index set of inactive constraints as follows

$$ \mathcal{N}_t := \left\{ k \in \{0, n-1\} : \frac{\partial E}{\partial y_k}(\mathbf{y}_t) = 0 \right\}. \qquad (14) $$
The restriction of a vector or a matrix to this index set corresponds to the projection on the vector space of inactive constraints. We define the subspace of the state space where the constraints are active as $A_t(\mathbf{y}_t) := \mathbb{R}^n \setminus N_t(\mathbf{y}_t)$, the associated index set being $A_t := \{0, \dots, n-1\} \setminus N_t$. We further define the convex-cones of positive variations in the active and nonactive subspaces as follows:

$$ \mathcal{N}_{t+}(\mathbf{y}_t) := \mathcal{N}_t(\mathbf{y}_t) \cap C_+, \quad A_{t+}(\mathbf{y}_t) := A_t(\mathbf{y}_t) \cap C_+. \quad (15) $$

**Second-order state-stability conditions.** For variations $\mathbf{z} \in A_{t+}(\mathbf{y}_t)$, $E'_t(\mathbf{y}_t)(\mathbf{z}) > 0$ and the first-order term in the series expansion (6) is sufficient to guarantee that the increment of the energy is strictly positive. This immediately follows from the definition of $A_{t+}(\mathbf{y}_t)$ and the requirement that $\mathbf{y}_t$ must respect the first-order optimality conditions (11). Vice-versa, for $\mathbf{z} \in A_{t+}(\mathbf{y}_t)$ the first-order term vanishes. To assure the local minimality condition (8) one must require the second derivative to be nonnegative for any non-vanishing admissible variation $\mathbf{z} = (\mathbf{v}, \boldsymbol{\beta})$ with $\boldsymbol{\beta} \ge 0$. This gives the following second order sufficient condition for the directional state-stability of an equilibrium $\mathbf{y}_t$

$$ 0 < \rho_+(t) := \min \left\{ \frac{\mathbf{H}_t(\mathbf{y}_t)\mathbf{z} \cdot \mathbf{z}}{\|\mathbf{z}\|^2}, \quad \mathbf{z} \in \mathcal{N}_{t+}(\mathbf{y}_t), \quad \|\mathbf{z}\| \neq 0 \right\}. \quad (16) $$

The loose inequality $\rho_+(t) \ge 0$ is a necessary condition for the directional stability. Hence, we can mark an equilibrium as state-stable (S-stable) if $\rho_+(t) > 0$ and as state-unstable (S-unstable) if $\rho_+(t) < 0$.

**Rate-problem and bifurcation conditions.** Taking the derivative of first-order conditions (12) with respect to $t$ gives the following rate problem to determine the rate $\dot{\mathbf{y}}_t$ at a given equilibrium $\mathbf{y}_t$

$$ \begin{gathered} \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u}^2} \dot{\mathbf{u}}_t + \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u} \partial \boldsymbol{\alpha}} \dot{\boldsymbol{\alpha}}_t + \frac{\partial \dot{E}_t(\mathbf{y}_t)}{\partial \mathbf{u}} = \mathbf{0}, \\ \dot{\boldsymbol{\alpha}}_t \ge \mathbf{0}, \quad \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u} \partial \boldsymbol{\alpha}} \dot{\mathbf{u}}_t + \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \boldsymbol{\alpha}^2} \dot{\boldsymbol{\alpha}}_t + \frac{\partial \dot{E}_t(\mathbf{y}_t)}{\partial \boldsymbol{\alpha}} \ge \mathbf{0}, \\ \dot{\boldsymbol{\alpha}}_t \cdot \left( \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \mathbf{u} \partial \boldsymbol{\alpha}} \dot{\mathbf{u}}_t + \frac{\partial^2 E_t(\mathbf{y}_t)}{\partial \boldsymbol{\alpha}^2} \dot{\boldsymbol{\alpha}}_t + \frac{\partial \dot{E}_t(\mathbf{y}_t)}{\partial \boldsymbol{\alpha}} \right) = \mathbf{0}. \end{gathered} $$

If this problem has a unique solution $\dot{\mathbf{y}}_t$, there is a single possible evolution path through the equilibrium state $\mathbf{y}_t$; vice-versa bifurcations from a fundamental path $\{\mathbf{y}_t\}_{t \in (t_0, t_1)}$ are possible when the rate problem admits multiple solutions. A sufficient condition for the uniqueness of the solution of the rate problem is that the Hessian matrix is positive definite on the vector space $\mathcal{N}(\mathbf{y}_t)$, i.e.:

$$ 0 < \rho(t) := \min \left\{ \frac{\mathbf{H}_t(\mathbf{y}_t)\mathbf{z} \cdot \mathbf{z}}{\|\mathbf{z}\|^2}, \quad \mathbf{z} \in \mathcal{N}_t(\mathbf{y}_t), \quad \|\mathbf{z}\| \neq 0 \right\}. \quad (17) $$

Indeed, when the above holds, the problem at hand is equivalent to the minimisation of a convex quadratic functional over a convex cone (Nocedal and Wright, 2006). Hence, a necessary condition for having a bifurcation at $\mathbf{y}_t$ is that $\rho(t) \le 0$. Introducing the reduced Hessian matrix $\tilde{\mathbf{H}}_t(\mathbf{y}_t)$ as the restriction of the full Hessian $\mathbf{H}_t(\mathbf{y}_t)$ to the degrees of freedom corresponding to inactive constraints and denoting by

$$ \{\lambda_t^{(k)}, \tilde{\mathbf{z}}_t^{(k)}\}_{k=1}^{\dim(\mathbb{N}_i)}, \quad \text{with} \quad \lambda_t^{(k)} \le \lambda_t^{(k+1)} \quad \text{and} \quad \tilde{\mathbf{z}}_t^{(k)} := (\mathbf{v}_t^{(k)}, \tilde{\boldsymbol{\beta}}_t^{(k)}), \quad (18) $$
its ordered pairs of real eigenvalues and eigenvectors, the bifurcation marker $\rho(t)$ coincides with the smallest eigenvalue $\lambda_t^{(1)}$. When $\lambda_t^{(1)} = 0$, a bifurcation is possible. The corresponding bifurcation mode is given by $z_i^{(1)} = (v_i^{(1)}, \beta_i^{(1)})$, where we denote by $z_i^{(k)} = (v_i^{(k)}, \beta_i^{(k)})$ the extension of $\tilde{z}_i^{(k)} = (v_i^{(k)}, \tilde{\beta}_i^{(k)})$ to $\mathbb{R}^n$ obtained by setting the components of $\beta_i^{(k)} = 0$ for all $k \in A_t$. If multiple eigenvalues vanish simultaneously, the possible bifurcation modes will be in the space generated by the corresponding eigenvectors.

## 2.4 Time-discrete evolution problem

The numerical solution of the evolution problem requires a discretisation of the time variable. We consider a finite-difference discretisation in $N$ steps $\{t_i\}_{i=1}^N$ so that the time derivative of a function $f(t)$ is approximated by the elemental Euler expansion $\dot{f}(t_i) \approx \frac{f(t_i)-f(t_{i-1})}{t_i-t_{i-1}}$. Numerical works on quasi-static variational approaches to fracture and phase-field fracture models seek the state at time step $t_i$ as the solution of the following minimisation problem for a discrete energy functional (5), that is

$$ y_i := (u_i, \alpha_i) \in \arg\min_{(u, \alpha) \in \mathcal{C}(\alpha_{i-1})} E_i(u, \alpha) \quad (19) $$

where, here and henceforth, for any function of time $f_t$, $f_i$ stands for evaluation at the discrete time $f_{t_i}$. We state the local minimality requirement (19) in the following directional form

Find $y_i \in \mathcal{C}(\alpha_{i-1})$ such that

$$ \forall y \in \mathcal{C}(\alpha_{i-1}), \quad \exists \epsilon > 0 : \quad E_i(y_i + \epsilon(y - y_i)) - E_i(y_i) \ge 0 \quad \forall \epsilon \in [0, \epsilon]. \quad (20) $$

**Time-discrete equilibrium paths.** Retaining only the first-order term of the Taylor expansion (6) for the energy increment and denoting by $z = (v, \beta)$ the variation $\epsilon(y - y_i)$, (20) implies that:

$$ E'_i(y_i)(z) = \frac{\partial E_i}{\partial u}(y_i) \cdot v + \frac{\partial E_i}{\partial \alpha}(y_i) \cdot \beta \ge 0, \quad \forall z = (v, \beta) : v \in \mathbb{R}^{n_u}, \beta \in \mathbb{R}^{n_\alpha}, \beta \ge \alpha_{i-1} - \alpha_i. \quad (21) $$

Differently from (11), admissible variations are not only in the form $=(\pm v, \beta)$ with $\beta \ge 0$. Variations in the form $(\pm v, \pm \beta)$ with $\|\beta\|$ sufficiently small are admissible provided that the damage variation is null wherever the damage does not evolve, i.e. $\beta \cdot (\alpha_{i-1} - \alpha_i) = 0$. Selecting a similar class of variations in the optimality condition (21) shows that, in addition to (11), the derivative of the energy with respect to the damage should vanish in the subspace where the damage evolves, i.e. $\frac{\partial E_i}{\partial \alpha}(y_i) \cdot (\alpha_i - \alpha_{i-1}) = 0$. Hence, we have the following proposition, which states the classical Karush-Kuhn-Tucker optimality conditions for a bound-constrained minimisation problem. We refer the reader to classical text on convex optimisation for a detailed proof, see e.g. (Nocedal and Wright, 2006).

**Proposition 1 (First-order optimality conditions).** A solution $y_i = (u_i, \alpha_i) \in \mathbb{R}^n$ of the local minimisation problem (19) must respect the following conditions

$$ \frac{\partial E_i}{\partial u}(y_i) = 0, \quad (22a) $$

$$ \alpha_i - \alpha_{i-1} \ge 0, \quad \frac{\partial E_i}{\partial \alpha}(y_i) \ge 0, \quad \frac{\partial E_i}{\partial \alpha}(y_i) \cdot (\alpha_i - \alpha_{i-1}) = 0. \quad (22b) $$
Conditions (22) can be interpreted as the discrete version of the equilibrium equation and the damage criterion, respectively. The set of inequalities (22b) states that damage cannot decrease (irreversibility), that the first derivative of the energy is nonnegative, and that, for each component of the vector, either the damage does not evolve or the first derivative vanishes. An equilibrium path in the time-discrete setting is defined as a sequence of equilibria ${\{y_i\}}_{i=0}^N$ verifying (22) at each time step.

**Incremental stability.** Given an equilibrium $y_i$ and the corresponding spaces of inactive and active constraints $N_i(y_i)$ and $A_i(y_i)$, any variation $z$ can be decomposed in its orthogonal projections $z_N \in N_i(y_i)$ and $z_A \in A_i(y_i)$ such that $z = z_N + z_A$. Because of the definition of these spaces, $E'_i(y_i)(z_N) = 0$ and $E'_i(y_i)(z_A) > 0$. Hence, the second-order Taylor expansion of the energy increment in local minimality conditions (20) reads as

$$E_i(y_i + \epsilon z) - E_i(y_i) = \epsilon E'_i(y_i)(z_A) + \frac{\epsilon^2}{2} (\mathbf{H}_i(y_i) \mathbf{z}_A \cdot \mathbf{z}_A + \mathbf{H}_i(y_i) \mathbf{z}_N \cdot \mathbf{z}_N + 2\mathbf{H}_i(y_i) \mathbf{z}_A \cdot \mathbf{z}_N) + o(\epsilon^2).$$

If $z_A \neq 0$, the first-order term is strictly positive hence the variation is positive for $\epsilon$ sufficiently small. If $z_A = 0$, then the sign of the variation is given by the sign of the quadratic term in $z_N$. Hence, we can conclude that the local minimality condition (20) is verified if first order conditions (22) are verified and the second derivative is (strictly) positive for all $z_N \in N_i(y_i)$. This result can be summarised in the following set of second order conditions for the solution of the unilaterally constrained minimisation problem (19), see Nocedal and Wright, 2006.

**Proposition 2 (Second-order conditions).** An admissible state $y_i = (\mathbf{u}_i, \alpha_i) \in \mathbb{R}^n$ is a solution of the local minimisation problem (19) if first-order conditions (22) are satisfied and if the second derivative is positive definite on the subspace of inactive constraints $N_i(y_i)$, that is

$$0 < \rho_i = \rho(t_i) := \min \left\{ \frac{\mathbf{H}_i(y_i)\mathbf{z} \cdot \mathbf{z}}{\|\mathbf{z}\|^2}, \ \mathbf{z} \in N_i(y_i), \ \|\mathbf{z}\| \neq 0 \right\}. \quad (23)$$

Conversely, a necessary condition for $y_i$ to be a minimum is the positive semi-definiteness of the quadratic form $\mathbf{H}_i$, i.e. the loose inequality $\rho_i \ge 0$.

In the following, we will denote these second-order conditions for local minimality in the time-discrete setting as *incremental stability* criteria for an equilibrium. We will call solutions $y_i$ of (22) an *incremental equilibrium*, and we will denote them as *incrementally stable* (I-Stable) when $\rho_i > 0$.

## 2.5 Comments

The comparison of the notions of evolution, stability, and bifurcation in the time-continuous and time-discrete evolution calls for several remarks.

* The system of inequalities given by the time-continuous first-order optimality conditions (11) and the energy balance (10) is equivalent to the system of time-discrete first-order optimality conditions (22) when replacing $\dot{\alpha}_t$ with its finite difference approximation $(\alpha_i - \alpha_{i-1})/(t_i - t_{i-1})$. Indeed, for smooth-in-time evolutions, the solution at each time step of the time-discrete first-order conditions (22) is a solution of the time-continuous first-order optimality condition and energy balance.
* The non-bifurcation criterion (17) for the time-continuous rate problem coincides with the incremental stability condition (23) of the time-discrete setting.

* Since $N_{t+}(y_i) \subset N_t(y_i)$,

$$ \rho_i = \rho(t_i) \le \rho_+(t_i) \qquad (24) $$

i.e. the incremental stability (or non-bifurcation) criterion, in the sense of (23)-(17), is a sufficient (but not necessary) condition for the state-stability (16).

* Assuming all the required time-regularity so that the time-discrete evolution has a well-defined smooth time-continuous limit as the time step goes to zero, the time-continuous limit will be a solution of the time-continuous evolution problem, respecting the time-continuous irreversibility, stability, and energy balance. The converse is not true, because a state-stable solution of the time-continuous evolution problem can be incrementally unstable. The example of the next section will provide an explicit illustration of this case.

# 3 A discrete model problem

## 3.1 Formulation

We discuss the solution of the evolution problem of a simple discrete system constituted by two identical damageable springs connected in series, see the schematic in Figure 1. Both springs, supposed undamaged at the beginning of the loading process, can undergo damage in the sense that their stiffness decreases when the corresponding internal damage variable increases. This is a minimal model problem, introduced in Lazarus et al., 2015, which unveils several fundamental issues appearing in the evolution problem including softening elements. We

<figure><img src="image_2.png" /><figcaption>Figure 1: The model 2-dof system composed of two identical damaging springs in series. The evolution problem associated to a monotonic traction experiment allows us to compute, in closed form, bifurcation and stability limits for homogeneous ($\alpha_1 = \alpha_2$) as well as all possible non-homogeneous ($\alpha_1 \neq \alpha_2$) evolution paths.</figcaption></figure>

assume that the system is fixed at one end and that a monotonically increasing displacement $\bar{u}_t = t$ is prescribed at the other end by a hard loading device. The parameter $t$ plays here the role of the time-like variable of rate-independent process. We denote by $u$ the unknown displacement of the point connecting the two springs at time $t$ and by $\alpha = (\alpha_1, \alpha_2) \in [0, 1] \times [0, 1]$ the damage level of each of the two springs. For this discrete structure, the equivalent of the total energy functional introduced in Section 2.1 is

$$ E_t(u, \alpha) = \frac{a_0}{2} (a(\alpha_1)u^2 + a(\alpha_2)(u-t)^2) + w_1(w(\alpha_1) + w(\alpha_2)) $$

where $a_0$ is the extensional stiffness of the undamaged spring and $w_1$ is the dissipated energy in the completely damaged state. We omit here the discrete equivalent of the gradient term in
(1). In order to make the computations explicit, we specify the constitutive model by setting

$$a(\alpha) := \frac{(1 - \alpha)}{(k - 1)\alpha + 1}, \quad w(\alpha) := \alpha, \qquad (25)$$

where $k > 1$ is a free parameter. The chosen one-parameter family of damage models is constructed such that the global one-dimensional homogeneous response in a traction experiment, as it will be seen in the following, gives raise to a piecewise-linear stress-strain curve during both the elastic and softening phase.

The first order optimality conditions, given by (22) in the time discrete setting and by (11) in the time continuous case, require that the first derivative of the energy with respect to $u$ vanishes at each $t$. This condition corresponds to mechanical equilibrium and allows to determine the mid-point displacement $u$ as a function of the damage variables and imposed load. We compute

$$\frac{\partial E_t}{\partial u} = a_0 a(\alpha_1)u - a_0 a(\alpha_2)(t - u) = 0 \quad \Rightarrow \quad u = \frac{a(\alpha_2)}{a(\alpha_1) + a(\alpha_2)}t.$$

We can rewrite the energy in the terms of $(\alpha_1, \alpha_2)$ only

$$E_t(\alpha) := \frac{a_0 t^2}{2} \frac{1}{s(\alpha_1) + s(\alpha_2)} + w_1(w(\alpha_1) + w(\alpha_2)),$$

where $s(\alpha) := 1/a(\alpha)$ is the compliance modulation function.

The solutions to the evolution problem must satisfy, at each time, the set of Karush-Kuhn-Tucker conditions given by (7), (10), and (11) in the time-continuous case, namely

$$\dot{\alpha}_1 \ge 0, \quad \frac{\partial E_t}{\partial \alpha_1} \ge 0, \quad \frac{\partial E_t}{\partial \alpha_1} \dot{\alpha}_1 = 0, \qquad (26a)$$

$$\dot{\alpha}_2 \ge 0, \quad \frac{\partial E_t}{\partial \alpha_2} \ge 0, \quad \frac{\partial E_t}{\partial \alpha_2} \dot{\alpha}_2 = 0, \qquad (26b)$$

with

$$\frac{\partial E_t}{\partial \alpha_1} = -\frac{a_0 t^2}{2} \frac{s'(\alpha_1)}{(s(\alpha_1) + s(\alpha_2))^2} + w_1 w'(\alpha_1) = \frac{s'(\alpha_1)}{2a_0} (\bar{\sigma}^2(\alpha_1) - \sigma_t^2(\alpha_1, \alpha_2)),$$

$$\frac{\partial E_t}{\partial \alpha_2} = -\frac{a_0 t^2}{2} \frac{s'(\alpha_2)}{(s(\alpha_1) + s(\alpha_2))^2} + w_1 w'(\alpha_2) = \frac{s'(\alpha_2)}{2a_0} (\bar{\sigma}^2(\alpha_2) - \sigma_t^2(\alpha_1, \alpha_2)),$$

where

$$\sigma_t(\alpha_1, \alpha_2) := \frac{a_0}{s(\alpha_1) + s(\alpha_2)} t, \quad \bar{\sigma}(\alpha) := \sqrt{\frac{2a_0 w_1 w'(\alpha)}{s'(\alpha)}},$$

define the stress and the maximum allowable stress as a function of the damage variables. For the stress-strain piecewise-linear model (25), we have

$$s(\alpha) = \frac{(k - 1)\alpha + 1}{(1 - \alpha)}, \quad s'(\alpha) = \frac{k}{(1 - \alpha)^2}, \quad w'(\alpha) = 1, \quad \bar{\sigma}(\alpha) = (1 - \alpha)\sigma_c, \quad \sigma_c = \sqrt{\frac{2a_0 w_1}{k}}.$$
## 3.2 Solutions to the evolution problem

Introducing the non-dimensional time variable $\tau$ defined by

$$ \tau = \frac{t}{t_c} \quad \text{with} \quad t_c = \sqrt{\frac{8w_1}{a_0k}}, $$

possible solutions of the evolution problem starting from the undamaged state $\alpha_1(0) = \alpha_2(0) = 0$ at $\tau = 0$, are as follows

* **Purely elastic:** For $\tau \in (0, \tau_c = 1)$, $\partial E_\tau / \partial \alpha_1 > 0$ and $\partial E_\tau / \partial \alpha_2 > 0$. The damage criterion (26) imposes that damage does not evolve, i.e. $\dot{\alpha}_1 = \dot{\alpha}_2 = 0$. The purely elastic solution starting from the undamaged state at $\tau = 0$ is the unique solution to the evolution problem and is such that

$$ \alpha_{1,\tau} = 0, \quad \alpha_{2,\tau} = 0, \quad \sigma_{\tau} = \sigma_c \tau, \quad \sigma_c := \bar{\sigma}(0) = \sqrt{2w_1a_0/k}. \quad (27) $$

At $\tau = \tau_c = 1$ the stress reaches the elastic limit $\sigma_c$ at which $\partial E_\tau / \partial \alpha_1 = \partial E_\tau / \partial \alpha_2 = 0$.

* **With homogeneous damage:** For $\tau > \tau_c$, starting from the elastic solution, one can look for solutions such that damage evolves at the same rate in both springs, i.e. with $\dot{\alpha}_1 = \dot{\alpha}_2 > 0$. The complementarity conditions in (26) give the equations $\partial E_\tau / \partial \alpha_1 = \partial E_\tau / \partial \alpha_2 = 0$. Assuming that both springs are undamaged at $\tau = \tau_c$, one finds

$$ \alpha_{1,\tau} = \alpha_{2,\tau} = \frac{\tau - 1}{k - 1}, \quad \sigma_{\tau} = \sigma_c \frac{k - \tau}{k - 1}. $$

These solutions are valid for $\tau \in (\tau_c, \tau_f := k)$. At $\tau = \tau_f$ the stress vanishes and the damage reaches its upper bound<fnref n="1" /> $\alpha = 1$.

* **With localised damage:** For $\tau \ge \tau_c$, starting from a solution with homogeneous damage $\alpha_1 = \alpha_2 = \alpha^* := (\tau^* - 1)(k - 1)$ at time $\tau^* > \tau_c$, one can look for solutions where only one of the damage variables keeps evolving, i.e. for solutions of (26) with $\dot{\alpha}_1 = 0$, $\partial E_\tau / \partial \alpha_2 = 0$, or with $\dot{\alpha}_2 = 0$, $\partial E_\tau / \partial \alpha_1 = 0$. Modulo a permutation of the indices (1, 2), these solutions are in the form

$$ \alpha_{1,\tau} = \alpha^* := \frac{\tau^* - 1}{k - 1}, \quad \alpha_{2,\tau} = 1 - \left(1 - \frac{\tau}{\tau_s}\right) c^*, \quad \frac{\sigma_\tau}{\sigma_c} = \left(1 - \frac{\tau}{\tau_s}\right) c^*. \quad (28) $$

Here,

$$ c^* = \frac{k(k - \tau^*)}{(k - 1)(k - 2\tau^*)} = \frac{k(1 - \alpha^*)}{(k - 1)(1 - 2\alpha^*) - 1}, \quad \text{and} \quad \tau_s := \frac{k}{2} = \frac{\tau_f}{2}, $$

can be interpreted as the (non-dimensional) rate at which damage increases and the critical load at ultimate failure, respectively.

Figure 2 shows the global stress-strain diagram and the evolution of potential energy corresponding to the three types of solutions above. Because of the irreversibility condition, branches with evolving damage are one-way paths. In Figure 2, the arrows indicate the admissible direction.

Identical solutions can be obtained in a time-discrete setting by imposing, at any discrete time $t_i$, first-order optimality conditions of the time-discrete evolution problem (22).

<footnote n="1">¹To simplify the presentation, we do not explicitly include the upper bound on damage in the formulation and assume a loading $\tau \in (0, \tau_f)$.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 2: Solutions of the discrete model problem in Figure 2. force-displacement diagram (top-left) and evolution of the potential energy (top-right). The solutions is purely elastic for $\tau < \tau_c = 1$ (in blue). For $\tau > \tau_c$ there is a fundamental branch where the damage is evolving equally in the two springs (in orange), showing a continuous family of possible bifurcated branches (in green) where damage increases only in one of the two springs. The homogeneous branch (in orange) is state-stable for $\tau < \tau_s = k/2$ (solid orange line) and state-unstable for $\tau > \tau_s$ (dashed orange line).</figcaption></figure>

## 3.3 Uniqueness and stability

For $\tau \in (0, \tau_c)$ the elastic solution is unique, all constraints are active and the solution is both state-stable and incrementally-stable because the first derivative of the energy is strictly positive.

For $\tau \in (\tau_c, \tau_s)$ the homogeneous solution is state-stable in the sense (16), although fails to be unique and is incrementally unstable in the sense (17)-(23). Indeed, the first derivative of the energy is zero and the stability is assessed by studying the sign of the Hessian matrix

$$ H_{\tau}^{(\text{hom})} = \begin{bmatrix} H_{11}^{(\text{hom})} & H_{12}^{(\text{hom})} \\ H_{12}^{(\text{hom})} & H_{11}^{(\text{hom})} \end{bmatrix}, \\ \text{with } H_{11}^{(\text{hom})} = w_1(k-1)\frac{k-2\tau}{\tau(k-\tau)}, \quad H_{12}^{(\text{hom})} = w_1(k-1)\frac{k}{\tau(k-\tau)}, \quad (29) $$

on $\mathcal{N}(\mathbf{y}_{\tau}) \equiv \mathbb{R}^2$ for the incremental-stability condition (23) and on $\mathcal{N}_+(\mathbf{y}_{\tau}) \equiv \mathbb{R}_+^2$ for the state-stability condition (16). After straightforward calculations, we can show that

$$ \rho^{(\text{hom})}(\tau) = H_{11}^{(\text{hom})} - H_{12}^{(\text{hom})} = -\frac{w_1(k-1)}{(k-\tau)}, \quad \rho_+^{(\text{hom})}(\tau) = H_{11}^{(\text{hom})} = w_1(k-1)\frac{k-2\tau}{\tau(k-\tau)}. \quad (30) $$

where $\rho^{(\text{hom})}$ and $\rho_+^{(\text{hom})}$ are respectively the markers on the homogeneous solution for the incremental-stability and state-stability conditions (23) and (16), and $1 < \tau^* < \tau_f$ is an arbitrary bifurcation load. Hence, the homogeneous solution is state-stable for $\tau \in (1, \tau_s)$ but state-unstable for $\tau \in (\tau_s, \tau_f)$, although it is incrementally-unstable and non unique for any $\tau \in (\tau_c, \tau_f)$. This is consistent with the existence of a continuous family of bifurcated branches (28) departing from the homogeneous solution in the interval $\tau \in (\tau_s, \tau_f)$. Figure 3 provides a graphical interpretation of the stability conditions, reporting the energy landscape in the $\alpha_1-\alpha_2$ plane. For two homogeneous states (marked by white points) obtained for $\tau < \tau_s$ (left) and $\tau > \tau_s$ (right), blue and green are regions with lower or larger energy than the current
<figure><img src="image_4.png" /><figcaption>Figure 3: Energy landscape in the plane $(\alpha_1, \alpha_2)$ for $\tau < \tau_s$ (left) and $\tau > \tau_s$ (right). The current state is denoted by a white point. In blue and green, regions with lower and higher energy than the current state, respectively. Shaded regions are inadmissible because of the irreversibility constraint. For $\tau < \tau_s$, the homogeneous solution is stable because there are no other admissible states with lower energy. For $\tau > \tau_s$ instead, there exist admissible inhomogeneous solutions with decreasing energy, hence the homogeneous state is unstable.</figcaption></figure>

solution, respectively. Shaded regions are inaccessible because of the irreversibility constraint. For $\tau < \tau_s$ the current state is a local minimum of the energy among all the admissible states and the solution is state-stable. For $\tau > \tau_s$, neighbouring states with lower energy are admissible, hence the solution is state-unstable. According to the incremental stability condition (23), the solution with $\tau < \tau_s$ is unstable, because with the irreversibility constraints given by the damage field at the previous time-step, there are admissible neighbouring lower energy states.

On bifurcated branches where only one damage variable can evolve, the incremental-stability and state-stability conditions coincide and depend on the sign of

$$ \rho^{(\text{loc})}(t) = \rho_{+}^{(\text{loc})}(t) = H_{11}^{(\text{loc})} = w_1(k-1) \frac{(k-2\tau^*)}{\tau(k-\tau^*)}, \quad (31) $$

which, being $k > 1$ and $\tau^* < \tau_f = k$, is positive for $\tau^* < \tau_s = k/2$ and negative for $\tau^* > k/2$. Hence, bifurcated branches originating at $\tau^* \in (1, k/2)$ are stable, whilst branches obtained for $\tau^* \in (k/2, k)$ are unstable, as marked by solid and dashed lines in Figure 2, respectively.

## 3.4 Discussion

Figure 2 displays three subtle phenomena typical of the evolution of systems with damage: (i) the presence of a continuous family of bifurcations along the homogeneous solution for $\tau \in (\tau_c, \tau_s)$; (ii) the presence of bifurcation points along the fundamental path notwithstanding its state-stability; (iii) distinct incremental-stability and state-stability thresholds. Similar behaviours are largely documented in the literature (Bigoni, 2012; Nguyen, 1994; Petryk, 1993) for other kind of inelastic systems. The paradigmatic example is the plastic
buckling of the Shanley's inelastic column (Shanley, 1947), where instabilities arise because of the interaction between the irreversible nature of perfect plasticity and geometrical nonlinearities. In the case of damage models presented here, geometrical nonlinearities are absent and instabilities are entirely due to the softening character of its constitutive law and the irreversibility constraint Benallal and Marigo, 2007.

A fundamental problem in this complex setting is to select among the infinite possible evolution paths those observable in physical experiments. The homogeneous solution for $\tau \in (\tau_c, \tau_s)$ in Figure 2 is an example of a state-stable solution of the time-continuous evolution problem, which is not a solution of the time-discrete evolution problem, being incrementally unstable. This solution will be hardly observable in an experiment because of the presence of bifurcated branches with lower energy. This is intuitive when looking at the values of the energy of each path reported in Figure 2-right. At the bifurcation point at the end of the elastic phase, there are two possible solution paths, the homogeneous and the localised solutions. Both paths are stable in the sense of the state-stability criterion (16). However, the localised solution is the path with the steepest energy descent, and thus preferable from the energetic point of view. On the basis of a similar reasoning, Petryk, 1993 introduced the concept of *path stability*, providing an in-depth analysis of the relation between the energetic stability concepts, the stability of evolution paths, and the dynamic definition of stability for inelastic systems, see also Fedelich and Ehrlacher, 1997. The definition of path stability is related to the property of the path to achieve, at each point, the steepest energy descent. The homogeneous path will be unstable according to this definition, whilst the localised solution will be stable. The incremental-stability condition (23), appears thus intimately related to path-stability because it selects among all admissible increments at each time, the one with minimal energy. Although state-stability is the appropriate definition for selecting *observable states*, independently of the evolution, it fails to be helpful to select the minimal energy direction at a bifurcation point and is unable to guarantee that a state-stable point can be reached as the result of an evolution through a minimal energy path.

In the following, we present and test numerical tools to test the *incremental stability* condition (23). This will allow us to detect bifurcation points and to select optimal evolution paths according to a minimal energy criterion on systems with many degrees of freedom encountered in the applications.

# 4 Numerical methods for stability and bifurcation analysis

Computational solvers for phase-field fracture models of the type (2) seek quasi-static evolutions as solutions of the time-discrete evolution problem (19). However, most of the solvers only impose first-order necessary conditions (22) without explicitly checking second-order local minimality conditions. The solution at the time step $t_i$ of the system of nonlinear variational inequalities (22) is typically performed iteratively by means of an energy-decreasing scheme, in its simplest incarnation, an alternate minimisation exploiting separate directional convexity of $u \mapsto E(u, \cdot)$ and $\alpha \mapsto E(\cdot, \alpha)$, as summarised in Algorithm 1 in the Appendix.

Here, we devise a numerical tool that enriches the standard alternate minimisation algorithm by performing the following tasks:

1. Given a solution $y_i$ of the first-order stability conditions (22), check the second-order time-discrete local minimality condition (23) to evaluate if the solution is incrementally-stable or unstable;
2. Construct stable evolution paths by introducing a continuation algorithm that bifurcates from unstable states, in the direction of minimal energy curvature, as soon as the evolution path ceases to be unique.

The numerical implementation of the stability assessment is summarised in Algorithm 2, whilst the continuation procedure is outlined in Algorithm 3. Both algorithms are reported in the Appendix to do not distract the reader. The basic ingredient of both algorithms is the solution of the minimisation problem appearing in the incremental stability condition (23). We compute stability in (23) by solving the eigenvalue problem (18) for the reduced Hessian $\tilde{H}_i$, the incremental stability indicator $\rho_i$ being computed as the smallest eigenvalue $\lambda_i^{(1)}$ of $\tilde{H}_i$. The stability-check Algorithm 2 takes as input a solution of the first-order optimality conditions and the solution at the previous time-step defining the irreversibility constraint. It returns the incremental-stability flag and the first eigenpair corresponding to the smallest eigenvalue of the reduced Hessian matrix. The continuation Algorithm 3 uses this information to provide a second-order stable solution $\mathbf{y}_i$ at the $i$-th time-step (given $\mathbf{y}_{i-1}$). First, it computes a tentative solution $\mathbf{y}_*$ by alternate minimisation (Algorithm 1) with $\mathbf{y}_{i-1}$ as initial guess. Hence, it checks its incremental-stability with Algorithm 2. If $\mathbf{y}_*$ is stable, it is accepted as solution at the $i$-th time-step. If not, the eigenvector $\mathbf{z}_*^{(1)}$ corresponding to the smallest negative eigenvalue $\lambda^{(1)}$ is used to generate a perturbed initial guess $\mathbf{y}_* + h \mathbf{z}_*^{(1)}$ for the alternate minimisation algorithm. The optimal amplitude $h_{\text{opt}}$ of the perturbation is selected through a bound constrained line-search procedure minimising the one-dimensional energy $E_i^{(1)}(h) := E(\mathbf{y}_* + h \mathbf{z}_*^{(1)})$ in the interval $h \in [h_{\text{min}}, h_{\text{max}}]$, where the bounds $h_{\text{min}}$ and $h_{\text{max}}$ are preliminarily computed to assure that $\alpha_{i-1} \le \alpha_* + h \beta_*^{(1)} \le 1$ for all $h \in [h_{\text{min}}, h_{\text{max}}]$. After checking its incremental-stability, the new solution of the alternate minimisation algorithm, if stable, is accepted as the solution at $i$-th time step. Possible failures to provide a stable solution happen when $h_{\text{min}} = h_{\text{max}}$ or when the newly calculated solution is unstable too. In those cases, several strategies could be conceived. For example, one can accept the unstable solution at the current time step, marking it as unstable, and proceed with the loading; if a stable solution is found at one of the next time-steps, this stable solution can be continued backward in time with a back-tracking strategy to provide a stable path at each loading step. As an alternative to a simple reiteration of the procedure, one can repeat Algorithm 3 using a different direction of negative curvature, as initial guess for the alternate minimisation, e.g. $\mathbf{z}_*^{(2)}$ or $\mathbf{z}_*^{(k,1)}$. We focus here on the presentation on the basic algorithm, leaving similar additional developments for future works.

In practice, the finite element discretisation of the damage energy functional (2) is performed using the finite-element framework FEniCS (Alnaes et al., 2015), which includes symbolic differentiation tools to calculate the second derivative of the energy and to automatically assemble finite element matrix and vector using PETSc (Balay et al., 2020) as the linear-algebra back end. PETSc provides also efficient large-scale solvers (SNES, TAO) for the solution of the elastic and the damage subproblems of the alternate minimisation Algorithm 1. The reduced Hessian matrix $\tilde{\mathbf{H}}_i$ is obtained by extracting from the full Hessian $\mathbf{H}_i$ the block matrix corresponding to the index set $N_i(\mathbf{y}_i)$, using PETSc's index-set (IS) objects and the related operators. The smallest eigenvalues and corresponding eigenvectors are computed using the large-scale eigenvalue solvers provided by SLEPc (Hernandez et al., 2005). Namely, we use the Krylov-Schur iterative algorithm with the shift-invert option and a small negative spectral shift (Roman et al., 2020). We also implement, as an alternative method to check the time-discrete stability condition, a direct evaluation of the inertia of the matrix $\tilde{\mathbf{H}}_i$ through a Cholevsky factorisation, which yields the number of negative, vanishing, and
<figure><img src="image_5.png" /><figcaption>Figure 4: Geometry and loading for the uniaxial traction test. We use unstructured simplicial mesh of uniform mesh size $h$. For a give internal length $l$, we typically set $h = l/5$.</figcaption></figure>

positive eigenvalues of the matrix. This information is returned directly by the parallel sparse direct solver MUMPS (Amestoy et al., 2019) when activating the specific options INFOG(12) and ICNTL(13)=1. If there is at least one negative eigenvalue, the state is marked as unstable. This achieves the task without explicitly computing the solution to the large and sparse second-order eigenvalue problem.

# 5 Numerical experiments

We illustrate the concepts and numerical tools introduced in this paper through the solution of two basic examples, focussing of the interplay of homogeneous and localised solutions and their stability properties. While we found that in most of the situations the alternate minimisation algorithm provides incrementally stable evolutions, we report in the following specific counterexamples showing that this is not true in general and that the proposed stability-check and continuation algorithms provide an useful addition to the existing numerical techniques. To this end, we present (i) the traction test of a slender bar and (ii) the multifissuration of a pre-stressed thin film. In these examples, we do not impose any Dirichlet boundary condition on the damage field, in order to allow for non-vanishing homogeneous solutions of the damage field.

For each of the following illustrative examples, at each time step, we determine the solution of first order necessary conditions (22) with the alternate minimisation Algorithm 1 and we perform the stability check with Algorithm 2 to verify the incremental stability of the computed critical states. This is done by checking the second order local minimality condition (23). Hence, we compare evolution obtained with or without activating the continuation Algorithm 3.

## 5.1 Uniaxial traction test

Let $\Omega = (0, L) \times (-H/2, H/2)$ be the reference configuration of a brittle bar of length $L$ and width $H$. The displacement in the axial direction is blocked on the left end and imposed to $\bar{u}_t = tL$ at the right end, with a free sliding condition, see Figure 4. To eliminate rigid body modes, we block both components of displacement at the lower-left corner. Conversely, we leave damage free on the entire boundary to allow both homogeneous states of damage as well as boundary cracks.

We consider the damage model described by the energy functional defined in (1)-(2) with the constitutive choice $a(\alpha) = (1 - \alpha)^2$ and $w(\alpha) = \alpha$, corresponding to a model which is frequently adopted in phase-field fracture, referred to as $\mathbf{AT}_1$ in (Tanné et al., 2018). Setting the elastic Lamé coefficients to reproduce plane-stress conditions, the energy functional (2)
reads as

$$ \mathcal{E}(u, \alpha) = E \int_{\Omega} (1 - \alpha)^2 \left( \frac{\nu \operatorname{tr}(e(u))^2}{2(1 - \nu^2)} + \frac{\|e(u)\|^2}{2(1 + \nu)} \right) dx + w_1 \int_{\Omega} (\alpha + \ell^2 |\nabla \alpha|^2) dx. \quad (32) $$

The uniaxial traction test has been largely discussed in previous analytical (Pham et al., 2011a) and numerical (Pham et al., 2011b) works, highlighting the competition between a solution with homogeneous damage and a solution where damage localises in a band transverse to the long axis. The width of the localisation band is proportional to the internal length $\ell$. The localised solution can be regarded as the phase-field approximation of a crack, for which the elastic energy and the stress vanish whilst the dissipated energy per unit length is $G_c = \frac{8}{3}w_1\ell$, see (Pham et al., 2011b; Tanné et al., 2018). A one-dimensional bar model allows for an analytical stability and bifurcation analysis of the homogeneous solution. As reported in (Pham et al., 2011a), the homogeneous solution is purely elastic for an imposed displacement $t < t_c$ and a stress $\sigma < \sigma_c = \sqrt{w_1E}$. After this first critical threshold, the damage level of the homogeneous solution monotonically grows with the loading $t$. Testing the infinite dimensional version of the non-bifurcation and the state-stability criterion (17)-(20), one can determine the loading intervals $t \in (0, t_b)$ and $t \in (0, t_s)$ with $t_s \ge t_b$ for which the homogeneous solution is locally unique (i.e. incrementally-stable) and state-stable, respectively. The critical loads are given by (see Pham et al., 2011a, Example 1 with $p = 0$, $q = 2$):

$$ t_c = \sqrt{\frac{w_1}{E}}, \quad \sigma_c = \sqrt{w_1 E}, \quad \frac{t_b}{t_c} = \min \left\{ 1, \pi \sqrt{\frac{2}{3}} \frac{\ell}{L} \right\}, \quad \frac{t_s}{t_c} = \min \left\{ 1, \frac{4\pi}{3} \sqrt{\frac{2}{3}} \frac{\ell}{L} \right\}. \quad (33) $$

With this parametrisation using the specific fracture energy $w_1$ instead of the fracture toughness $G_c$ to represent the dissipative properties of the material, the critical load for the onset of damage is independent of the internal length $\ell$. However, the stability threshold critically depends upon the ratio between the internal length and the size of the specimen. The response is characterized by a size effect distinguishing two different regimes, that of long and of short bars, in relative terms with regards to the material internal length. We show that these two regimes feature different qualitative properties in terms of stability of the homogeneous response.

We perform a series of simulations varying the internal length $\ell$, with $\ell \in [L/5, L]$, setting the aspect ratio $H/L = 1/10$ and the Poisson ratio $\nu = 0.3$. The dependence upon other material parameters is eliminated by rescaling the displacement by $u_0 = L\sqrt{w_1/E}$. For each simulation, we monotonically increase the loading from an undamaged state at $t = 0$. Figure 5 shows the graph of the parametric stability diagram as a function of the relative length $L/\ell$ and the loading $t$. This diagram serves as a verification experiment for the numerical approach against the analytic solution of (Pham et al., 2011a). Black curves represent the analytical state-stability and non-bifurcation thresholds in Equation (33). These curves partition the state space in regions where, according to the analytical prediction, the state is: (i) purely elastic, in blue; (ii) with non vanishing damage and both incrementally- and stable-stable, in orange; (iii) incrementally-unstable and state-stable, in white; (iv) both incrementally- and state-unstable, in green. Each horizontal set of dots corresponds to an irreversible evolution computed for a given value of $L/\ell$ and each time step is represented with a dot: blue dots for purely elastic states and orange dots for incrementally-stable states, i.e. such that $\rho_i > 0$
<figure><img src="image_6.png" /><figcaption>Figure 5: Stability diagram for the uniaxial traction test of Figure 4 when increasing the loading $t$ as a function of the relative length of the bar $L/\ell$. Each horizontal array of symbols represents the result of a numerical experiment at the different loading steps for a given value of $L/\ell$. Dots indicate stable states and crosses indicate incrementally unstable states ($\rho_i < 0$). The vertical line indicates the elastic limit load $t_c$, the thick curve is the analytic limit of stability $t_s$, and the dash-dotted curve represents the incremental-stability threshold $t_b$ given in Equation (33). These lines partition the diagram in four regions, whose colours indicate the stability properties of the homogeneous solution according to the analytical predictions: stable states in blue, both incrementally- and state-stable in orange, both incrementally- and state-unstable states in green, incrementally-unstable but state-stable states in white.</figcaption></figure>

within the numerical tolerance<fnref n="2" />. The simulation is halted upon encountering an incrementally unstable state with $\rho_i < 0$, represented with a cross. Overall, the numerical results are in very good agreement with the theoretical predictions.

Figure 6 reports the evolution of the stress, the energy, and the smallest eigenvalue of the Hessian matrix for the case of a relatively short bar ($L/\ell = 1.7$). The left column shows the solution obtained with by running the plain alternate minimisation algorithm at each time-step. The presence of negative eigenvalues (in red in the bottom diagram) indicates that the alternate minimisation returns an incrementally unstable solution for $t > t_b$. The corresponding solution is a stationary point of the energy which verifies first-order optimality conditions, but fails to be a local energy minimizer in the sense of (19). This is a clear issue of the plain alternate minimisation algorithm, commented also in (Bourdin, 2007), because unstable solutions are physically unobservable. The diagrams in the right column report the results of the same simulation obtained when activating the continuation Algorithm 3. In this case, when the stability check detects an unstable state, we switch from the current evolution path towards a stable state at the current time step. A new alternate minimisation solver is launched, seeded by a perturbation obtained from the eigenmodes as detailed in Algorithm 3. The proposed continuation algorithm improves the results of the plain alternate minimisation, being able to retrieve an evolution made of incrementally-stable states.

<footnote n="2">²in practice, we get the smallest eigenvalue $\rho_i$ from the solver with a numerical relative tolerance $rtol = 10^{-8}$.</footnote>
Figure 7 shows the analogous results for the case of a long bar, where the theoretical bifurcation threshold coincides with both the elastic limit and the state-stability threshold ($t_s = t_b = t_c$). This case is more interesting for phase-field fracture applications, where typically $L/\ell \gg 1$. We report here the case $L/\ell = 5$, the qualitative behaviour being similar for larger values of $L/\ell$. As for short bars, the plain alternate minimisation returns unstable states, with a jump in the response, corresponding to a localisation of damage, only for a loading significantly larger than the theoretical stability threshold. Vice-versa, activating the continuation algorithm, a bifurcation into the direction of minimal energy curvature allows the algorithm to follow a stable evolution, as soon as the damage criterion is met.

Figure 8 summarises the key results, showing the evolution of the damage field in time and space along the central axis for short and long bars, with and without the continuation Algorithm. For the given boundary conditions (free damage at the ends), localisations and instabilities are characterised by the appearance of one single boundary crack which can indifferently be at the left or right end.

The mechanism of loss of stability is different in the case of short and long bars. For the former (Figure 6), the instability happens at a loading strictly larger than the elastic limit, when the smallest eigenvalue smoothly crosses zero in correspondence to a non-vanishing homogeneous damage level. For long bars (Figure 7), the instability happens immediately after the elastic limit, when one or more negative eigenvalues suddenly appear as the irreversibility constraints become inactive and the corresponding space $\mathcal{N}(\mathbf{y}_i)$ in (23) suddenly changes.

Figure 9 illustrates some details of the stability-check and continuation algorithms for the case of the long bar, for which several negative eigenvalues are detected immediately after $t = t_c$. Figure 9-left shows the damage distribution along the bar longitudinal axis for the two eigenmodes $\beta_*^{(1)}$, $\beta_*^{(2)}$ associated with the two smallest negative eigenvalues. They correspond to damage localisations at a boundary (in grey), at the centre of the bar (in black). The boundary localisation $\beta_*^{(1)}$ is associated with the smallest eigenvalue and is energetically preferable, when admitted by boundary conditions. In this specific case, its corresponding eigenspace is of multiplicity two, because of the possibility of localisations at both ends, due to the symmetry of the system. Figure 9-right shows the optimal amplitude of the perturbation for the damage along the minimal curvature direction $\beta_*^{(1)}$. This is determined through the constrained line-search formulated in Algorithm 3. The inset shows the variation of energy $E_i^{(1)}(h) := E_i(\mathbf{y}_* + h\mathbf{z}_*^{(1)})$ along the direction pointed by the first eigenmode, within the corresponding admissible lower and upper bounds obtained from the constrained line-search procedure. This perturbation technique allows us to construct evolutions composed by states which are incrementally stable, in a way that is more pertinent from the physical point of view. Still, the solution to the evolution problem is not necessarily unique. For example, we are able to reach a different local energy minima perturbing the unstable homogeneous state with the second eigenmode, leading to a fully developed crack in the centre of the specimen (not reported here). As it can be checked analytically, the (optimal) profile of the boundary crack is half the symmetric profile associated to a macroscopic crack in the bulk. Consequently, the total energy of the state with one bulk crack in the bulk is twice that of the bar with a boundary crack.
<figure><img src="image_7.png" /><figcaption>Figure 6: Traction experiment for a relatively short bar ($L/ℓ = 1.7$). The plots report the numerical evolution of the average axial stress $\sigma$ across a cross-section (top), of the energy (middle), and of the smallest eigenvalues (bottom) as a function of the loading, for $t \in (0, 3)$. The left column is obtained exploiting the plain alternate minimisation Algorithm 1, whilst the right panel shows the solution provided by the continuation Algorithm 3. Eigenvalues and stability of states are computed with Algorithm 2. The regions shaded in light blue, orange, and green, indicate respectively purely elastic, homogeneously damaging and bifurcated (localised) regions, according to the analytical stability analysis; the white region covers homogeneous states which are state-stable yet incrementally-unstable. In the bottom row we mark in light green positive eigenvalues and in red negative ones, reporting also the maximum value of damage (blue curve, right axis).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: Traction experiment for a relatively long bar ($L/\ell = 5.0$). Legend as in Figure 6. For a long bar, the state of homogeneous damage is unstable (state-wise and incrementally) as soon as the damage criterion is met. Multiple negative eigenvalues (cf. left panel, bottom) correspond to potential directions of energy decrease. In the right panel, the continuation Algorithm 3 delivers a second-order stable evolution featuring the emergence of one sudden crack at $t = t_c = t_b = t_s$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>(a) Short bar - plain alternate minimisation</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Short bar - Continuation algorithm.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Long bar - plain alternate minimisation</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(d) Long bar - Continuation algorithm.</figcaption></figure>

Figure 8: Bar in uniaxial traction: space-time evolution of damage along the longitudinal axis of the bar for the short and long cases, both with and without continuation Algorithm 3. Only the continuation algorithm is able to retrieve a stable evolution (see Figures 6-7). Colours indicate the damage level from 0 (sound material, violet) to 1 (fully cracked, red).
<figure><img src="image_13.png" /><figcaption>Figure 9: Perturbations and energy landscape for the traction problem of a long bar of Figure 8d immediately after the elastic limit $t_c$. Left: first two normalised eigenmodes associated to negative eigenvalues. Right: optimal perturbation along the first eigenmode used for seeding the alternate minimisation Algorithm 1 at the occurrence of an unstable solution, according to Algorithm 3. The inset shows the energy $E_i^{(1)}(h) := E_i(y_* + hz_*^{(1)})$ along the direction of the first mode, used to determine the optimal amplitude subject to the constraint on $h: \alpha_{i-1} \le \alpha_* + h\beta_*^{(1)} \le 1$. We represent damage along the beam axis ($y=0$), the damage field being invariant in the transversal $y$-direction.</figcaption></figure>

## 5.2 Multifissuration in thin films

We now consider a two-dimensional model of a thin film bonded to a substrate by the means of a compliant layer, modelling this system as an effective two-dimensional membrane resting onto a linearly elastic foundation. The film is brittle in the sense that it undergoes an irreversible damaging process under tensile stress, eventually leading to fracture. In this case, the interplay between tensile membrane deformations, possibly released by transverse cracking, and the energy stored in the foundation, produces a rich phenomenology leading to the emergence of robust and variegated fracture patterns (Xia and Hutchinson, 2000). Indeed, differently from the one-dimensional bar in traction, the system at hand features an additional characteristic length scale of purely elastic origin which accounts for the interaction between the membrane and the underlying substrate. We focus on the material regime in which the elastic length scale is small compared to the size of the domain and large with respect to the damage length. In this regime, the loss of stability of the sound state results in the emergence of fracture patterns characterised by robust morphological properties (Leon Baldelli et al., 2014).

Let $\Omega_{3d} = \Omega \times [-T/2, T/2] \subset \mathbb{R}^3$ be the reference configuration of a thin film of diameter $L := |\Omega|^{1/2}$, mid-plane $\omega \subset \mathbb{R}^2$, and thickness $T$. The film is in plane-stress conditions, as indeed can be inferred via an asymptotic dimension-reduction procedure, and undergoes in-plane displacements $u: \omega \mapsto \mathbb{R}^2$ which are elastically constrained on a rigid substrate. Following (Leon Baldelli et al., 2014; Xia and Hutchinson, 2000), the elastic interaction between the thin film and the substrate is modelled as a linear elastic foundation of stiffness $K > 0$, which introduces an additional elastic energy density $\frac{K}{2}|u|^2$. As common in thin film applications, we suppose that the membrane is loaded by a given inelastic deformation tensor $\bar{e}_t$, that can model, for example, the differential thermal expansion between the thin-film and the substrate. Adding to the plane-stress total energy functional (32) the contribution due to the presence of the substrate, we associate to $(u, \alpha) \in H^1(\Omega, \mathbb{R}^2) \times H^1(\Omega)$ the following
non-dimensional energy per unit of thickness

$$ \mathcal{F}(u, \alpha) = E \int_{\Omega} \left( (1-\alpha)^2 \left( \frac{\nu \operatorname{tr}(e(u) - \bar{e}_t)^2}{2(1-\nu^2)} + \frac{|e(u) - \bar{e}_t|^2}{2(1+\nu)} \right) + \frac{1}{\ell_e^2} \frac{|u|^2}{2} \right) dx \\ + w_1 \int_{\Omega} (\alpha + \ell^2 |\nabla \alpha|^2) dx, \quad (34) $$

where $\ell_e := \sqrt{ET/K}$ is the relative weight of the elastic energy stored in the membrane relative to the elastic foundation. We focus on the stability properties of the homogeneous solution when the elastic internal length is small compared to the size of the structure, namely $\ell_e \ll L$, corresponding to a stiff interaction between the film and the substrate. Heuristically, in this regime the displacement field $u$ is close to zero in the entire domain (possibly up to a boundary layer of size $\ell_e$) because the foundation term in the energy has a large multiplicative coefficient. As a consequence, the stress is approximately constant and the damage criterion is met almost uniformly (up to a boundary layer). This, coupled with the existence of the additional length scale $\ell$, is responsible for the appearance of structured crack patterns beyond the elastic limit. As for (32), without loss of generality, one can set $E = 1$, $w_1 = 1$. This is equivalent to rescaling the displacement field by the scaling factor $u_0 = L\sqrt{w_1/E}$. Hence, the relevant dimensionless material parameters appearing in (34) are the relative elastic length scale $\ell_e/L$, the relative internal length for the damage field $\ell/L$, and the Poisson ratio $\nu$.

We report a first numerical experiment for a film in the form a thin strip $\Omega \equiv (-L/2, L/2) \times (-H/2, H/2)$, as in Figure 4, with $L = 6$ and $H = 0.1$, by setting material parameters to $\ell/L = 0.1$, $\ell_e/L = 0.3$, $\nu = 0$. The thin-film strip has fixed, homogeneous, horizontal displacements at the left and right ends $x = \pm L/2$, is free on the top and bottom boundaries, $y = \pm H/2$, whilst no Dirichlet boundary conditions are imposed on the damage field, which consequently satisfies natural homogeneous Neumann conditions. To focus on the almost one-dimensional regime and simplify the interpretation of the results, we consider the case of an uniaxial loading $\bar{e}_t = t\hat{e}_x \otimes \hat{e}_x$, where $\hat{e}_x$ stands for the unit vector in the axial direction of the bar. With this almost one-dimensional setting, the solution for the damage field is invariant with respect to the transverse direction and depends only on the axial variable $x \in (-L/2, L/2)$. Figure 10 compares the results of the plain alternate minimization Algorithm 1 (left) with those of the proposed continuation Algorithm 3 (right) for the evolution problem obtained when monotonically increasing the loading $t$. In each column, from top to bottom, we report the evolution in time of the space distribution of the damage field along the axis of the strip, the energies, and the minimal eigenvalue of the reduced Hessian matrix used as a stability marker. The plain alternate minimisation Algorithm 1 (left) delivers a uniform state of damage beyond the critical load $t_c$, which is unstable for $t > t_c = \sqrt{w_1/E}$, from which multiple clustered negative eigenvalues appears. The associated eigenvectors are functions oscillating along the axis of the thin-film strip with a specific wavelength, of the order of $\ell_e$. Figure 11-left shows the first of such modes provided by the stability check algorithm. This mode is used by the continuation Algorithm to perturb the homogeneous solution in the direction of energy descent. Following this path, it retrieves a stable solution which shows equally spaced damage localisations, as shown in Figure 10-right. This localisation pattern can be regarded as the phase-field approximation of a periodic crack array, which is typical of thin-film systems (Leon Baldelli et al., 2013; Xia and Hutchinson, 2000). Interestingly the spacing of this crack pattern is related to the wave-length of the unstable eigenmodes.
<figure><img src="image_14.png" /><figcaption>Figure 10: Thin film strip with $L = 6$, $H = 0.1$, $\ell = 0.1$, $\ell_e = 0.3$. Solution of the evolution problem when monotonically increasing the pre-stress $t$, obtained with plain alternate minimization (left column) and with the proposed path-selection Algorithm (right column). Top: Damage field $\alpha$ along the central axis as a function of the space variable $x$ and the time $t$. Centre: energies. Bottom: minimal eigenvalues resulting from the stability analysis. The critical load in the figures corresponds to the elastic limits $t_c = \sqrt{w_1/E}$. The plain alternate minimisation retrieves an unstable homogeneously damaged solution for $t > t_c$, whilst the continuation algorithms triggers a sudden jump of the solution toward the stable branch at $t = t_c$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 11: Left: Detail of the bifurcation at the onset of damage, showing the unstable homogeneous state $\alpha_*$, the irreversibility lower bound $\alpha_{i-1}$ and the optimally perturbed state $\alpha + h_{opt}\beta^{(1)}$, see Algorithm 3. Right: Fully localised solution for the damage field $\alpha$ obtained at the end of the simulation.</figcaption></figure>

The two-dimensional counterpart of the thin strip shows the emergence of complex two-dimensional quasi-periodic crack patterns, stemming from a sound and homogeneous state of damage. We consider a thin film occupying a disk of radius $R$ loaded by a uniform equiaxial $\bar{e}_t = tI$. The displacement and the damage are left free on the boundary. The brittle system is in its undamaged state for $t < t_c$ where, here, $t_c := \frac{(1-\nu)}{2}\sqrt{w_1/E}$. For the given boundary conditions, a solution with homogeneous non-vanishing damage respects the first-order optimality conditions for $t > t_c$. However, the stability of such state is conditional upon the values of material parameters. In particular, the homogeneous solution is unstable in the regime $\ell_e \ll R$ for $t \ge t_c$, in the sense that there exists a family of negative eigenvalues solving (23) as soon as the damage criterion is attained. This is a challenging scenario from the numerical standpoint because of the homogeneity of the solution, the existence of a large number of local minima whose features span several length scales, and the large dimension of the space generated by eigenvectors associated to negative eigenvalues. We report in Figure 12 the result of a simulation with $R = 1$, $\ell_e = 0.15$, $\ell = 0.03$, $\nu = 0.3$, and a uniform mesh size $h = \ell/3$. For $t = t_c$, the stability-check Algorithm detects 189 negative eigenvalues, computing the inertia of the reduced Hessian. The solution of the full second order eigenvalue problem returns the smallest eigenvalues and the corresponding eigenvectors, see the spectrum in Figure 12-left. The continuation algorithm 3 thus seeks a stable state departing, for $t = t_c$, from the homogeneous critical point obtained imposing first order optimality conditions. We represent in Figure 12 (top-right) an elevated plot of the damage component of the first mode at $t \simeq t_c$, which is associated to the smallest negative eigenvalue $-6.7 \times 10^{-5}$. Seeding a new first order alternate minimisation with an optimal perturbation of the homogeneous state in the direction of the first eigenmode leads to an energy decreasing path and a new, stable, damaged state. The damage field at the end of the simulation ($t \simeq 1.7t_c$) is shown in Figure 12 (bottom-right). Damage, structured in bands (of the size of $\ell$), forms a quasi-periodic network of intersecting fractures whose features are reminiscent of the eigenmodes at the bifurcation point. We leave the quantitative analysis of fracture patterns in their relation to effective unstable modes to a future work.
<figure><img src="image_16.png" /><figcaption>Figure 12: Left. Energy diagram (top) and spectrum of the second derivative (bottom), showing a stable irreversible evolution departing from the undamaged state. The small disks on the top show few snapshots of the damage field during the evolution ($\alpha = 1$ in black and $\alpha = 0$ in light grey). At the critical bifurcation point $t_c$, the lowest negative eigenmode (top-right) induces damage localisation on a scale correlated to the internal elastic length $l_e$. The continuation Algorithm 3 returns a stable evolution with a spectrum bounded away from zero, along the imposed increasing load, showing further fracture nucleation and branching. At the end of the loading program, the damage field shows a two-dimensional fracture pattern (bottom-right). The unstructured computational mesh, with approximately 100K degrees of freedom, is indicated in transparency. Parameters for this computation are $R = 1$, $l_e = 0.15$, $\ell = 0.03$, $E = 1$, $\nu = 0.3$.</figcaption></figure>

# 6 Conclusions

We have discussed the numerical detection of instabilities and possible bifurcations in quasi-static evolutions of a class of softening damage models that is used in phase-field approach to fracture. Our approach is based on the variational formulation of the evolution problem as a rate-independent process. The stability of the solution is defined as an energy minimality requirement of the damage energy functional under the irreversibility constraint on the damage variable. We distinguished between the notion of *state-stability* and *incremental-stability*, associated to the bound-constrained energy minimality condition in the time-continuous and the time-discrete formulations of the evolution problem, respectively. We have shown that the incremental-stability is a stronger requirement than the state-stability condition. For smooth evolutions, the incremental-stability condition coincides with the non-bifurcation criterion of the time-continuous rate-problem. The different concepts were illustrated through the analytical solution of a two degrees-of-freedom model problem. This example explicitly highlights the possible presence of admissible energy decreasing evolution paths departing from state-stable
solutions. This led us to suggest the incremental-stability condition as a more appropriate
criterion to select the evolution paths that are actually observable in physical experiments,
echoing previous works of (Petryk, 1993) for plasticity.

Hence, we discussed the numerical techniques to practically test the increment-stability by
checking the second-order optimality condition of a bound-constrained minimisation problem
for the total energy functional. This requires the numerical solution of a coupled eigenvalue
problem in the displacement and damage variables for a reduced Hessian matrix, which is
obtained by eliminating the degrees-of-freedom associated to active irreversibility constraints
on damage. The sign of the smallest eigenvalue determines the incremental stability condition.
Eigenvectors associated to negative eigenvalues are directions with negative curvature of the
energy and represent possible modes of instability. We have shown by explicit examples
that the standard alternate minimisation algorithm, testing only the first-order minimality
condition, can converge to incrementally unstable solutions. Hence, we proposed a basic path
selection algorithm that filters out unstable paths. When detecting unstable states, it re-
initialises the alternate minimisation solver with a perturbed seed, obtained by a line-search
procedure along the direction of minimal energy curvature.

The numerical framework has been tested on a classical traction bar problem and on a
more complex thin film multifissuration problem. The traction bar problem, for which explicit
analytical solutions are available, serves as a verification test-case. We have shown that the
presented algorithm allows the detection of unstable evolution paths retrieved with standard
first-order algorithms. Finally, the thin film problem illustrates the interest of the proposed
stability and bifurcation analysis on more complex problems, where the unstable eigenmodes
also provide useful information on the emerging crack pattern and on its mechanism of for-
mation.

We believe that the proposed techniques and concepts are useful for the solution of prac-
tical crack nucleation problems in a computational setting. Further developments should
include the improvement of the numerical stability check and path selection algorithms, cou-
pling the present algorithm to more advanced path following techniques (see e.g. Cochelin
et al., 2007; Haslinger et al., 2018; Lorentz and Godard, 2011).
# 7 Appendix

We list here the three algorithms introduced in the article for the solution of the first-order optimality condition (Algorithm 1), the numerical stability and bifurcation analysis (Algorithm 2), and for the continuation procedure (Algorithm 3).

**Algorithm 1:** Alternate minimisation

**Data:** The irreversibility constraint $\alpha_{i-1}$ and the initial guess $\alpha^0$

**Result:** Solution $\mathbf{y}_i = (\mathbf{u}_i, \alpha_i)$ of the first-order optimality conditions (22) at the time step $t_i$

$k \leftarrow 1$

while not converged do

$\mathbf{u}^k \leftarrow \arg \min_{\mathbf{u}} E_i(\mathbf{u}, \alpha^{k-1})$

$\alpha^k \leftarrow \arg \min_{\alpha \ge \alpha_{i-1}} E_i(\mathbf{u}^k, \alpha)$

$k \leftarrow k + 1$

end

$\mathbf{y}_i = (\mathbf{u}_i, \alpha_i) \leftarrow (\mathbf{u}^k, \alpha^k)$

**Algorithm 2:** Incremental stability check (I-StabilityCheck)

**Data:** The critical state $\mathbf{y}_i$ satisfying the first-order optimality conditions (22) and the irreversibility constraint $\alpha_{i-1}$

**Result:** Time discrete stability of the state $\mathbf{y}_i$ (I-Stable) and, if unstable, the first $n_{\text{eig}}$ smallest eigenpairs $\{\lambda_i^{(k)}, \mathbf{z}_i^{(k)}\}_{k=1}^{n_{\text{eig}}}$.

Compute the inactive index set $\mathcal{N}_i(\mathbf{y}_i)$

if $\mathcal{N}_i(\mathbf{y}_i) = \emptyset$ then

$\mathbf{y}_i$ is purely elastic I-Stable $\leftarrow$ True

else

/* assemble the Hessian $H_i$ and extract the reduced Hessian $\tilde{H}_i$ */

/* calculate the number of negative ($n_{\text{neg}}$) and vanishing ($n_0$) eigenvalues of $\tilde{H}_i$ via MUMPS' Cholesky factorisation and getInertia operator in the PETSc interface */

$n_{\text{neg}}, n_0, n_{\text{pos}} \leftarrow \text{getInertia}(\tilde{H}_i(\mathbf{y}_i))$;

if $n_{\text{neg}} = n_0 = 0$ then

$\text{I-Stable} \leftarrow$ True

else if $n_{\text{neg}} > 0$ or $n_0 > 0$ then

$\text{I-Stable} \leftarrow$ False;

/* solve the reduced eigenvalue problem (18) extending eigenvectors by zero to the full space $\mathbb{R}^n */

$\{\lambda_i^{(k)}, \mathbf{z}_i^{(k)}\}_{k=1}^{n_{\text{eig}}} \leftarrow \text{eigen}(\tilde{H}_i(\mathbf{y}_i))$

end
**Algorithm 3:** Time-discrete continuation algorithm seeking a second-order stable evolution.

**Result:** $\mathbf{y}_i = (\mathbf{u}_i, \boldsymbol{\alpha}_i)$: A solution respecting the second-order optimality condition

**Data:** $\mathbf{y}_{i-1} = (\mathbf{u}_{i-1}, \boldsymbol{\alpha}_{i-1})$: The solution at the previous time step

/* Run the alternate minimisation with $\boldsymbol{\alpha}_{i-1}$ as irreversibility constraint and initial guess */
$\mathbf{y}_* \leftarrow \text{AltMin}(\boldsymbol{\alpha}_{i-1}, \boldsymbol{\alpha}_{i-1})

/* Calculate the time-discrete stability of $\mathbf{y}_i$ with $\boldsymbol{\alpha}_{i-1}$ as irreversibility constraint */
I-Stable, $\{\lambda_*^{(k)}, \mathbf{z}_*^{(k)}\}_{k=1}^{n_{\text{eig}}} \leftarrow \text{I-StabilityCheck}(\mathbf{y}_*, \boldsymbol{\alpha}_{i-1})

if I-Stable = True then
$\mathbf{y}_i \leftarrow \mathbf{y}_* / *$ accept $\mathbf{y}_*$ as solution at the time step $i$ */

else

/* Find the optimal amplitude of the perturbation by solving the constrained line-search problem in the direction of minimal curvature $\mathbf{z}_*^{(1)} = (\mathbf{v}_*^{(1)}, \boldsymbol{\beta}_*^{(1)})$ */
repeat
  $h_{\text{opt}} \leftarrow \arg \min_h \{E_i^{(1)}(h) := E_i(\mathbf{y}_* + h \mathbf{z}_*^{(1)}), h: \boldsymbol{\alpha}_{i-1} \le \boldsymbol{\alpha}_* + h \boldsymbol{\beta}_*^{(1)} \le 1\}$

if $h_{\text{opt}} \ne 0$ then
  /* Run the alternate minimisation with $\boldsymbol{\alpha}_{i-1}$ as irreversibility constraint and $\boldsymbol{\alpha}_* + h_{\text{opt}} \boldsymbol{\beta}_*^{(1)}$ as initial guess */
  $\mathbf{y}_*^{(1)} \leftarrow \text{AltMin}(\boldsymbol{\alpha}_{i-1}, \boldsymbol{\alpha}_* + h_{\text{opt}} \boldsymbol{\beta}_*^{(1)})

I-Stable, $\{\lambda_*^{(1,k)}, \mathbf{z}_*^{(1,k)}\}_{k=1}^{n_{\text{eig}}} \leftarrow \text{I-StabilityCheck}(\mathbf{y}_*^{(1)}, \boldsymbol{\alpha}_{i-1})$

if I-Stable = True then
  $\mathbf{y}_i \leftarrow \mathbf{y}_*^{(1)} / *$ accept $\mathbf{y}_*^{(1)}$ as solution at the time step $i$ */

else

/* Unable to get a stable solution, issue a warning and return the minimal energy solution (alternatively, restart the procedure with $\mathbf{z}_*^{(1,k)}$) */
  $\mathbf{y}_i \leftarrow \arg \min(E_i(\mathbf{y}_*), E_i(\mathbf{y}_*^{(1)}))$

else

/* unable to get a stable solution, issue a warning and return
   $\mathbf{y}_*$ (alternatively, restart the same procedure with $\mathbf{y}_*^{(2)}$ or $\mathbf{z}_*^{(1,k)}$, or use an infeasible solution as initial guess) */
  $\mathbf{y}_i \leftarrow \mathbf{y}_*$

until I-Stable = True;

end

$i \leftarrow i + 1 / *$ update the time step */;
References

Alnaes, M. et al. (2015). "The FEniCS Project Version 1.5". en. In: *Archive of Numerical Software* 3.100.

Ambrosio, L. and V. M. Tortorelli (1992). "On the approximation of Free Discontinuity Problems". In: *Bollettino dell'Unione Matematica Italiana* 7.6-B, pp. 105-123.

Amestoy, P., A. Buttari, J.-Y. L'Excellent, and T. Mary (2019). "Performance and Scalability of the Block Low-Rank Multifrontal Factorization on Multicore Architectures". In: *ACM Transactions on Mathematical Software* 45 (1), 2:1-2:26.

Amor, H., J.-J. Marigo, and C. Maurini (2009). "Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments". In: *Journal of the Mechanics and Physics of Solids* 57.8, pp. 1209-1229.

Balay, S. et al. (2020). *PETSc Users Manual*. Tech. rep. ANL-95/11 - Revision 3.13. Argonne National Laboratory.

Bazant, Z. P. and M. R. Tabbara (1992). "Bifurcation and stability of structures with interacting propagating cracks". In: *International Journal of Fracture* 53.3, pp. 273-289.

Beaurain, J. (2011). "Recherche de solutions bifurquées et étude de leur stabilité dans les problèmes d'endommagement". PhD thesis. Université Pierre et Marie Curie-Paris VI.

Beaurain, J., K. Kazymyrenko, and J.-J. Marigo (2011). "Méthode numérique pour l'étude de la stabilité de modèles d'endommagement à gradient". In: *CSMA 2011*.

Benallal, A. and J.-J. Marigo (2007). "Bifurcation and stability issues in gradient theories with softening". In: *Modelling and Simulation in Materials Science and Engineering* 15, pp. 283-295.

Benallal, A., R. Billardon, and G. Geymonat (1993). "Bifurcation and Localization in Rate-Independent Materials. Some General Considerations". In: *Bifurcation and Stability of Dissipative Systems*. Ed. by Q. S. Nguyen. Vienna: Springer Vienna, pp. 1-44.

Bigoni, D. (2012). *Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability*. Cambridge University Press.

Bleyer, J. and R. Alessi (2018). "Phase-field modeling of anisotropic brittle fracture including several damage mechanisms". In: *Computer Methods in Applied Mechanics and Engineering* 336, pp. 213-236.

Bourdin, B. (2007). "Numerical implementation of the variational formulation for quasi-static brittle fracture". In: *Interfaces and Free Boundaries* 9, pp. 411-430.

Bourdin, B., G.-A. Francfort, and J.-J. Marigo (2000). "Numerical experiments in revisited brittle fracture". In: *Journal of the Mechanics and Physics of Solids* 48, pp. 787-826.

Bourdin, B., J.-J. Marigo, C. Maurini, and P. Sicsic (2014). "Morphogenesis and Propagation of Complex Cracks Induced by Thermal Shocks". In: *Physical Review Letters* 112, p. 014301.

Cochelin, B., N. Damil, and M. Potier-Ferry (2007). *Méthode asymptotique numérique*. Hermes Lavoissier.

Del Piero, G., G. Lancioni, and R. March (2007). "A variational model for fracture mechanics: numerical experiments". In: *Journal of the Mechanics and Physics of Solids* 55, pp. 2513-2537.

Fedelich, B. and A. Ehrlacher (1997). "An analysis of stability of equilibrium and of quasi-static transformations on the basis of the dissipation function". In: *European Journal of Mechanics A* 16 (5), pp. 833-855.

Francfort, G. and J.-J. Marigo (1998). "Revisiting brittle fracture as an energy minimization problem". In: *Journal of the Mechanics and Physics of Solids* 46.8, pp. 1319-1342.
Freddi, F. and G. Royer-Carfagni (2010). "Regularized variational theories of fracture : a unified approach". In: *Journal of the Mechanics and Physics of Solids* 58, pp. 1154–1174.

Halphen, B. and Q. S. Nguyen (1975). "Generalized Standard Materials". In: *Journal de Mécanique* 14.1, pp. 39–63.

Haslinger, J., V. Janovský, R. Kučera, and K. Motyčková (2018). "Nonsmooth continuation of parameter dependent static contact problems with Coulomb friction". In: *Mathematics and Computers in Simulation* 145, pp. 62–78.

Hernandez, V., J. E. Roman, and V. Vidal (2005). "SLEPc: A Scalable and Flexible Toolkit for the Solution of Eigenvalue Problems". In: *ACM Transactions on Mathematical Software* 31.3, pp. 351–362.

Hill, R. (1957). "On uniqueness and stability in the theory of finite elastic strain". In: *Journal of the Mechanics and Physics of Solids* 5.4, pp. 229–241.

Hutchinson, J. (1974). "Plastic Bukling". In: *Advances in Applied Mechanics* 14, pp. 67–144.

Lazarus, A., C. Maurini, and S. Neukirch (2015). "Stability of discretized nonlinear elastic systems". In: *Extremely Deformable Structures*. Ed. by D. Bigoni. Vienna: Springer Vienna, pp. 1–53.

Leon Baldelli, A., J.-F. Babadjian, B. Bourdin, D. Henao, and C. Maurini (2014). "A variational model for fracture and debonding of thin films under in-plane loadings". In: *Journal of the Mechanics and Physics of Solids* 70.0, pp. 320–348.

Leon Baldelli, A., B. Bourdin, J.-J. Marigo, and C. Maurini (2013). "Fracture and debonding of a thin film on a stiff substrate: Analytical and numerical solutions of a one-dimensional variational model". In: *Continuum Mechanics and Thermodynamics* 25.2-4, pp. 243–268.

Li, B., C. Peco, D. Millán, I. Arias, and M. Arroyo (2015). "Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy". In: *International Journal for Numerical Methods in Engineering* 102.3-4, pp. 711–727.

Lorentz, E. and V. Godard (2011). "Gradient damage models: Toward full-scale computations". In: *Computer Methods in Applied Mechanics and Engineering* 200.21–22, pp. 1927–1944.

Marigo, J.-J., C. Maurini, and K. Pham (2016). "An overview of the modelling of fracture by gradient damage models". en. In: *Meccanica* 51.12, pp. 3107–3128.

Mielke, A. and T. Roubíček (2015). *Rate-Independent Systems: Theory and Application*. Springer, p. 660.

Nguyen, Q. S. (1987). "Bifurcation and post-bifurcation analysis in plasticity and brittle fracture". In: *Journal of the Mechanics and Physics of Solids* 35.3, pp. 303–324.

— (1994). "Bifurcation and Stability in Dissipative Media (Plasticity, Friction, Fracture)". In: *Applied Mechanics Reviews* 47.1, pp. 1–31.

Nguyen, Q. S. (2000). *Stability and nonlinear solid mechanics*. John Wiley.

Nocedal, J. and S. Wright (2006). *Numerical Optimization*. 2nd ed. Springer series in operations research. Springer.

Petryk, H. (1993). "Theory of Bifurcation and Instability in Time-Independent Plasticity". In: *Bifurcation and Stability of Dissipative Systems*. Ed. by Q. S. Nguyen. Vienna: Springer Vienna, pp. 95–152.

Petryk, H. and K. Thermann (1992). "On discretized plasticity problems with bifurcations". In: *International Journal of Solids and Structures* 29.6, pp. 745–765.

Pham, K. and J.-J. Marigo (2012). "Damage localization and rupture with gradient damage models". In: *Frattura ed Integrita Strutturale* 19, pp. 5–19.
Pham, K., J.-J. Marigo, and C. Maurini (2011a). "The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models". In: *Journal of the Mechanics and Physics of Solids* 59.6, pp. 1163-1190.

Pham, K., H. Amor, J.-J. Marigo, and C. Maurini (2011b). "Gradient damage models and their use to approximate brittle fracture". In: *International Journal of Damage Mechanics* 20.4, pp. 618-652.

Pham, K. and J.-J. Marigo (2010). "Approche variationnelle de l'endommagement : II. Les modèles à gradient". In: *Comptes Rendus Mecanique* 338.4, pp. 199-206.

Potier-Ferry, M. (1985). "Towards a catastrophe theory for the mechanics of plasticity and fracture". In: *International Journal of Engineering Science* 23.8, pp. 821-837.

Roman, J. E., C. Campos, E. Romero, and A. Tomas (2020). *SLEPc Users Manual*. Tech. rep. DSIC-II/24/02 - Revision 3.13. D. Sistemes Informàtics i Computació, Universitat Politècnica de València.

Shanley, F. R. (1947). "Inelastic column theory". In: *Journal of the Aeronautical Sciences* 14.5.

Sicsic, P., J.-J. Marigo, and C. Maurini (2014). "Initiation of a periodic array of cracks in the thermal shock problem: A gradient damage modeling". In: *Journal of the Mechanics and Physics of Solids* 63.0, pp. 256-284.

Tanné, E., T. Li, B. Bourdin, J.-J. Marigo, and C. Maurini (2018). "Crack nucleation in variational phase-field models of brittle fracture". In: *Journal of the Mechanics and Physics of Solids* 110, pp. 80-99.

Xia, Z. and J. Hutchinson (2000). "Crack patterns in thin films". In: *Journal of the Mechanics and Physics of Solids* 48, pp. 1107-1131.