# On the formulation and implementation of mixed mode I and mode II extrinsic cohesive zone models with contact and friction

N.A. Collins-Craft<sup>a,b</sup> and V. Acary<sup>a</sup>

<sup>a</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000, Grenoble, France

<sup>b</sup>Laboratoire Navier, Ecole Nationale des Ponts et Chaussées, Institut Polytechnique de Paris, Université Gustave Eiffel, CNRS, Champs-sur-Marne, France

August 22, 2025

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Formulation of extrinsic cohesive zone models with contact and friction</td><td>4</td></tr><tr><td>2.1</td><td>State variables, powers and principle of virtual power</td><td>4</td></tr><tr><td>2.2</td><td>A nonsmooth thermo-mechanics potential</td><td>5</td></tr><tr><td>2.3</td><td>A linear evolution of the cohesion: triangle law</td><td>9</td></tr><tr><td>3</td><td>Nonsmooth elasto-dynamics of finite-dimensional systems</td><td>12</td></tr><tr><td>3.1</td><td>Finite-dimensional systems via space-discretisation</td><td>12</td></tr><tr><td>3.2</td><td>Nonsmooth dynamics and impacts</td><td>12</td></tr><tr><td>4</td><td>Numerical time integration</td><td>14</td></tr><tr><td>4.1</td><td>Principles of the time integration scheme</td><td>14</td></tr><tr><td>4.2</td><td>The discrete linear complementarity problem</td><td>15</td></tr><tr><td>4.3</td><td>Existence of the solution of the discrete LCP</td><td>17</td></tr><tr><td>4.4</td><td>Discrete energy balance</td><td>19</td></tr><tr><td>4.5</td><td>Quasi-static scheme</td><td>20</td></tr><tr><td>5</td><td>Numerical simulations</td><td>21</td></tr><tr><td>5.1</td><td>Quasi-static scalar case with elastic springs</td><td>22</td></tr><tr><td>5.1.1</td><td>The soft spring system</td><td>22</td></tr><tr><td>5.1.2</td><td>The stiff spring system</td><td>23</td></tr><tr><td>5.2</td><td>Dynamic scalar case with elastic springs</td><td>24</td></tr><tr><td>5.2.1</td><td>The soft spring system</td><td>24</td></tr><tr><td>5.2.2</td><td>The stiff spring system</td><td>25</td></tr><tr><td>5.3</td><td>The edge-cracked block</td><td>26</td></tr><tr><td>6</td><td>Conclusions</td><td>31</td></tr><tr><td>A</td><td>Reformulation of a bi-dimensional Coulomb-like inclusion into a complementarity problem</td><td>31</td></tr><tr><td>B</td><td>Linear complementarity problem development</td><td>33</td></tr><tr><td>B.1</td><td>Dynamic linear complementarity problem development</td><td>33</td></tr><tr><td>B.2</td><td>Quasi-static linear complementarity problem development</td><td>35</td></tr><tr><td>C</td><td>Finite element implementation details</td><td>36</td></tr></table>

**Notation** For vectors and tensors, we use the following notation:

$$ \|x\|^2 = \|x\|^2 = \underbrace{x_i x^i}_{\text{indicial notation}} = \underbrace{x \cdot x}_{\text{tensor notation}} = \underbrace{x^\top x}_{\text{vector notation}}. \quad (1) $$
# Abstract

An extrinsic cohesive zone model for mixed mode I and mode II fracture that encompasses contact and Coulomb friction is developed in the framework of nonsmooth mechanics. The model is extended to include the effects of dynamics with impact and sliding, and is discretised in time so that it can be written as a linear complementarity problem (LCP). The LCP is proved to have a solution, subject to a condition on the size of the time-step. Finally, we study the behaviour of the LCP system numerically, by observing the response of a simple test geometry to rapid loading, and observe the numerical method reproduces analytically predicted and experimentally observed behaviours, without requiring impractically small time-steps.

# 1 Introduction

Fracture mechanics is of wide interest in both theoretical and applied mechanics, as it plays an important role in phenomena ranging from geological scales such as earthquake rupture (Ferry et al., 2025; Okubo et al., 2019) and avalanche initiation (Bergfeld et al., 2021), to microscopic scales such as the behaviour of grains of sandstone (Jiang et al., 2021) under mechanical loading or metal alloys (Auth et al., 2022) in high temperature environments.

The classical description of linear elastic fracture mechanics (LEFM) was given by Griffith (1921), and while this model remains fundamental to how fracture is understood today, it suffers from significant problems. Most prominently, the stress field diverges at the crack tip, rendering the model nonphysical in an area that is referred to as the “fracture process zone”. Cohesive zone models (CZMs) were introduced by Dugdale (1960) and Barenblatt (1962), which served to regularise LEFM models by ensuring finite tractions across the incipient crack surface. The CZMs are described in terms of a traction-separation law, and the total area under the graph of this law is equivalent to the critical fracture energy $G_c$ of LEFM. It is convenient to describe the evolution of the interface in terms of the intensity of the cohesion, which we will denote $\beta \in [0, 1]$. In the following, we will take $\beta = 1$ to indicate a perfectly intact surface, while $\beta = 0$ represents a completely broken interface.

Cohesive zones are extremely difficult to observe experimentally, due to their generally small size (that can be on the order of nanometres (Azab et al., 2020)) and the extremely rapid motion that is characteristic of fracture processes under all but the most idealised conditions. Careful experimentation has allowed some direct observations to be made (Célarié et al., 2003; Guilloteau et al., 1996) of classical fractures, while the discovery that frictional ruptures behave in ways that are analogous to cracks allows another potential source of insight into the phenomenology of cohesive zones (Berman et al., 2020; Gvirtzman and Fineberg, 2021). Given the difficulty that experimentalists have had in obtaining clear observations of the fracture process zone, and hence in calibrating CZMs, typically calibration must be performed using sophisticated back-analysis of substantial amounts of data observed at the structural level (Réthoré and Estevez, 2013; Vargas et al., 2020), and the particular values of the material parameters such as $G_c$, the critical traction $\sigma_c$ (at which softening commences) and the critical length $\delta_c$ (the point at which the graph of the traction-separation law is zero) can vary substantially, depending on the details of the chosen model. Further complications can arise when microstructural heterogeneity is present, which will change the deformation of the crack front (Lebihain et al., 2022; Roch et al., 2023) and possibly require a larger set of parameters in order to be accurately modelled.

Cohesive zone models all ultimately derive from one of two families, intrinsic or extrinsic cohesive zone models. Intrinsic models feature an initial elasticity, with the critical traction being obtained when the separation reaches the hardening length $\delta_h$ (Falk et al., 2001). These elements are typically inserted into the mesh prior to simulation (whence their name, as they are intrinsic to the mesh), which makes code parallelisation relatively straightforward (Nguyen, 2014). However, the cohesive elements act like nonlinear springs due to their initial elastic branch, inducing a softening of the apparent stiffness of the structure. This phenomenon is referred to as artificial compliance, and it becomes more severe with increasing numbers of cohesive elements (that is to say, the more we refine the mesh, the worse our results) (Falk et al., 2001). This difficulty can be somewhat alleviated by increasing the initial elasticity of the cohesive law, but this in turn imposes unwanted costs on our simulation. For quasi-static analyses (that typically use implicit integration), the stiffness results in very unstable numerical integrations, while in dynamic analyses (that typically use explicit integration), the stiffness imposes such severe restrictions on the stable time-step size that the method becomes essentially unusable (Nguyen, 2014). In both cases, having an elasticity associated with the interface also allows for nonphysical results such as interface interpenetration, and hence negative displacement jumps, in the case where the two surfaces are pushed back into contact (Acary and Monerie, 2006). Finally, it is not actually possible to measure the stiffness across an interface before that interface exists, and as such whatever stiffness is determined is necessarily an arbitrary numerical property rather than a real physical parameter, as is hinted at by the ease with which modellers are willing to modify it to avoid artificial compliance.

Extrinsic models are characterised by an initial rigidity, and thus start to decohere only once the critical traction across the surface is exceeded (Kubair and Geubelle, 2003; Seagraves and Radovitzky, 2009). Due to their rigidity, this family of models do not effect the elasticity of the bulk, but they are typically inserted into the simulation on-the-fly (whence their name, as they are extrinsic to the mesh) (Zhou and Molinari, 2004). Given the absence of artificial compliance, extrinsic models are typically used for dynamic analyses (Dureisseix et al., 2024), but great care must be taken to ensure the time-continuity of the system. Failure to do so typically results in spurious oscillations of the system or an inability to converge in time (Papoulia et al.,
2003; Sam et al., 2005), but sufficiently sophisticated formulations are able to guarantee correct behaviour (Cazes et al., 2013). While extrinsic models do not demonstrate the pathological mesh-dependency that is characteristic of intrinsic models, the large number of elements and remeshing required to fully resolve the system for arbitrary crack paths is generally impractical, and so simulations are often not fully converged. The use of more sophisticated finite element techniques such as extended finite element (Moës and Belytschko, 2002) or discontinuous Galerkin (Versino et al., 2015) may lead to truly mesh-independent results, but come with their own implementation difficulties.

Regardless of how highly-refined the mesh is, cohesive zone models are still subject to difficulties related to ill-posedness (Foulk, 2010). In particular, in quasi-static simulations, solution jumps can appear in “soft” systems where the system stores more energy elastically than the cohesive zone is able to release by fully decohering (Acary and Monerie, 2006). In order to try and maintain well-posedness, researchers have either turned to sophisticated finite element schemes (Samimi et al., 2011) or viscous regularisation (Chaboche et al., 2001), however in our previous work (Collins-Craft et al., 2022) we were able to show that for pure mode I fracture and contact problems, working in dynamics (with a small enough time-step) is sufficient to maintain a well-posed problem.

While (carefully implemented) extrinsic models demonstrate desirable behaviour under monotonic loading, in more realistic loading cases such as those generated by multiple impacts or complex stress waves, cohesive elements may experience partial decohesion before being subject to unloading and reloading cycles. Almost all existing extrinsic cohesive zone models possess an elastic unload-reload branch (e.g. see Bybordiani and Dias-da-Costa (2021), Camacho and Ortiz (1996), Parrinello (2020), Parrinello and Borino (2020), Sam et al. (2005), and Snozzi and Molinari (2013)), and such models can be described as being “shifted intrinsic models”, as they have the same underlying mathematical structure as intrinsic models (Kubair and Geubelle, 2003). These shifted intrinsic models can demonstrate all the same pathologies as classical intrinsic models, as the elasticity of the unload-reload branch induces the same artificial compliance, and allows the same nonphysical behaviour such as interface interpenetration, and can have an arbitrarily large value for small amounts of decohesion. As such, there is a strong motivation to formulate extrinsic cohesive zone models that completely eliminate this unload-reload elasticity, and guarantee physically correct behaviour regardless of the loading history.

The correct mathematical framework in which to develop such extrinsic cohesive zone models is convex analysis, and in particular we follow the works of Jean-Jacques Moreau, who both developed and applied this mathematical framework to create nonsmooth mechanics (Moreau, 1970, 1974, 1986). This formulation requires a careful specification of the energy and pseudopotential of dissipation using indicator functions of convex sets, which are nondifferentiable functions. Using this framework allows the development of thermodynamically admissible evolution laws for the system that can include unilateral constraints on the internal variables and their rates (Halphen and Nguyen, 1975; Houlsby, 2019; Marigo, 1981). This formulation has been extended to cohesive zone modelling by Michel Frémond (Frémond, 1988, 2002, 2012a,b), which has in turn been extended in various ways for intrinsic cohesive zone models (Acary and Monerie, 2006; Chaboche et al., 2001; Monerie and Acary, 2001; Nkoumbou Kaptchouang et al., 2021; Perales et al., 2010; Raous et al., 1999), but which has only been minimally exploited for the development of extrinsic cohesive zone models (Jean et al., 2001; Talon and Curnier, 2003).

In our previous work (Collins-Craft et al., 2022), we used the nonsmooth mechanics framework to specify a mode I (opening mode) cohesive zone model that also included contact. Using this formulation, the model could be straightforwardly constrained to only admit physical solutions (i.e. the cohesion variable $\beta \in [0, 1]$ is strictly enforced, interface interpenetration is strictly forbidden), and in addition succeeded in eliminating the unload-reload elasticity from the formulation. The problem was then reformulated in terms of dynamics, and the discrete-in-space-and-time problem was able to be written as a linear complementarity problem. This problem formulation is particularly favourable, as it enabled a proof of the well-posedness of the problem, and is able to be numerically resolved in a very efficient manner, with comparatively large time-steps. The aim of this paper is to extend the mode I model to mixed modes I and II (opening and sliding modes) fracturing, and include contact and Coulomb friction in the formulation. The presence of tangential cohesion and friction substantially complicates the analysis and implementation of the model, and demands the exploitation of further mathematical results beyond those required for the mode I model.

**Novelty of the contribution and outline of the article.** The novelty of our work is that we formulate an extrinsic cohesive zone model that:

1. is based on nonsmooth thermo-mechanics principles, takes into account inertia and the inequality constraints on state variables and their rates such as unilateral contact and irreversibility and has a straightforward unload-reload behaviour that avoids any of the problems of shifted intrinsic model structures, and

2. accounts for both pure mode I and pure mode II fracture as well as mixed mode fracture without any mathematical singularities when the displacement jumps are zero, and

3. includes the contact and friction problem within the same mathematical framework as the cohesive zone problem, and a numerical algorithm that benefits from the following properties:
1. an implicit time-stepping scheme that is consistent with the nonsmooth contact dynamics approach that uses the Moreau-Jean scheme, that is dissipative in discrete time in the presence of friction and impacts, and energy-preserving in their absence, and which provides a stable numerical scheme at reasonably large time-steps, and

2. a linear complementarity problem formulation for the space-and-time-discretised problem with a proof of the existence of the solution that solves all of the constraints in an implicit manner, that in practice avoids solution jumps that occur in quasi-statics for models that are not regularised by viscosity or higher-order deformations, and

3. a formulation as a monolithic complementarity problem that allows the exploitation of efficient and robust algorithms developed by the mathematical programming community for this class of problems.

Finally, we demonstrate the practical interest of our approach by applying the model to some pertinent example systems and observe that the results match the expected physical behaviour.

## 2 Formulation of extrinsic cohesive zone models with contact and friction

First, we define the state variables of our model, and the corresponding powers associated with them. Then, using the principle of virtual power, we express the equilibrium equations and boundary conditions for our system. We then specify a particular constitutive model, study the continuous-time energy balance, and demonstrate its analytical solution for a single contact point.

### 2.1 State variables, powers and principle of virtual power

We consider a body $\Omega \in \mathbb{R}^3$. The current position is defined by the vector $\mathbf{x}$ and the initial position by $\mathbf{X}$. We consider that the mechanical properties and the geometry of the body do not vary with $x_3$ (i.e. we are in the plane strain or plane stress condition), and so we consider only the behaviour of the system with respect to $\mathbf{x} = [x_1, x_2]$ and treat the system as if it were two dimensional i.e. $\mathbf{x} \in \mathbb{R}^2$. We then define a vector of the displacement $\mathbf{u}(\mathbf{x}) = \mathbf{x} - \mathbf{X}$, from which we obtain the velocity $\mathbf{v}(\mathbf{x}) = \dot{\mathbf{u}}(\mathbf{x})$. We consider that the body is initially undamaged, and hence at a point $\mathbf{x} \in \Omega$, the displacement $\mathbf{u}(\mathbf{x})$ and the velocity $\mathbf{v}(\mathbf{x})$ are continuously differentiable functions of $\mathbf{x}$. When an interface is created by the process of fracture, two material points $x_l$ and $x_r$ are defined by splitting the bodies assuming that they correspond to the material point $\mathbf{X}$ initially. For ease of notation, we denote $x_l$ by $\mathbf{x}$. The displacement jump that characterises the crack is defined by the difference in the position of the material point that was at $\mathbf{X}$ initially, with respect to $\mathbf{x}$ and $x_r$, that is $[\mathbf{u}(\mathbf{x}, x_r)] = \mathbf{u}(x_r) - \mathbf{u}(\mathbf{x})$. We define an orthonormal local frame on the interface defined by $(\mathbf{x}, \mathbf{n}, t)$ where $\mathbf{n} \in \mathbb{R}^2$ is the normal unit vector from $\mathbf{x}$ towards $x_r$ and the vector $t \in \mathbb{R}^2$ gives the tangential direction along the fracture surface. In this work we consider the normal displacement (jump), defined by $u_N(\mathbf{x}, x_r) = [\mathbf{u}(\mathbf{x}, x_r)] \cdot \mathbf{n} \in \mathbb{R}$, and the tangential displacement (jump) in the along-crack direction, defined by $u_T(\mathbf{x}, x_r) = [\mathbf{u}(\mathbf{x}, x_r)] \cdot t \in \mathbb{R}$. The relative normal velocity (jump in space) is given by $v_N = \dot{u}_N$, while the tangential velocity (jump in space) is given by $v_T = \dot{u}_T$.

To describe the state of the cohesion we introduce the cohesion variable $\beta(\mathbf{x}) \in [0, 1]$, using the notation introduced by Frémond for describing the intensity of cohesion. For a point $\mathbf{x}$ on the interface, the power of the cohesion for a surface $\Gamma$ is defined by

$$ \mathcal{P}_{\text{coh}} = \int_{\Gamma} \dot{\beta} A \, dx, \qquad (2) $$

where we have introduced $A$ which is the thermodynamic dual force (driving force) associated with $\beta$. Now, over the same surface, the power of contact is given by

$$ \mathcal{P}_{\text{con}} = \int_{\Gamma} (v_N r_N + v_T r_T) \, dx, \qquad (3) $$

where $r_N$ is the normal reaction force related to the stress $\sigma(\mathbf{x})$ at the interface by $r_N = -\sigma \cdot \mathbf{n} \cdot \mathbf{n}$, and $r_T$ is the corresponding tangential reaction force $r_T = -\sigma \cdot \mathbf{n} \cdot \mathbf{t}$.

For the material in $\Omega$, the power of the external, internal and inertial forces are respectively given by

$$ \mathcal{P}_{\text{ext}} = \int_{\Omega} \mathbf{v} \cdot \mathbf{f} \, dx - \int_{\Gamma_N} \mathbf{v} \cdot \boldsymbol{\tau} \, dx - \int_{\Gamma} \dot{\beta} \Theta \, dx, \qquad (4) $$

$$ \mathcal{P}_{\text{int}} = - \int_{\Omega} \boldsymbol{\sigma} : \dot{\boldsymbol{\varepsilon}} \, dx + \int_{\Gamma} v_N r_N \, dx + \int_{\Gamma} v_T r_T \, dx + \int_{\Gamma} \dot{\beta} A \, dx, \qquad (5) $$

$$ \mathcal{P}_{\text{acc}} = \int_{\Omega} \rho \mathbf{v} \cdot \dot{\mathbf{v}} \, dx, \qquad (6) $$

where $\mathbf{f}$ is the body force in $\Omega$, $\boldsymbol{\tau}$ is the surface traction on $\Gamma_N$ (i.e. the region of the surface where the Neumann boundary condition is applied), $\Theta$ is an external force that does work on the cohesion (such forces can occur due to thermal or chemical
effects) that we take to be identically zero in this work, $\varepsilon$ is the strain in $\Omega$, $\rho$ is the density and $\dot{v}$ is the acceleration.

The principle of virtual power states that for any virtual velocities $\bar{v}$, $\dot{\bar{\varepsilon}}$ and $\dot{\bar{\beta}}$, we have

$$ \bar{P}_{\text{acc}} = \bar{P}_{\text{ext}} + \bar{P}_{\text{int}}, $$

$$ \int_{\Omega} \rho \bar{v} \cdot \dot{v} \, dx = \int_{\Omega} \bar{v} \cdot f \, dx - \int_{\Gamma_{\mathcal{N}}} \bar{v} \cdot \tau \, dx - \int_{\Gamma} \dot{\bar{\beta}} \Theta \, dx - \int_{\Omega} \boldsymbol{\sigma} : \dot{\bar{\varepsilon}} \, dx + \int_{\Gamma} (\bar{v}_{\mathcal{N}} r_{\mathcal{N}} + \bar{v}_{\mathcal{T}} r_{\mathcal{T}}) \, dx + \int_{\Gamma} \dot{\bar{\beta}} A \, dx. \quad (7) $$

For further details, a rigorous mathematical treatment of this principle may be found in Frémond (1988). Given sufficient smoothness assumptions on the fields, the equations describing the equilibrium and boundary conditions of the system are given by localisation as:

$$ \left\{ \begin{array}{ll} \nabla \cdot \boldsymbol{\sigma} + f = \rho \dot{v} & \text{in } \Omega, \\ A = \Theta = 0 & \text{on } \Gamma, \\ \boldsymbol{\tau} = \boldsymbol{\sigma} \cdot \boldsymbol{n} & \text{on } \Gamma_{\mathcal{N}}, \\ u = u_{\mathcal{D}} & \text{on } \Gamma_{\mathcal{D}}, \\ r_{\mathcal{N}} = -\boldsymbol{\sigma} \cdot \boldsymbol{n} \cdot \boldsymbol{n} & \text{on } \Gamma, \\ r_{\mathcal{T}} = -\boldsymbol{\sigma} \cdot \boldsymbol{n} \cdot \boldsymbol{t} & \text{on } \Gamma. \end{array} \right. \quad (8) $$

where $\Gamma_{\mathcal{D}}$ is the region of the boundary on which the Dirichlet boundary condition is applied and $u_{\mathcal{D}}$ is the prescribed displacement.

## 2.2 A nonsmooth thermo-mechanics potential

In defining our model, we begin by defining a ratio between the critical shear traction to the critical opening traction. Then, we write the free energy potential for the bulk and surface, from which we derive the state laws. By specifying a dissipation pseudo-potential, we obtain the laws describing the irreversible behaviour, and thus we are able to describe the complete model as a complementarity problem. Finally, we then describe the energy balance of the system, from which we obtain an expression for the fracture energy.

**Critical traction ratio.** We start by firstly introducing the ratio of the critical shear traction to the critical opening traction $\gamma = \frac{\sigma_{c,\Pi}}{\sigma_{c,i}}$ as a convenient expression of the anisotropy of the cohesive strength (Xu and Needleman, 1994). From hereon, we denote $\sigma_{c,i}$ by $\sigma_c$ for the sake of notational simplicity.

**Free energy and reversible state laws.** We consider the free energy of the system, from which we obtain the laws describing the reversible behaviour of the system. The total free energy $\Psi$ of the system is the sum of the free energy in the bulk with the free energy of the surface:

$$ \Psi = \int_{\Omega} \Psi_{e}(\varepsilon) \, dx + \int_{\Gamma} \Psi_{s}(u_{\mathcal{N}}, u_{\mathcal{T}}, \beta) \, dx, \quad (9) $$

where $\Psi_e$ and $\Psi_s$ are the volume and surface free energies, respectively. As we consider an elasto-brittle system in this work, we will assume that all strain is elastic. Firstly, the stresses may be obtained by assuming a classical linear elastic potential for the bulk:

$$ \Psi_e(\varepsilon) = \frac{1}{2} \varepsilon : E : \varepsilon, \quad (10) $$

$$ \boldsymbol{\sigma}(\varepsilon) = \frac{\partial \Psi_e(\varepsilon)}{\partial \varepsilon} = E : \varepsilon, \quad (11) $$

where $E$ is a fourth order stiffness tensor. The stress-like variables for the surface are similarly derived from the surface potential by

$$ \left\{ \begin{array}{l} -r_{\mathcal{N}}^{r} \in \partial u_{\mathcal{N}} \Psi_s(u_{\mathcal{N}}, u_{\mathcal{T}}, \beta), \\ -r_{\mathcal{T}}^{r} \in \partial u_{\mathcal{T}} \Psi_s(u_{\mathcal{N}}, u_{\mathcal{T}}, \beta), \\ -A^{r} \in \partial \beta \Psi_s(u_{\mathcal{N}}, u_{\mathcal{T}}, \beta), \end{array} \right. \quad (12) $$

where $r_{\mathcal{N}}^r$ is the reversible part of the normal reaction force, $r_{\mathcal{T}}^r$ is the reversible part of the tangential reaction force, and $\partial u_{\mathcal{N}}$, $\partial u_{\mathcal{T}}$ and $\partial \beta$ indicate the subdifferentials with respect to $u_{\mathcal{N}}$, $u_{\mathcal{T}}$ and $\beta$ of a convex (in each individual variable) but nonsmooth potential. The first assumption in our model is that the normal displacement is constrained to be positive, i.e. $u_{\mathcal{N}} \ge 0$, which
is enforced as a unilateral constraint, while the cohesion variable is constrained to be $0 \le \beta \le 1$. Under these constraints, a rather generic form of the surface free energy may be specified by

$$ \Psi_s(u_N, u_T, \beta) = \psi(u_N, u_T, \beta) + \mathcal{I}_{\mathbb{R}^+}(u_N) + \mathcal{I}_{[0,1]}(\beta), \quad (13) $$

where $\mathcal{I}_C$ is the indicator function of a convex set $C$. In order to obtain an extrinsic CZM, the model must have no tangent stiffness when $u_N \ge 0$ or $|u_T| \ge 0$. This constrains the possible forms of (13) to those functions that fulfil the conditions

$$ \frac{\partial^2 \psi}{\partial u_N^2} = 0 \text{ and } \left. \frac{\partial^2 \psi}{\partial u_T^2} \right|_{u_T \neq 0} = 0. \text{ One such free energy is given by} $$

$$ \Psi_s(u_N, u_T, \beta) = \beta \sigma_c u_N + \beta \sigma_c \gamma |u_T| + w f(\beta) + \mathcal{I}_{\mathbb{R}^+}(u_N) + \mathcal{I}_{[0,1]}(\beta), \quad (14) $$

where $w > 0$ is the surface free energy which is released by decohesion and $f(\beta)$ is a function that parametrises the evolution of $\beta$ as decohesion progresses. The state laws in their specific form are then obtained by applying (12) to (14):

$$ \left\{ \begin{array}{l} -(r_N^r + \beta \sigma_c) \in \partial \mathcal{I}_{\mathbb{R}^+}(u_N), \\ -r_T^r \in \beta \sigma_c \gamma \partial |u_T|, \\ -(A^r + \sigma_c u_N + \sigma_c \gamma |u_T| + w f'(\beta)) \in \partial \mathcal{I}_{[0,1]}(\beta). \end{array} \right. \quad (15) $$

The surface free energy $\Psi_s$ is not in general a convex function of its arguments $(u_N, u_T, \beta)$. However, $\Psi_s$ is separately convex with respect to $u_N$ and $u_T$. So long as we respect the condition $f''(\beta) \ge 0$, $\Psi_s$ will remain convex with respect to $\beta$. However, as we show in Remark 2, this does not imply convexity of $\Psi_s$. The continuous part of the free energy in (14) is composed of three terms. The first two terms $\beta \sigma_c u_N$ and $\beta \sigma_c \gamma |u_T|$ account for the potential energy of the normal and tangential cohesive forces $\beta \sigma_c$ and $\beta \sigma_c \gamma$ in the displacement field given by $u_N$ and $u_T$. The third term $w f(\beta)$ accounts for the surface free energy released by decohesion at a given $\beta$. We also place certain constraints on the form of the function $f$. When the interface is intact and $\beta = 1$, we require that $f(1) = 0$, as there cannot be any free energy on a surface that does not yet exist. When the interface is fully broken and $\beta = 0$, the cohesive free energy $w$ must have been completely released, and hence $f(0) = 1$. Finally, the release of the cohesive free energy with $\beta$ must be monotonic, and hence $f'(\beta) \le 0$. Given these conditions, and provided the minimum of $f$ is attained at $\beta = 1$, $f$ will be suitable for our purposes. In this work, we propose a triangle cohesive law that fulfils these requirements, although other shapes that satisfy these requirements are also possible.

**Dissipation pseudo-potential and irreversible processes.** In order to respect the second law of thermodynamics, we define a dissipation function describing the irreversible part of the decohesion process under isothermal conditions:

$$ \mathcal{D} = -\mathcal{P}_{\text{int}} - \int_{\Omega} \dot{\Psi}_e(\varepsilon) \, dx - \int_{\Gamma} \dot{\Psi}_s(u_N, u_T, \beta). \quad (16) $$

This function must be non-negative for all admissible values of the state variables. In order to be able to compute the time derivatives of a nonsmooth and nonconvex pseudo-potential, we must make certain assumptions. The functions $u_N(t)$, $u_T(t)$ and $\beta(t)$ are assumed to be absolutely continuous, and hence they have a derivative almost everywhere but not necessarily at any given point. As absolutely continuous functions have bounded variations, the left and right derivatives exist for all of these functions. Assuming $\dot{\Psi}_s$ is convex in $u_N$, $u_T$ and in $\beta$ separately, and $\dot{\Psi}_s(u_N, u_T, \beta)$ is also an absolutely continuous function for $u_N \ge 0$, $|u_T| \ge 0$, $\beta \in [0, 1]$, a result in Frémond (2002, Appendix A.1.9) provides us with the following inequality:

$$ \dot{\Psi}_s(u_N, u_T, \beta) \le -v_N r_N^r - v_T r_T^r - \dot{\beta} A^r, \quad (17) $$

for any $r_N^r, r_T^r$ and $\dot{\beta}$ that satisfies the inclusions in (15). Now, substituting this inequality into the definition of the dissipation rate given by (16)

$$ \mathcal{D} \ge \int_{\Omega} \sigma : \dot{\varepsilon} \, dx - \int_{\Gamma} (v_N r_N + v_T r_T + \dot{\beta} A) \, dx - \int_{\Omega} \frac{\partial \Psi_e(\varepsilon)}{\partial \varepsilon} \dot{\varepsilon} \, dx + \int_{\Gamma} (v_N r_N^r + v_T r_T^r + \dot{\beta} A^r) \, dx. \quad (18) $$

In (18), the terms integrated over the bulk $\Omega$ cancel out. Considering only the surface $\Gamma$, we obtain

$$ \mathcal{D} \ge \int_{\Gamma} (-v_N r_N - v_T r_T - \dot{\beta} A + v_N r_N^r + v_T r_T^r + \dot{\beta} A^r) \, dx. \quad (19) $$

Since the laws of thermodynamics strictly require that $\mathcal{D} \ge 0$, we will assume that

$$ -v_N r_N - v_T r_T - \dot{\beta} A + v_N r_N^r + v_T r_T^r + \dot{\beta} A^r \ge 0, \text{ or equivalently, } -v_N r_N^{ir} + v_T r_T^{ir} - \dot{\beta} A^{ir} \ge 0, \quad (20) $$

where we have exploited the standard decompositions $r_N = r_N^r + r_N^{ir}$, $r_T = r_T^r + r_T^{ir}$ and $A = A^r + A^{ir}$ for the irreversible parts of the stress-like variables. We may guarantee that this inequality will be respected by specifying a proper closed convex pseudo-potential of dissipation $\Phi(v_N, v_T, \beta)$, from which the irreversible stress-like variables may be derived:

$$ \left\{ \begin{array}{l} -r_N^{ir} \in \partial_{v_N} \Phi(v_N, v_T, \dot{\beta}), \\ -r_T^{ir} \in \partial_{v_T} \Phi(v_N, v_T, \dot{\beta}), \\ -A^{ir} \in \partial_{\dot{\beta}} \Phi(v_N, v_T, \dot{\beta}). \end{array} \right. \quad (21) $$
For reasons of simplicity, we will assume that the dissipation process of fracture depends exclusively on the rate of $\beta$, and that it does so in a linear way, resulting in a rate-independent dissipative behaviour. However, we also wish to include the effects of frictional contact in our model. To do so, we specify the pseudo-potential of dissipation as

$$ \Phi(v_N, v_T, \dot{\beta}) = \mathcal{I}_{\mathbb{R}^-}(\dot{\beta}) + \mu(r_N + \beta\sigma_c)|v_T|, \qquad (22) $$

where the last term is the standard pseudo-potential associated with the Coulomb friction and $\mu$ is the dynamic coefficient of friction. We assume that the static and dynamic coefficients of friction are identical. This model of dissipation imposes that the evolution of $\beta$ must decrease with time, i.e. $\dot{\beta} \le 0$. The dissipative laws are thus

$$ \left\{ \begin{array}{l} -r_N^{\text{ir}} = 0, \\ -r_T^{\text{ir}} \in \mu(r_N + \beta\sigma_c)\partial|v_T|, \\ -A^{\text{ir}} \in \partial\mathcal{I}_{\mathbb{R}^-}(\dot{\beta}). \end{array} \right. \qquad (23) $$

We may also write the last line of (23) as

$$ \dot{\beta} \in \partial \mathcal{I}_{\mathbb{R}^+}(-A^{\text{ir}}), \qquad (24) $$

which allows us to interpret that $A^{\text{ir}}$ is the force that drives the evolution of $\dot{\beta}$. In (8), $A = 0$, so as a consequence $A^r = -A^{\text{ir}}$, and hence we may conclude

$$ \dot{\beta} \in \partial \mathcal{I}_{\mathbb{R}^+}(A^r). \qquad (25) $$

**Remark 1.** The dissipation pseudo-potential in (22) can be straightforwardly generalised to generate models that are rate-dependent for the decohesion process, the frictional sliding, or both, by making the dissipation nonlinear in the relevant dissipative variable (see Acary and Monerie (2006) for further details).

**Remark 2 (Nonconvexity of $\Psi_s$).** Even if we consider that $f(\beta)$ is convex, $\Psi_s$ is not a priori a convex function. Let us consider here that $\beta \in (0, 1)$ and $u_N \ge 0$. The surface free energy reduces to

$$ \Psi_s(u_N, u_T, \beta) = \psi(u_N, u_T, \beta) = \beta\sigma_c u_N + \beta\sigma_c\gamma|u_T| + w f(\beta), \qquad (26) $$

from (15) and the assumption that $\beta \in (0, 1)$. Let us consider that we have an evolution with $\dot{\beta} < 0$, then $-A^{\text{ir}} = A^r = 0$ and from (30), we get

$$ \sigma_c u_N + \sigma_c \gamma |u_T| = -w f'(\beta). \qquad (27) $$

which is the equation that gives the decohesion curve that relates $\beta$ to the displacement jump. In this case, the surface free energy reduces to

$$ \Psi_s(u_N, u_T, \beta) = w(f(\beta) - \beta f'(\beta)) = g(\beta). \qquad (28) $$

Computing the second derivative of the function $g(\beta)$ yields

$$ g''(\beta) = -w(f''(\beta) + \beta f'''(\beta)). \qquad (29) $$

Assuming that $f$ is strictly convex ($f''(\beta) > 0$) and at least $C^3$, we can remark that $g''(\beta) < 0$ for small values of $\beta$. We conclude that $\Psi_s$ can be nonconvex while being separately convex in $\beta$.

**Complete extrinsic cohesive zone model.** We start by noting that since $r_N^{\text{ir}} = 0$, we must have $r_N = r_N^r$. Thus, the complete model of the interface is given by

$$ \left\{ \begin{array}{l} \dot{\beta} \in \partial \mathcal{I}_{\mathbb{R}^+}(A^r), \\ -(r_N + \beta\sigma_c) \in \partial \mathcal{I}_{\mathbb{R}^+}(u_N), \\ -r_T^r \in \beta\sigma_c\gamma\partial|u_T|, \\ -r_T^{\text{ir}} \in \mu(r_N + \beta\sigma_c)\partial|v_T|, \\ -(A^r + \sigma_c u_N + \sigma_c\gamma|u_T| + w f'(\beta)) \in \partial \mathcal{I}_{[0,1]}(\beta). \end{array} \right. \qquad (30) $$

**Formulation of the model using slack variables** Introducing the slack variables $\xi$, $\lambda$ and denoting the contact forces by $r_N^c = r_N + \beta\sigma_c$, we may write

$$ \left\{ \begin{array}{l} \dot{\beta} = -\lambda, \quad A^r + \sigma_c u_N + \sigma_c \gamma |u_T| + w f'(\beta) = \xi, \quad r_N^c = r_N + \beta \sigma_c, \quad -r_N^c \in \mathbb{N}_{\mathbb{R}^+}(u_N), \\ -\xi \in \mathbb{N}_{[0,1]}(\beta), \quad -\lambda \in \mathbb{N}_{\mathbb{R}^+}(A^r), \quad -u_T \in \mathbb{N}_{[-\beta \sigma_c \gamma, \beta \sigma_c \gamma]}(r_T^r), \quad -v_T \in \mathbb{N}_{[-\mu r_N^c, \mu r_N^c]}(r_T^{\text{ir}}). \end{array} \right. \qquad (31) $$

In our previous work (Collins-Craft et al., 2022), we showed that while the $\xi$ constraint is seemingly two-sided, one of the sides is redundant i.e. $\beta \le 1$ does not need to be enforced, only $\beta \ge 0$. Hence, we may in fact write $-\xi \in \mathbb{N}_{\mathbb{R}^+}(\beta)$, which simplifies the problem somewhat:

$$ \left\{ \begin{array}{l} \dot{\beta} = -\lambda, \quad A^r + \sigma_c u_N + \sigma_c \gamma |u_T| + w f'(\beta) = \xi, \quad r_N^c = r_N + \beta \sigma_c, \quad 0 \le r_N^c \perp u_N \ge 0, \\ 0 \le \xi \perp \beta \ge 0, \quad 0 \le \lambda \perp A^r \ge 0, \quad -u_T \in \mathbb{N}_{[-\beta \sigma_c \gamma, \beta \sigma_c \gamma]}(r_T^r), \quad -v_T \in \mathbb{N}_{[-\mu r_N^c, \mu r_N^c]}(r_T^{\text{ir}}). \end{array} \right. \qquad (32) $$
**Energy balance** With the chosen constitutive laws, the power of internal forces can be written as

$$ 
\begin{aligned}
\mathcal{P}_{\text{int}} &= -\int_{\Omega} \boldsymbol{\varepsilon} : \mathbf{E} : \dot{\boldsymbol{\varepsilon}} \, dx + \int_{\Gamma} r_{\text{N}} v_{\text{N}} \, dx + \int_{\Gamma} r_{\text{T}} v_{\text{T}} \, dx + \int_{\Gamma} A \dot{\boldsymbol{\beta}} \, dx, \\
&= -\frac{d}{dt} \left( \int_{\Omega} \boldsymbol{\varepsilon} : \mathbf{E} : \boldsymbol{\varepsilon} \, dx \right) + \int_{\Gamma} \left( r_{\text{N}} v_{\text{N}} + r_{\text{T}} v_{\text{T}} + A \dot{\boldsymbol{\beta}} \right) \, dx, \\
&= -\dot{\mathcal{U}} + \mathcal{P}_{\text{s,int}},
\end{aligned}
\quad (33) $$

where $\mathcal{U}$ is the potential elastic strain energy, and $\mathcal{P}_{\text{s,int}}$ the power of surface internal forces defined by $\mathcal{P}_{\text{s,int}} = \int_{\Gamma} p_{\text{s,int}} \, dx$ with $p_{\text{s,int}} = r_{\text{N}}v_{\text{N}} + r_{\text{T}}v_{\text{T}} + A\dot{\boldsymbol{\beta}}$. Now, let us expand $p_{\text{s,int}}$:

$$ 
\begin{aligned}
p_{\text{s,int}} &= r_{\text{N}}v_{\text{N}} + r_{\text{T}}v_{\text{T}} + A\dot{\boldsymbol{\beta}}, \\
&= v_{\text{N}}(r_{\text{N}}^{\text{ir}} + r_{\text{N}}^{\text{r}}) + v_{\text{T}}(r_{\text{T}}^{\text{ir}} + r_{\text{T}}^{\text{r}}) + \dot{\boldsymbol{\beta}}(A^{\text{ir}} + A^{\text{r}}), \\
&= v_{\text{N}}r_{\text{N}}^{\text{r}} + v_{\text{T}}r_{\text{T}}^{\text{r}} + \dot{\boldsymbol{\beta}}A^{\text{r}} + v_{\text{N}}r_{\text{N}}^{\text{ir}} + v_{\text{T}}r_{\text{T}}^{\text{ir}} + \dot{\boldsymbol{\beta}}A^{\text{ir}}, \\
&= v_{\text{N}}r_{\text{N}}^{\text{r}} + v_{\text{T}}r_{\text{T}}^{\text{r}} + \dot{\boldsymbol{\beta}}A^{\text{r}} + v_{\text{T}}r_{\text{T}}^{\text{ir}},
\end{aligned}
\quad (34) $$

where in order to pass from the second-last to last line we have exploited that $r_{\text{N}}^{\text{ir}} = 0$ and $\dot{\boldsymbol{\beta}}A^{\text{ir}} = 0$ using (23). The last term of (34) is the frictional power

$$ v_{\text{T}}r_{\text{T}}^{\text{ir}} = -\mu r_{\text{N}}^{c}|v_{\text{T}}| = -\mu(r_{\text{N}} + \beta\sigma_{c})|v_{\text{T}}|, \quad (35) $$

since $-r_{\text{T}}^{\text{ir}} \in \mu(r_{\text{N}} + \beta\sigma_{c}) \operatorname{sgn}(v_{\text{T}})$. The idea is now to relate the power of the surface internal forces, minus the frictional power, to the rate of the free energy. When the free energy is smooth, the following chain rule is used

$$ \dot{\psi} = v_{\text{N}} \frac{\partial \psi}{\partial u_{\text{N}}} + v_{\text{T}} \frac{\partial \psi}{\partial u_{\text{T}}} + \dot{\beta} \frac{\partial \psi}{\partial \dot{\beta}}. \quad (36) $$

Since in our case, the free energy is nonsmooth we will assume that the following chain rule is valid

$$ \dot{\psi} = v_{\text{N}} \partial_{u_{\text{N}}} \psi + v_{\text{T}} \partial_{u_{\text{T}}} \psi + \dot{\beta} \partial_{\beta} \psi, \quad (37) $$

almost everywhere. Demonstrating this rule would require some technical mathematical developments (see for instance Bolte and Pauwels (2021)) that are outside the scope of this article. Using the constitutive relation for the reversible processes in (12), we obtain

$$ \dot{\psi} = -v_{\text{N}}r_{\text{N}}^{\text{r}} - v_{\text{T}}r_{\text{T}}^{\text{r}} - \dot{\beta}A^{\text{r}}. \quad (38) $$

From (34) and (38), the relation between the power of surface internal forces and the rate of the free energy is given by

$$ p_{\text{s,int}} = r_{\text{N}}v_{\text{N}} + r_{\text{T}}v_{\text{T}} + A\dot{\boldsymbol{\beta}} = v_{\text{T}}r_{\text{T}}^{\text{ir}} - \dot{\psi} \quad \text{almost everywhere.} \quad (39) $$

Hence, we are able to conclude that the power of the surface internal forces is almost everywhere

$$ \mathcal{P}_{\text{s,int}} = \int_{\Gamma} (v_{\text{T}}r_{\text{T}}^{\text{ir}} - \dot{\psi}) \, dx. \quad (40) $$

Then, we can apply the principle of virtual power for the velocities of the system to write the balance equation for the kinetic energy $\mathcal{K}$, almost everywhere given by $\dot{\mathcal{K}} = \mathcal{P}_{\text{ext}} + \mathcal{P}_{\text{int}}$. Then, substituting in (33) and (40) we arrive at the expression

$$ \dot{\mathcal{K}} + \dot{\mathcal{U}} + \int_{\Gamma} (\dot{\psi} - v_{\text{T}}r_{\text{T}}^{\text{ir}}) \, dx = \mathcal{P}_{\text{ext}}. \quad (41) $$

We apply the first law of thermodynamics and assume isothermal conditions, which results in the energy balance $\dot{\mathcal{K}} + \dot{\mathcal{E}} = \mathcal{P}_{\text{ext}}$, where $\mathcal{E}$ is the internal energy. It is then straightforward to substitute the kinetic energy balance and (33) to arrive at

$$ \dot{\mathcal{E}} = -\mathcal{P}_{\text{int}} = \int_{\Gamma} (\dot{\psi} - v_{\text{T}}r_{\text{T}}^{\text{ir}}) \, dx + \dot{\mathcal{U}}. \quad (42) $$

Let us define the fracture energy as the opposite of the work done by the cohesion forces

$$ \mathcal{G} = \int \int_{\Gamma} \beta \sigma_{c} v_{\text{N}} - v_{\text{T}} r_{\text{T}}^{\text{r}} \, dx \, dt. \quad (43) $$
To relate the fracture energy (43) to the rate of change of the free energy $\dot{\psi}$, we expand further the equation (38) using the constitutive relations

$$ 
\begin{aligned}
\dot{\psi} &= -v_N r_N^r - v_T r_T^r - \dot{\beta} A^r, \\
&= -v_N (r_N^c - \beta \sigma_c) - v_T r_T^r + \dot{\beta} A^{ir}, \\
&= -v_N r_N^c + v_N \beta \sigma_c - v_T r_T^r,
\end{aligned}
\quad (44) $$

since $\dot{\beta}A^{ir} = 0$ using the last line of (23). Since $-r_N^c \in N_{R+}(u_N)$, it is also possible to show that $v_N r_N^c = 0$ almost everywhere and we conclude that $\dot{\psi} = v_N \beta \sigma_c - v_T r_T^r$ and

$$ \mathcal{G} = \iint_{\Gamma} \dot{\psi} \, dx \, dt \quad (45) $$

The dissipated energy by friction is denoted by

$$ \mathcal{F} = \iint_{\Gamma} -v_T r_T^{ir} \, dx \, dt. \quad (46) $$

as the opposite of the work done by the friction force. We can then write the incremental energy balance of the system as

$$ \Delta \mathcal{K} + \Delta \mathcal{U} + \Delta \mathcal{G} + \Delta \mathcal{F} = \int_{t_1}^{t_2} \mathcal{P}_{\text{ext}} \, dt, \quad \text{and} \quad \Delta \mathcal{E} = \Delta \mathcal{U} + \Delta \mathcal{G}. \quad (47) $$

**Remark 3.** The relation (45) allows us to conclude that to pass from an intact interface with $u_N = 0$, $u_T = 0$ and $\beta = 1$ to a broken interface with $\beta = 0$ on $\Gamma$, we must have $\Delta\mathcal{G} = \int_{\Gamma} w \, dx$.

## 2.3 A linear evolution of the cohesion: triangle law

In order to obtain a specific cohesive zone model from our general framework, we specify the exact form of the surface potential. This potential must fulfil certain conditions, namely that when the surface does not exist, i.e. $u_N = 0$, $u_T = 0$ and $\beta = 1$, $\Psi_s = 0$, and that when the surface is fully decohered, i.e. $\beta = 0$, $\Psi_s = w$. We assume that this decohesion is complete when the displacement jump has surpassed a critical length $\delta_c > 0$, which can be different in each direction. Choosing

$$ 
\begin{cases}
w = \frac{1}{2}\sigma_{c,i}\delta_{c,i} = \frac{1}{2}\sigma_{c,ii}\delta_{c,ii} = \frac{1}{2}\sigma_{c,i}\gamma\delta_{c,ii} = G_c, \\
f(\beta) = (\beta - 1)^2,
\end{cases}
\quad (48) $$

fulfils the required conditions for the energy while ensuring that the potential remains convex in $\beta$. As the fracture energy is constant in both directions and $\sigma_{c,i} = \sigma_{c,ii} = \gamma$, we also have $\delta_{c,ii} = \delta_{c,i}/\gamma$.

**Remark 4.** Our formulation for (48) is essentially determined by our conceptual framework, where we consider fracture in elasto-brittle systems to be a process of transforming bulk strain energy to surface energy and kinetic energy. As such, the energy of a surface cannot be determined by the type of fracture that formed it (i.e. mode I, II or III), and thus there can be only a single value of $G_c$ for the material. This does not prevent different values of the critical traction in each direction, but the critical opening distance must in turn adjust to maintain a constant $G_c$. Such a formulation is common with other authors that adopt a similar energetic picture to us e.g Talon and Curnier (2003) or Lorentz (2008). It is certainly possible to have different values of $G_c$ depending on the fracture mode, such as in the seminal paper of Camanho et al. (2003) or more recent works in the same framework (Venzal et al., 2020). Differing values of $G_c$ depending on the mode implies that significant plasticity is being captured by the cohesive zone model.

While $\beta \in [0, 1]$ and $\dot{\beta} < 0$, $A^r = 0$ and $\xi = 0$, meaning that after substituting the form (48) into the expression for $A^r$ in (32), we obtain $\beta = 1 - \frac{u_N + \gamma|u_T|}{\delta_{c,N}}$. Substituting the displacement-cohesion law into the free energy (14), we get:

$$ \Psi_s(u_N, u_T) = \sigma_c(u_N + \gamma|u_T|) - \frac{\sigma_c}{\delta_{c,N}}(u_N + \gamma|u_T|)^2 + \mathcal{I}_{R+}(u_N) + \mathcal{I}_{[0,1]} \left( 1 - \frac{u_N + \gamma|u_T|}{\delta_{c,N}} \right). \quad (49) $$

The linear evolution of $\beta$, obtained as a consequence of (48), is depicted in Figure 1.
<figure><img src="image_1.png" /><figcaption>Figure 1: (a) The shape of the surface potential $\Psi_s$ with respect to $\beta$. (b) The consequent linear evolution of $\beta$ with $u_T$, assuming $u_N = 0$. (c) The shape of the surface potential $\Psi_s$ with respect to $u_T$, assuming $u_N = 0$, where the initial slope is equivalent to the value of $\sigma_c\gamma$, and the final value is $w$.</figcaption></figure>

The model is now specified for a triangle law by inserting (48) in (32):

$$ \left\{ \begin{array}{l} \dot{\beta} = -\lambda, \quad \sigma_c \delta_{c,N}(\beta - 1) + A^r + \sigma_c u_N + \sigma_c \gamma |u_T| = \xi, \quad r_N^c = r_N + \beta \sigma_c, \quad 0 \le r_N^c \perp u_N \ge 0, \\ 0 \le \xi \perp \beta \ge 0, \quad 0 \le \lambda \perp A^r \ge 0, \quad -u_T \in N[-\beta \sigma_c \gamma, \beta \sigma_c \gamma](r_T^r), \quad -v_T \in N[-\mu r_N^c, \mu r_N^c](r_T^{ir}). \end{array} \right. \qquad (50) $$

**Analytical expressions for an experiment with a given driven tangential displacement** The model in (50) is sufficiently explicit to study analytically. Let us assume that $\beta(0) = \beta_0 \le 1$, and that $v_T$ is a function of time $t$ given by the following piecewise linear function:

$$ v_T(t) = \begin{cases} \frac{1}{2} & \text{for } 0 \le t < 1, \\ -\frac{1}{2} & \text{for } 1 \le t < 2, \\ \frac{1}{2} & \text{for } 2 \le t, \end{cases} \qquad (51) $$

which produces the tangential displacements

$$ u_T(t) = \begin{cases} \frac{1}{2}t & \text{for } 0 \le t < 1, \\ 1 - \frac{1}{2}t & \text{for } 1 \le t < 2, \\ -1 + \frac{1}{2}t & \text{for } 2 \le t, \end{cases} \qquad (52) $$

assuming $u_T(0) = 0$, and where the displacement $u_T$ is measured in mm, and the time $t$ in ms. The time integration of the model described in (50) leads to the following piecewise linear response, where we use $\sigma_c = 0.5$ MPa, $\delta_{c,N} = 1$ mm, $\gamma = 1$ and a friction coefficient $\mu = 0.5$. We assume that the evolution is continuous, and that $u_N = 0$ i.e. that the system is always in contact.

* **First loading phase** $0 \le t < 1$

Since $v_T(t) > 0$ for $t \in [0, 1)$, the dissipative frictional force is $r_T^{ir}(t) = -\mu r_N^c(t)$. Similarly, since $u_T(t) \ge 0$, the tangential cohesive force is $r_T^r(t) = -\beta(t)\sigma_c\gamma$. Let us assume that $\beta(t) > 0$ for $t \in [0, 1)$. We deduce that $\xi(t) = 0$. Let us note that $A^r(0) = -\sigma_c\delta_{c,N}(\beta_0 - 1) - \sigma_c\gamma|u_T(0)| = 0$. Let us assume that $\lambda(t) = 0$, $t \in [0, \varepsilon]$, $\varepsilon > 0$ or equivalently $\dot{\beta}(t) = 0$, $t \in [0, \varepsilon]$, $\varepsilon > 0$. In that case, we get $\dot{A}^r(t) = -\sigma_c\gamma v_T(t) < 0$ and then $A^r(\varepsilon) < 0$, for $\varepsilon > 0$ which is impossible as $A^r$ is constrained to be non-negative. Let us try with $\dot{\beta}(t) < 0$, then $A^r(t) = 0$ and $\beta(t) = 1 - \frac{\gamma|u_T(t)|}{\delta_{c,N}}$ and $\dot{\beta}(t) = -\frac{\gamma v_T(t)}{\delta_{c,N}} < 0$. Since $\beta(1) = 1 - \frac{\gamma}{2\delta_{c,N}} > 0$, this is the only consistent solution for $t \in [0, 1)$.

* **Unloading phase** $1 \le t < 2$

Since $v_T(t) < 0$ for $t \in [1, 2)$, the dissipative frictional force is $r_T^{ir}(t) = \mu r_N^c(t)$. Similarly, since $u_T(t) \ge 0$, the tangential cohesive force is $r_T^r(t) = -\beta(t)\sigma_c\gamma$. Let us assume that $\beta(t) > 0$ for $t \in [1, 2)$ and $\xi(t) = 0$. Let us assume that $\lambda(t) = 0$, $t \in [1, 2)$ or equivalently $\dot{\beta}(t) = 0$. In that case, we get $\dot{A}^r(t) = -\sigma_c\gamma v_T(t) > 0$. Hence $A^r(t) = -\sigma_c\delta_{c,N}(\beta(1) - 1) - \sigma_c\gamma|u_T(t)| > 0$. This solution satisfies the complementarity condition up to $t = 2$.

* **Second loading phase** $2 \le t$

Since $v_T(t) > 0$ for $t \ge 2$, the dissipative frictional force is $r_T^{ir}(t) = -\mu r_N^c(t)$. Similarly, since $u_T(t) \ge 0$, the tangential cohesive force is $r_T^r(t) = -\beta(t)\sigma_c\gamma$. Let us assume that $\beta(t) > 0$ for $t \in [t_1, t_2)$ and $\xi(t) = 0$. Let us assume that $\lambda(t) = 0$, $t \in [2, 2 + \varepsilon]$, $\varepsilon > 0$ or equivalently $\dot{\beta}(t) = 0$. In that case, we get $\dot{A}^r(t) = -\sigma_c\gamma v_T(t) < 0$ and then $A^r(t) = -\sigma_c\delta_{c,N}(\beta(1) - 1) - \sigma_c\gamma|u_T(t)|$ which is positive for $t < 3$. For $t \ge 3$, the only possible solution of the complementarity leads to $\dot{\beta}(t) = -\frac{\gamma v_T(t)}{\delta_{c,N}} < 0$ and $A^r(t) = 0$. The cohesion variable $\beta$ is then $\beta(t) = 1 - \frac{\gamma|u_T(t)|}{\delta_{c,N}}$ which is positive for $t < 4$. For $t \ge 4$, the solution is $\beta(t) = 0$, $A^r(t) = 0$ and $\xi(t) = \sigma_c\gamma u_T(t) - \sigma_c\delta_{c,N}$.
As we are interested in the tangential response, we will simulate the system as if some normal pressure is being applied by setting $r_N(t) = 1 + \beta(t)\sigma_c$ MPa, in order to consistently obtain a frictional force. The solution of this experiment is depicted in Figure 2.

<figure><img src="image_2.png" /><figcaption>Figure 2: Illustration of the extrinsic cohesize law with a linear evolution of cohesion. (a) The displacement $u_T$ as a function of time $t$. (b) The cohesion $\beta$ as a function of time $t$. (c) The slack variable $\lambda$ as a function of time $t$. (d) The slack variable $\xi$ as a function of time $t$. (e) The thermodynamic driving force $A^r$ as a function of time. (f) The cohesion $\beta$ as a function of displacement $u_T$. (g) The slack variable $r_N^c$ as a function of time $t$. (h) The reversible tangential force $r_T^r$ and the irreversible tangential force $r_T^{ir}$ as a function of time $t$.</figcaption></figure>

**Remark 5.** When the interface has fully decohered, we retrieve a standard unilateral contact model with Coulomb friction. This model separates clearly cohesion and friction, while treating them within the same framework. While the presence of cohesion can affect the normal confining force, and exerts an additional restoring force to move the system back towards the $u_T = 0$ state, the frictional force in and of itself depends only on the direction of sliding. This is in contrast to other models wherein the frictional resistance depends not only the sliding direction, but also on the value of $u_T$ (e.g. see Raous et al. (1999) or Parrinello and Borino (2020) as two examples among many).
# 3 Nonsmooth elasto-dynamics of finite-dimensional systems

Now, we will extend our model to consider bodies with finite numbers of degrees of freedom, which may possess masses and stiffnesses, and to which external forces may be applied. These degrees of freedom may form part of a cohesive zone, but equally may not. We also consider the dynamic interaction of multiple bodies via the formulation of impact laws for the system.

## 3.1 Finite-dimensional systems via space-discretisation

We start by considering a finite-dimensional model of a linear elastic mechanical system, of the sort which may arise after a space-discretisation using the finite element method (FEM). Let us note by $u \in \mathbb{R}^n$ the displacements of the system and $v = \dot{u}$ the velocity. Starting from the principle of virtual power (7), the equilibrium equation can be written as

$$M\dot{v} + Ku = F, \quad \dot{u} = v, \qquad (53)$$

where $M \in \mathbb{R}^{n \times n}$ is the mass matrix, assumed to be symmetric positive definite, $K \in \mathbb{R}^{n \times n}$ is the symmetric positive semi-definite structural stiffness matrix and $F \in \mathbb{R}^n$ is the external applied force.

Let us now add that the cohesive zone is applied on a finite set of cohesive contact points that are indexed by $\alpha \in [1, m]$. We will work in a system of small perturbations, and hence the local normal displacements at contact $u_N = \text{col}(u_N^\alpha, \alpha \in [1, m])$ are related to the displacements $u$ by a linear relation written as

$$u_N = H_N u + b_N, \qquad (54)$$

where $H_N \in \mathbb{R}^{m \times n}$ is a selection and interpolation matrix for the degrees of freedom and $b_N \in \mathbb{R}^m$ is a linear adjustment to express misalignment between nodes. In a similar fashion, the local tangential displacements at contact $u_T = \text{col}(u_T^\alpha, \alpha \in [1, m])$ are related to the global displacements $u$ by

$$u_T = H_T u + b_T, \qquad (55)$$

where $H_T \in \mathbb{R}^{m \times n}$ is a selection and interpolation matrix for the degrees of freedom and $b_T \in \mathbb{R}^m$ is another linear adjustment expressing misalignment between nodes. Collecting all variables at contact in the same way ($x = \text{col}(x^\alpha, \alpha \in [1, m])$), the equilibrium equations of the system are given by

$$\begin{cases} M\dot{v} + Ku = F + H_N^\top S r_N + H_T^\top S r_T, \\ \dot{u} = v, \quad u_N = H_N u + b_N, \quad H_T u + b_T, \end{cases} \qquad (56)$$

coupled with (50) to describe the cohesive zone behaviour. The matrix $S \in \mathbb{R}^{m \times m}$ is a diagonal matrix that contains the tributary area of each cohesive zone node after space-discretisation of the interface. For the sake of simplicity, we assume that the parameters of the cohesive zone model do not depend on $\alpha$, but this can be straightforwardly extended.

## 3.2 Nonsmooth dynamics and impacts

In the presence of unilateral contacts, solutions of finite-dimensional dynamical systems with a mass matrix with finite masses associated with degrees of freedom may exhibit jumps in velocities. In this context, the nonsmooth dynamics must be carefully treated to obtain a consistent time-discretisation (Moreau, 1999). The velocity is considered as a function of bounded variations and we denote by $v^-$ and $v^+$ the left and right limits of $v$ with respect to time. We introduce the equations of motion of a discrete (finite-dimensional) mechanical system with reactions at contacts written in terms of differential measures by

$$\begin{cases} M dv + Ku dt = F dt + H_N^\top di_N + H_T^\top di_T, \\ \dot{u} = v, \end{cases} \qquad (57)$$

where $dv$ is the differential measure associated with the velocity $v$, $di_N$ is the measure of the normal reaction at the contact, and $di_T$ is the measure of the tangential reaction at the contact. Almost everywhere (with respect to $dt$), the system is non-impulsive and is equivalent to the dynamics in (56). It follows that the densities of $di_T$ and $di_T$ with respect to $dt$ are given by

$$\frac{di_N}{dt} = S r_N = S(r_N^c - \beta \sigma_c), \quad \text{and} \quad \frac{di_T}{dt} = S r_T = S(r_T^r + r_T^{ir}), \quad \text{dt-almost everywhere.} \qquad (58)$$

To simplify the system in view of numerical time integration, we split the contribution of the cohesion and the contribution of unilateral contact and friction. To this aim, we introduce two measures $di_N^c$ and $di_T^c$ such that the density with respect to $dt$ satisfies

$$\frac{di_N^c}{dt} = S r_N^c, \quad \text{and} \quad \frac{di_T^c}{dt} = S r_T^r, \quad \text{dt-almost everywhere.} \qquad (59)$$
This yields

$$ \frac{d\dot{i}_N}{dt} = -S\beta\sigma_c + \frac{d\dot{i}_N^c}{dt}, \quad \text{and} \quad \frac{d\dot{i}_T}{dt} = Sr_T^r + \frac{d\dot{i}_T^c}{dt}, \quad \text{dt-almost everywhere.} \qquad (60) $$

For the impulsive part of the motion, we assume that the densities of $di_N$ (respectively $di_T$) and $di_N^c$ (respectively $di_T^c$) with respect to a discrete measure $dr_N^c$ (a sum of Dirac atoms $\sum_i \delta_{t_i}$) are equal and denoted by $p_N$ (respectively $p_T$), that is

$$ p_N = \frac{di_N^c}{dr_N^c} = \frac{di_N}{dr_N^c}, \quad \text{and} \quad p_T = \frac{di_T^c}{dr_T^c} = \frac{di_T}{dr_T^c}, \quad \text{dr_N^c-almost everywhere.} \qquad (61) $$

The relation (60) and (61) are contained in

$$ di_N = di_N^c - S\beta\sigma_c dt \quad \text{and} \quad di_T = di_T^c + Sr_T^r dt. \qquad (62) $$

Doing so, we implicitly assume that the cohesion forces ($-S\beta\sigma_c$ and $Sr_T^r$) are Lebesgue integrable functions and the corresponding measure $di - di^c$ does not contain Dirac atoms and is just represented by $-S\beta\sigma dt$ and $Sr_T^r dt$. Since the evolution of the cohesion variable (and hence the driving force) is governed by the displacements $u_N$ and $u_T$, that are both assumed to be absolutely continuous in time, we can expect to have the assumed regularity for $\beta$ and the cohesion forces.

**Additional impact laws.** Since the contact and friction measures $dp_N$ and $dp_T$ may contain atoms, some further modelling choices must be added to properly define $p_N$ and $p_T$. For the normal part of a contact $\alpha$, the Moreau impact law is chosen as

$$ 0 \le p_N^\alpha \perp v_N^{\alpha,+} + ev_N^{\alpha,-} \ge 0 \text{ if } u_N^\alpha \le 0, \text{ else } p_N^\alpha = 0, \qquad (63) $$

where $e \in [0, 1]$ is the Newton coefficient of restitution. If an impact occurs, $p_N^\alpha > 0$ implies the Newton impact law $v_N^{\alpha,+} = -ev_N^{\alpha,-}$. If not, the velocity is continuous and we have the constraints $(1+e)v_N^\alpha \ge 0$. Almost everywhere, we have for the normal part

$$ 0 \le r_N^\alpha + \beta^\alpha \sigma_c \perp u_N^\alpha \ge 0, \text{ that is } 0 \le r_N^{c,\alpha} \perp u_N^\alpha \ge 0. \qquad (64) $$

The Signorini condition (64) is implied by the following condition at the velocity level:

$$ 0 \le r_N^{c,\alpha} \perp v_N^\alpha \ge 0 \text{ if } u_N^\alpha \le 0, \text{ else } r_N^{c,\alpha} = 0, \qquad (65) $$

that we will keep for the formulation. The relations (63) and (65) are contained in the following measure complementarity condition:

$$ 0 \le di_N^{c,\alpha} \perp v_N^{\alpha,+} + ev_N^{\alpha,-} \ge 0 \text{ if } u_N^\alpha \le 0, \text{ else } di_N^{c,\alpha} = 0. \qquad (66) $$

For the tangential part $p_T^\alpha$, we choose the Coulomb friction impact model of Frémond (Frémond, 2017) to ensure the energy dissipation at impact when coupled with a Newton impact law (see Acary and Collins-Craft (2025) for details):

$$ -p_T^\alpha \in \mu p_N^\alpha \operatorname{sgn}(\frac{1}{2}(v_T^{\alpha,+} + v_T^{\alpha,-})). \qquad (67) $$

Doing so, the power dissipated by the frictional impulse is always negative:

$$ \frac{1}{2}(v_T^{\alpha,+} + v_T^{\alpha,-})p_T^\alpha \le 0. \qquad (68) $$

Almost everywhere, the Coulomb friction model is given by

$$ -r_T^{\alpha,\text{ir}} \in \mu(r_N^\alpha + \beta^\alpha \sigma_c) \operatorname{sgn}(v_T^\alpha) = \mu r_N^{c,\alpha} \operatorname{sgn}(v_T^\alpha). \qquad (69) $$

The relations (67) and (69) are contained in the following measure relation:

$$ -di_T^{c,\alpha} \in \mu di_N^{c,\alpha} \operatorname{sgn}(\frac{1}{2}(v_T^{\alpha,+} + v_T^{\alpha,-})). \qquad (70) $$

**Remark 6.** Rigorously, almost everywhere with respect to dt, the relations (66) and (70) are

$$ \left\{ \begin{array}{l} 0 \le Sr_N^{c,\alpha} \perp v_N^\alpha \ge 0 \text{ if } u_N^\alpha \le 0, \text{ else } Sr_N^{c,\alpha} = 0, \\ -Sr_T^{\alpha,\text{ir}} \in \mu Sr_N^{c,\alpha} \operatorname{sgn}(v_T^\alpha). \end{array} \right. \qquad (71) $$

As for the unilateral contact and Coulomb friction case, this relation involve cones and can be rescaled by $S$ ending up with (65) and (69).
To keep the notation compact, the index set $I = \{\alpha, u_N^\alpha \le 0\}$ is defined and the following notation $di_N^c = \text{col}(di_N^{c,\alpha}, \alpha \in I)$, $di_N^c = \text{col}(di_N^{c,\alpha}, \alpha \in I)$, $v_N = \text{col}(v_N^\alpha, \alpha \in I)$, $\bar{H}_N = \text{row}(H_{N,\alpha\bullet}, \alpha \in I)$, $\bar{H}_T = \text{row}(H_{T,\alpha\bullet}, \alpha \in I)$ is adopted. Thus, we can write the equations of the system in a nonsmooth setting using differential measures as

$$ 
\left\{ 
\begin{array}{l}
M dv + Ku dt = F dt + \bar{H}_N^\top dp_N + \bar{H}_T^\top dp_T - H_N^\top S\beta\sigma_c dt + H_T^\top Sr_T^r dt, \\
\dot{u} = v, \quad u_N = H_N u + b_N, \quad u_T = H_T u + b_T, \quad v_N = \bar{H}_N v, \quad v_T = \bar{H}_T v, \\
\dot{\beta} = -\lambda, \quad A^r + \sigma_c u_N + \sigma_c \gamma |u_T| + \sigma_c \delta_{c,N}(\beta - 1) = \xi, \\
0 \le \xi \perp \beta \ge 0, \quad 0 \le \lambda \perp A^r \ge 0, \quad -r_T^r = \beta \sigma_c \gamma \text{sgn}(u_T), \\
0 \le di_N^c \perp v_N^+ + ev_N^- \ge 0, \quad -di_T^c \in \mu di_N^c \text{sgn}(\frac{1}{2}(v_T^+ + v_T^-)).
\end{array}
\right.
\qquad (72) $$

**Remark 7.** While it would appear that the index sets outlined above would exclude contact points on surfaces not formed by fracture from consideration in our model, in fact such points can easily be included in the matrices $H_N$ and $H_T$ by considering them as cohesive zones where the value of $\beta$ has been set to zero. We then recover the standard model of unilateral contact and friction.

# 4 Numerical time integration

Here we present the time-discretisation of the system that allows us to write a fully-discretised system suitable for numerical integration. In particular, we demonstrate that we are able to write the discrete system as a linear complementarity problem, and that the solution of this problem exists. We then finally demonstrate that the discrete energy balance is in general dissipative and energy-preserving in the absence of contact and friction.

## 4.1 Principles of the time integration scheme

We choose a time-integration scheme that is based on the Moreau-Jean scheme adapted for the Frémond impact law for friction (Acary and Collins-Craft, 2025), which is in turn based on the classical Moreau-Jean scheme (Acary and Brogliato, 2008; Jean, 1999; Jean and Moreau, 1992; Moreau, 1999), that is widely used in contact dynamics. We consider a time interval $t \in [t_0, t_{\text{Max}}]$, which is discretised by the sequence $\{t_0, t_1, \dots, t_k, t_{k+1}, \dots, t_{\text{Max}}\}$. For the impulsive terms which appear in (72), the measure of the time interval $(k, k+1]$ is kept as a primary unknown:

$$ p_{N,k,k+1} \approx di_N^c((k, k+1)) = \int_{(k,k+1]} di_N^c \quad \text{and} \quad p_{T,k,k+1} \approx di_T^c((k, k+1)) = \int_{(k,k+1]} di_T^c. \qquad (73) $$

We approximate all of the Lebesque integrable terms using a $\theta$-method given by

$$ \int_{t_k}^{t_{k+1}} x(t) dt \approx hx_{k+\theta}, \qquad (74) $$

where we use the notation $x_{k+\theta} = \theta x_k + (1-\theta)x_{k+1}$ with $\theta \in (0, 1]$ and $h$ is the size of the time-step. For the cohesive reaction force, we have

$$ \int_{(k,k+1]} di_N = \int_{(k,k+1]} di_N^c - S\sigma_c \int_{t_k}^{t_{k+1}} \beta dt \approx p_{N,k,k+1} - hS\sigma_c\beta_{k+\theta}, \qquad (75) $$

and

$$ \int_{(k,k+1]} di_T = \int_{(k,k+1]} di_T^c + S \int_{t_k}^{t_{k+1}} r_T^r dt \approx p_{T,k,k+1} + hSr_{T,k,k+1}^r, \qquad (76) $$

This yields for the first line of (72)

$$ M(v_{k+1} - v_k) + hK u_{k+\theta} = hF_{k+\theta} - h\sigma_c H_N^\top S\beta_{k+\theta} + hH_T^\top r_{T,k+\theta}^r + \bar{H}_N^\top p_{N,k,k+1} + \bar{H}_T^\top p_{T,k,k+1}. \qquad (77) $$

The velocities and displacements are approximated by $v_{k+1} \approx v(t_{k+1})$ and $u_{k+1} \approx u(t_{k+1})$, and by application of the $\theta$-method are linked by

$$ u_{k+1} = u_k + h v_{k+\theta}. \qquad (78) $$

For the contact and friction part (the fifth line of (72)), we choose our contact detection distance $\tilde{u}_{N,k}$ to simply be the relative displacement at the previous time step i.e.

$$ \tilde{u}_{N,k} = u_{N,k}. \qquad (79) $$

Evidently, other discrete approximations of this distance are possible and are feasible implementations of the system, we simply must choose one in order to have a condition that activates the contact and allows the normal and frictional percussion to
be nonzero. In the following, we redefine the index set in discrete time $I_k = \{\alpha, \tilde{u}_{N,k}^{\alpha} \le 0, v_{N,k} < 0\}$ and the following compact notation $p_{N,k,k+1} = \text{col}(p_{N,k,k+1}^{\alpha}, \alpha \in I_k)$, $p_{\tau,k,k+1} = \text{col}(p_{\tau,k,k+1}^{\alpha}, \alpha \in I_k)$, $v_{N,k+1} = \text{col}(v_{N,k+1}^{\alpha}, \alpha \in I_k)$, $\tilde{H}_N = \text{row}(H_{N,\alpha\bullet}, \alpha \in I_k)$, $\tilde{H}_{\tau} = \text{row}(H_{\tau,\alpha\bullet}, \alpha \in I_k)$. The local relative normal and tangential velocities and displacements are obtained from the global velocities and displacements by applying the relevant selection matrices:

$$v_{N,k+1} = \tilde{H}_N v_{k+1}, \quad v_{\tau,k+1} = \tilde{H}_{\tau} v_{k+1}, \qquad (80)$$

and

$$u_{N,k+1} = H_N u_{k+1} + b_{N,k+1} \quad \text{and} \quad u_{\tau,k+1} = H_{\tau} u_{k+1} + b_{\tau,k+1}. \qquad (81)$$

The $\beta$ variable is integrated following the same principle

$$\beta_{k+1} = \beta_k - h\lambda_{k+\theta}, \qquad (82)$$

and the cohesive part of the model (third and fourth lines of (72)) are expressed with the $k + \theta$ approximations:

$$\left\{ \begin{array}{l} A^r_{k+\theta} + \sigma_c u_{N,k+\theta} + \sigma_c \gamma |u_{\tau,k+\theta}| + \sigma_c \delta_{c,N} (\beta_{k+\theta} - 1) = \xi_{k+\theta}, \\ 0 \le \xi_{k+\theta} \perp \beta_{k+\theta} \ge 0, \quad 0 \le \lambda_{k+\theta} \perp A^r_{k+\theta} \ge 0, \\ -r^r_{\tau,k+\theta} = \sigma_c \gamma \beta_{k+\theta} \operatorname{sgn}(u_{\tau,k+\theta}). \end{array} \right. \qquad (83)$$

The discretisation of the contact and friction are given by

$$\left\{ \begin{array}{l} 0 \le p_{N,k,k+1} \perp v_{N,k+1} + ev_{N,k} \ge 0, \\ -p_{\tau,k,k+1} \in \mu p_{N,k,k+1} \operatorname{sgn}(v_{k+\theta}), \end{array} \right. \qquad (84)$$

which can be expressed in terms of $v_{N,k+\theta}$ as

$$\left\{ \begin{array}{l} 0 \le p_{N,k,k+1} \perp \theta v_{N,k+\theta} + (\theta(1+e) - 1)v_{N,k} \ge 0, \\ -p_{\tau,k,k+1} \in \mu p_{N,k,k+1} \operatorname{sgn}(v_{k+\theta}). \end{array} \right. \qquad (85)$$

Following the principles above, the time-stepping scheme for the full elasto-dynamic cohesive-frictional-contact problem is written as follows:

$$\left\{ \begin{array}{l} M(v_{k+1} - v_k) + hK u_{k+\theta} = hF_{k+\theta} - h\sigma_c H_N^\top S \beta_{k+\theta} + hH_T^\top r_{\tau,k+\theta}^r + \tilde{H}_N^\top p_{N,k,k+1} + \tilde{H}_T^\top p_{\tau,k,k+1}, \\ u_{k+1} = u_k + h v_{k+\theta}, \\ u_{N,k+\theta} = H_N u_{k+\theta} + b_{N,k+\theta}, \quad u_{\tau,k+\theta} = H_{\tau} u_{k+\theta} + b_{\tau,k+\theta}, \\ v_{N,k+\theta} = \tilde{H}_N v_{k+\theta}, \quad v_{\tau,k+\theta} = \tilde{H}_{\tau} v_{k+\theta}, \\ \beta_{k+1} = \beta_k - h\lambda_{k+\theta}, \\ \sigma_c \delta_{c,N} (\beta_{k+\theta} - \mathbf{1}) + \sigma_c u_{N,k+\theta} + \sigma_c \gamma |u_{\tau,k+\theta}| + A^r_{k+1} = \xi_{k+\theta}, \\ 0 \le A^r_{k+\theta} \perp \lambda_{k+\theta} \ge 0, \\ 0 \le \beta_{k+\theta} \perp \xi_{k+\theta} \ge 0, \\ -r^r_{\tau,k+\theta} = \beta_{k+\theta} \sigma_c \gamma \operatorname{sgn}(u_{\tau,k+\theta}), \\ 0 \le p_{N,k,k+1} \perp \theta v_{N,k+\theta} + (\theta(1+e) - 1)v_{N,k} \ge 0, \\ -p_{\tau,k,k+1} \in \mu p_{N,k,k+1} \operatorname{sgn}(v_{k+\theta}), \end{array} \right. \qquad (86)$$

where $\mathbf{1}$ represents a vector of ones.

## 4.2 The discrete linear complementarity problem

We choose to formulate our model in the form of a Linear Complementarity Problem (LCP) (Cottle et al., 2009) defined for a matrix $L$ and a vector $q$ as

$$\left\{ \begin{array}{l} w = Lz + q, \\ 0 \le w \perp z \ge 0, \end{array} \right. \qquad (87)$$

and denoted by LCP(L, q). To obtain an LCP, the terms involving the sign functions need to be expressed as a complementarity condition. To this end, we will use the cone reformulations developed in Lemma 2 and Lemma 3. Lemma 2 will be used for friction. For tangential cohesion, we will use the reformulation of Lemma 3, which avoids the addition of spurious solutions.
Lemma 3 requires the introduction of a regularisation variable $\eta_q$ that very slightly increases the cohesion threshold. The system

$$ egin{cases} 0 \le \beta_{k+\theta} \perp \xi_{k+\theta} \ge 0, \\ -r_{\tau,k+\theta}^r = (\sigma_c \gamma \beta_{k+\theta} + \eta_q) \operatorname{sgn}(u_{\tau,k+\theta}), \\ 0 \le p_{N,k,k+1} \perp \theta v_{N,k+\theta} + (\theta(1+e)-1)v_{N,k} \ge 0, \\ -p_{\tau,k,k+1} \in \mu p_{N,k,k+1} \operatorname{sgn}(v_{k+\theta}), \end{cases} \quad (88) $$

is replaced by

$$ egin{cases} r_{\tau,k+\theta}^r = D \hat{r}_{\tau,k+\theta}^r, \\ 0 \le S\beta_{k+\theta} \perp \xi_{k+\theta} \ge 0, \\ 0 \le \hat{r}_{\tau,k+\theta}^r \perp 1\chi_{k+1} + D^\top u_{\tau,k+\theta} \ge 0, \\ 0 \le \chi_{k+\theta} \perp (\sigma_c \gamma S\beta_{k+\theta} + \eta_q) - 1^\top \hat{r}_{\tau,k+\theta}^r \ge 0, \\ p_{\tau,k+\theta}^r = D \hat{p}_{\tau,k+\theta}, \\ 0 \le p_{N,k,k+1} \perp v_{N,k+\theta} + (\theta(1+e)-1)v_{N,k} \ge 0, \\ 0 \le \hat{p}_{\tau,k,k+1} \perp 1\zeta_{k+\theta} + D^\top v_{\tau,k+\theta} \ge 0, \\ 0 \le \zeta_{k+\theta} \perp \mu p_{N,k,k+1} - 1^\top \hat{p}_{\tau,k,k+1} \ge 0. \end{cases} \quad (89) $$

We introduced two new slack variables $\zeta_{k+\theta} \in \mathbb{R}_m^+$ and $\chi_{k+\theta} \in \mathbb{R}_m^+$ and the matrices $D \in \mathbb{R}^{m \times 2m}$ and $1 \in \mathbb{R}^{2m \times m}$ defined by

$$ D = \begin{bmatrix} 1 & -1 & 0 & 0 & \dots & 0 & 0 \\ 0 & 0 & 1 & -1 & \dots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & \dots & 1 & -1 \end{bmatrix}, \quad \text{and} \quad 1 = \begin{bmatrix} 1 & 0 & \dots & 0 \\ 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 \\ 0 & 0 & \dots & 1 \end{bmatrix}. \quad (90) $$

We must also consider our expression for the slack variable $\xi_{k+\theta} = A_{k+\theta}^r + \sigma_c \delta_{c,N} (\beta_{k+\theta} - 1) + \sigma_c u_{N,k+\theta} + \sigma_c \gamma |u_{\tau,k+\theta}|$. Our expression for $\xi_{k+\theta}$ features a dependence on an absolute value $|u_{\tau,k+\theta}|$. This does not pose any problems for the trailing terms, but must be eliminated for the terms dependant on the entries of $z$ in order to obtain an LCP. As shown in Lemma 3, we have $|u_{\tau,k+\theta}| = \chi_{k+\theta}$.

Now, in order to obtain a workable formulation to express as an LCP, we first consider the entries in the complementarity variable vectors $w$ and $z$ that are given by:

$$ w = \begin{bmatrix} h\theta\lambda_{k+\theta} \\ \xi_{k+\theta} \\ \theta v_{N,k+\theta} + \theta(\theta(1+e)-1)v_{N,k} \\ 1\theta\zeta_{k+\theta} + \theta D^\top v_{\tau,k+\theta} \\ \mu p_{N,k,k+1} - 1^\top \hat{p}_{\tau,k,k+1} \\ 1\chi_{k+\theta} + D^\top u_{\tau,k+\theta} \\ \sigma_c \gamma S\beta_{k+\theta} - 1^\top \hat{r}_{\tau,k+\theta}^r \end{bmatrix}, \quad z = \begin{bmatrix} SA_{k+\theta}^r \\ S\beta_{k+\theta} \\ p_{N,k,k+1} \\ \hat{p}_{\tau,k,k+1} \\ \theta\zeta_{k+\theta} \\ \hat{r}_{\tau,k+\theta}^r \\ \chi_{k+\theta} \end{bmatrix}, \quad (91) $$

where all of the cohesive zone variables should be understood as vectors. As in Collins-Craft et al. (2022), we take $h\theta\lambda_{k+\theta}$ as the complementarity variable (rather than $\theta\lambda_{k+\theta}$) in order to avoid an ill-conditioned matrix as $h \to 0$. Some of the variables have been rescaled by $\theta > 0$ to facilitate the proof of existence of solutions (see §4.3). We detail the algebraic manipulations required to arrive at the LCP formulation in Appendix B.1, and present here only the final result. Once we arrive at the formulation of the LCP, we add a small numerical regularisation term $\eta_L$ to the main diagonal of the $L$ matrix to facilitate the operation of the pivot solver when the pivot values are very close to zero. In the most general case of multiple cohesive zones with boundary
conditions enforced, the LCP($L, q$) is defined by

$$ L = \begin{bmatrix} \mathbf{0}^{m \times m} & -S^{-1} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ S^{-1} & \sigma_c(\delta_{c,N}S^{-1} - h^2\theta^2\sigma_c U_{NN}) & h\theta^2\sigma_c V_{NN}^{\top} & h\theta^2\sigma_c V_{TN}^{\top} D & \mathbf{0}^{m \times m} & h^2\theta^2\sigma_c U_{NT} D & \sigma_c\gamma I^{m \times m} \\ \mathbf{0}^{m \times m} & -h\theta^2\sigma_c V_{NN} & \theta^2 W_{NN} & \theta^2 W_{NT} D & \mathbf{0}^{m \times m} & h\theta^2 V_{NT} D & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & -h\theta^2\sigma_c D^{\top} V_{TN} & \theta^2 D^{\top} W_{TN} & \theta^2 D^{\top} W_{TT} D & \mathbf{1} & h\theta^2 D^{\top} V_{TT} D & \mathbf{0}^{2m \times m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \mu I^{m \times m} & -\mathbf{1}^{\top} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & -h^2\theta^2\sigma_c D^{\top} U_{TN} & h\theta^2 D^{\top} V_{NT}^{\top} & h\theta^2 D^{\top} V_{TT}^{\top} D & \mathbf{0}^{2m \times m} & h^2\theta^2 D^{\top} U_{TT} D & \mathbf{1} \\ \mathbf{0}^{m \times m} & \sigma_c\gamma I^{m \times m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & -\mathbf{1}^{\top} & \mathbf{0}^{m \times m} \end{bmatrix}, $$

$$ q = \begin{bmatrix} \beta_k \\ \sigma_c (q_{u_N} - \delta_{c,N} \mathbf{1}^m) \\ \theta \bar{H}_N \hat{M}^{-1} \hat{i}_{k,k+1} + \theta(\theta(1+e) - 1) \bar{H}_N v_k \\ \theta D^{\top} \bar{H}_T \hat{M}^{-1} \hat{i}_{k,k+1} \\ \mathbf{0}^m \\ D^{\top} q_{u_T} \\ \eta_q \mathbf{1}^m \end{bmatrix}. \quad (92) $$

where $\mathbf{0}^{a \times b}$ corresponds to a matrix of zeros with $a$ rows and $b$ columns, $\eta_q$ is the numerical regularisation factor applied to $q$, and the $U$, $V$ and $W$ matrices are defined in Appendix B.1.

## 4.3 Existence of the solution of the discrete LCP

We want to show that the system LCP($L, q$) has at least one solution. As the system contains Coulomb friction, uniqueness cannot be guaranteed for any value of the coefficient of friction.

**Assumption 1.** The time-step $h$ is chosen small enough that $\sigma_c(\delta_{c,N}S^{-1} - h^2\theta^2\sigma_c U_{NN})$ is positive definite.

Since $\sigma_c\delta_{c,N}S^{-1}$ is a positive definite matrix, a time-step exists such that Assumption 1 holds.

**Lemma 1.** Under Assumption 1, $L$ is copositive on the positive orthant.

**Proof** Let us recall that $L$ is copositive on the positive orthant if $z^\top L z \ge 0$ for all $z \ge 0$. Since $z^\top L z = \frac{1}{2} z^\top (L + L^\top) z$, we first compute the symmetric part of $L$:

$$ \frac{1}{2} (L + L^{\top}) = \begin{bmatrix} \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ \mathbf{0}^{m \times m} & \sigma_c(\delta_{c,N}S^{-1} - h^2\theta^2\sigma_c U_{NN}) & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \sigma_c\gamma I^{m \times m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \theta^2 W_{NN} & \theta^2 W_{NT} D & \frac{1}{2}\mu I^{m \times m} & h\theta^2 V_{NT} D & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & \mathbf{0}^{2m \times m} & \theta^2 D^{\top} W_{TN} & \theta^2 D^{\top} W_{TT} D & \mathbf{0}^{2m \times m} & h\theta^2 D^{\top} V_{TT} D & \mathbf{0}^{2m \times m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \frac{1}{2}\mu I & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & \mathbf{0}^{2m \times m} & h\theta^2 D^{\top} V_{NT}^{\top} & h\theta^2 D^{\top} V_{TT}^{\top} D & \mathbf{0}^{2m \times m} & h^2\theta^2 D^{\top} U_{TT} D & \mathbf{0}^{2m \times m} \\ \mathbf{0}^{m \times m} & \sigma_c\gamma I^{m \times m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \end{bmatrix}. \quad (93) $$

By splitting the matrix, it is possible to obtain the following relation:

$$ \frac{1}{2} z^{\top} (L + L^{\top}) z = z_2^{\top} (\sigma_c (\delta_{c,N} S^{-1} - h^2 \theta^2 \sigma_c U_{NN}) z_2 + \theta^2 \begin{bmatrix} z_3 \\ z_4 \\ z_6 \end{bmatrix}^{\top} \begin{bmatrix} W_{NN} & W_{NT} D & h V_{NT} D \\ D^{\top} W_{TN} & D^{\top} W_{TT} D & h D^{\top} V_{TT} D \\ h D^{\top} V_{NT}^{\top} & h D^{\top} V_{TT}^{\top} D & h^2 D^{\top} U_{TT} D \end{bmatrix} \begin{bmatrix} z_3 \\ z_4 \\ z_6 \end{bmatrix} + 2 \sigma_c \gamma z_7^{\top} z_2 + \mu z_3^{\top} z_5. \quad (94) $$

Under Assumption 1, the first term is nonnegative

$$ z_2^{\top} [\sigma_c(\delta_{c,N}S^{-1} - h^2\theta^2\sigma_c U_{NN})] z_2 \ge 0, \quad (95) $$
and

$$2\sigma_c\gamma z_7^\top z_2 + \mu z_3^\top z_5 \ge 0 \text{ for all } z \ge 0. \quad (96)$$

For the remaining term in (94), let us rewrite the matrix as

$$\begin{bmatrix} W_{NN} & W_{NT}D & hV_{NT}D \\ D^\top W_{TN} & D^\top W_{TT}D & hD^\top V_{TT}D \\ hD^\top V_{NT}^\top & hD^\top V_{TT}^\top D & h^2D^\top U_{TT}D \end{bmatrix} = \begin{bmatrix} I & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times 2m} \\ \mathbf{0}^{m \times m} & D & \mathbf{0}^{m \times 2m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & hD \end{bmatrix}^\top \begin{bmatrix} \bar{H}_N \\ \bar{H}_T \\ H_T \end{bmatrix} \hat{M}^{-1} \begin{bmatrix} I & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times 2m} \\ \mathbf{0}^{m \times m} & D & \mathbf{0}^{m \times 2m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & hD \end{bmatrix}. \quad (97)$$

The matrix on the left-hand side of (97) has the form $A\hat{M}^{-1}A^\top$. Since $\hat{M}^{-1}$ is a positive definite symmetric matrix, this matrix is also a positive semi-definite matrix, and we can conclude that

$$\theta^2 \begin{bmatrix} z_3 \\ z_4 \\ z_6 \end{bmatrix}^\top \begin{bmatrix} W_{NN} & W_{NT}D & h\theta V_{NT}D \\ D^\top W_{TN} & D^\top W_{TT}D & hD^\top V_{TT}D \\ hD^\top V_{NT}^\top & hD^\top V_{TT}^\top D & h^2D^\top U_{TT}D \end{bmatrix} \begin{bmatrix} z_3 \\ z_4 \\ z_6 \end{bmatrix} \ge 0. \quad (98)$$

Each component of the sum in (94) is nonnegative, and so $L$ is copositive on the positive orthant. $\square$

To prove the existence of solutions, we need to add an assumption on the surjectivity of $H_T$.

**Assumption 2.** The matrix $H_T$ is surjective, i.e, $\forall b \in \mathbb{R}^m, \exists a \in \mathbb{R}^n$ such that $b = H_T a$.

If the constraints are not redundant, i.e if the matrix has full row rank then Assumption 2 is satisfied.

**Proposition 1.** If Assumption 1 and 2 hold then the LCP($L, q$) has a solution.

**Proof** Let us first compute the set of solutions of the homogeneous LCP($L, 0$), also called the kernel of the LCP, given by

$$\mathcal{K}(L) = \{z \ge 0, Lz \ge 0, z^\top Lz = 0\}. \quad (99)$$

Using (94) and (97), the relation $z^\top Lz = 0$ with Assumption 1 is equivalent to

$$z_2 = 0, \bar{H}_N^\top z_3 = 0, \bar{H}_T^\top D z_4 = 0, H_T^\top D z_6 = 0, \quad (100)$$

for $\theta > 0$. The relations $z \ge 0, Lz \ge 0$ further implies

$$z_1 \ge 0, z_5 \ge 0, z_7 \ge 0. \quad (101)$$

The set of solutions of the homogeneous LCP($L, 0$) is then given by

$$\mathcal{K}(L) = \{z \mid z_1 \ge 0, z_2 = 0, \bar{H}_N^\top z_3 = 0, \bar{H}_T^\top D z_4 = 0, z_5 \ge 0, H_T^\top D z_6 = 0, z_7 \ge 0\}. \quad (102)$$

For $z \in \mathcal{K}(L)$, we have

$$\begin{aligned} q^\top z &= \beta_k^\top z_1 + [\theta \bar{H}_N \hat{M}^{-1} \hat{i}_{k,k+1} + \theta(\theta(1+e)-1) \bar{H}_N v_k]^\top z_3 + [\theta D^\top \bar{H}_T \hat{M}^{-1} \hat{i}_{k,k+1}]^\top z_4 + [D^\top q_{u_T}]^\top z_6 + \eta_q \mathbf{1}^\top z_7, \\ &= \beta_k^\top z_1 + [\theta \hat{M}^{-1} \hat{i}_{k,k+1} + \theta(\theta(1+e)-1) v_k]^\top \bar{H}_N^\top z_3 + [\theta \hat{M}^{-1} \hat{i}_{k,k+1}]^\top \bar{H}_T^\top D z_4 + [q_{u_T}]^\top D z_6 + \eta_q \mathbf{1}^\top z_7. \end{aligned} \quad (103)$$

From (100) we know that $\bar{H}_N^\top z_3 = 0$ and $\bar{H}_T^\top D z_4 = 0$ which allows us to simplify (103) as

$$q^\top z = \beta_k^\top z_1 + [q_{u_T}]^\top D z_6 + \eta_q \mathbf{1}^\top z_7. \quad (104)$$

Using the expression of $q_{u_T}$, we can further simplify to

$$q^\top z = \beta_k^\top z_1 + [h\theta \hat{M}^{-1} \hat{i}_{k,k+1} + u_k]^\top H_T^\top D z_6 + b_{T,k+\theta}^\top D z_6 = b_{T,k+\theta}^\top D z_6 + \eta_q \mathbf{1}^\top z_7, \quad (105)$$

since $H_T^\top D z_6 = 0$. Using Assumption 2, there exists an $a$ such that $b_{T,k+\theta}^\top = H_T a$, then we get

$$q^\top z = a^\top H_T^\top D z_6 + \beta_k^\top z_1 + \eta_q \mathbf{1}^\top z_7 = \beta_k^\top z_1 + \eta_q \mathbf{1}^\top z_7 \quad (106)$$

As $\eta_q > 0$ and $z_7 \ge 0, \eta_q \mathbf{1}^\top z_7 \ge 0$ and $\beta_k \ge 0$, we can conclude

$$q^\top z \ge 0 \text{ for all } z \in \mathcal{K}(L). \quad (107)$$

Since $L$ is a copositive matrix, the condition (107) implies that the LCP($L, q$) has a solution by virtue of Theorem 3.8.6 in Cottle et al. (2009). $\square$
## 4.4 Discrete energy balance

In this section, we will show that the time-stepping scheme described in (86) satisfies the discrete energy balance, even in the case of impacts. We will start by considering the momentum balance given in (72), and multiplying by $\frac{1}{2} (v^+ + v^-)^\top$:

$$ \frac{1}{2} (v^+ + v^-)^\top M dv + \frac{1}{2} (v^+ + v^-)^\top Ku dt = \frac{1}{2} (v^+ + v^-)^\top F dt + \frac{1}{2} (v^+ + v^-)^\top H_N^\top (di_N^c - S\sigma_c\beta dt) \\ + \frac{1}{2} (v^+ + v^-)^\top H_T^\top (Sr_T^r dt + di_T^c). \quad (108) $$

As $v^+ dt = v^- dt = v dt$, and $M$ and $K$ are symmetric matrices, we are able to rewrite (108) as

$$ d\left(\frac{1}{2}v^\top Mv\right) + d\left(\frac{1}{2}u^\top Ku\right) = v^\top F dt + \left[\frac{1}{2}H_N(v^+ + v^-)\right]^\top di_N^c - v_N^\top S\sigma_c\beta dt \\ + v_T^\top Sr_T^r dt + \left[\frac{1}{2}H_T(v^+ + v^-)\right]^\top di_T^c. \quad (109) $$

As the space-discretised kinetic and strain energies are given by $K = \frac{1}{2}v^\top Mv$ and $U = \frac{1}{2}u^\top Ku$, we can write

$$ dK + dU = v^\top F dt + \frac{1}{2} (v_N^+ + v_N^-)^\top di_N^c - v_N^\top S\sigma_c\beta dt + v_T^\top Sr_T^r dt + \frac{1}{2} (v_T^+ + v_T^-)^\top di_T^c. \quad (110) $$

We now specify the space-discretised version of the continuous fracture energy given in (45):

$$ dG = (v_N^\top \sigma_c S\beta - v_T^\top Sr_T^r) dt, \quad (111) $$

which lets us write the energy balance of the space-discretised system:

$$ dK + dU + dG = v^\top F dt + \frac{1}{2} (v_N^+ + v_N^-)^\top di_N^c + \frac{1}{2} (v_T^+ + v_T^-)^\top di_T^c. \quad (112) $$

Let us define the total energy of the system as $T = K + U + G$, and integrate (112) over a time interval of $(t_1, t_2]$. We thus obtain the incremental energy balance:

$$ \Delta T = T^+(t_2) - T^-(t_1) = W_{\text{ext}} + W_{\text{impact}} + W_{\text{friction}}, \quad (113) $$

where the work done by the external forces, contact/impact and friction on the time interval $(t_1, t_2]$ is given by

$$ W_{\text{ext}} = \int_{t_1}^{t_2} v^\top F dt, \quad (114) $$

$$ W_{\text{impact}} = \int_{(t_1, t_2]} \frac{1}{2} (v_N^+ + v_N^-)^\top di_N^c, \quad (115) $$

$$ W_{\text{friction}} = \int_{(t_1, t_2]} \frac{1}{2} (v_T^+ + v_T^-)^\top di_T^c. \quad (116) $$

The contact law specified in (63) has been shown to dissipate energy i.e. $W_{\text{impact}} \le 0$ (Acary, 2016). The use of the Frémond impact law for friction specified in (67) allows us to conclude that $W_{\text{friction}} \le 0$ (see (Acary and Collins-Craft, 2025) for details). We can in turn conclude that

$$ \Delta T - W_{\text{ext}} = W_{\text{impact}} + W_{\text{friction}} \le 0. \quad (117) $$

The space-discretised incremental energy balance is thus either conservative or dissipative, that is to say that no spurious energy is added to the system.

We now show that the equivalent incremental energy balance is satisfied by the space-and-time-discretised scheme. Following the method in Acary (2016), the variation of the sum of the kinetic and elastic energy over the time-step is given by

$$ \Delta K_{k,k+1} + \Delta U_{k,k+1} = \left(\frac{1}{2} - \theta\right) \left(\|v_{k+1} - v_k\|_M^2 + \|u_{k+1} - u_k\|_K^2\right) \\ + h v_{k+\theta}^\top F_{k+\theta} + v_{N,k+\theta}^\top p_{N,k,k+1} - h v_{N,k+\theta}^\top \sigma_c S \beta_{k+\theta} \\ + h v_{T,k+\theta}^\top S r_{T,k+\theta}^r + v_{T,k+\theta}^\top p_{T,k,k+1}, \quad (118) $$
where $|| \cdot ||_M^2$ and $|| \cdot ||_K^2$ are the 2-norm and 2-seminorm induced by the matrices $M$ and $K$ respectively. Let us define the discrete approximations of the increment of fracture energy, the increment of external work, the increment of contact work and the increment of frictional work:

$$ \Delta G_{k,k+1} = h \left( v_{N,k+\theta}^{\top} \sigma_c S \beta_{k+\theta} - v_{T,k+\theta}^{\top} r_{T,k+\theta}^{\mathrm{r}} \right) \approx \int \left( v_{N}^{\top} \sigma_c S \beta - v_{T}^{\top} S r_{T}^{\mathrm{r}} \right) dt, \quad (119) $$

$$ \Delta W_{\mathrm{ext},k,k+1} = h v_{k+\theta}^{\top} F_{k+\theta} \approx \int_{t_k}^{t_{k+1}} v^{\top} F dt, \quad (120) $$

$$ \Delta W_{\mathrm{impact},k,k+1} = v_{N,k+\theta}^{\top} p_{N,k,k+1}, \quad (121) $$

$$ \Delta W_{\mathrm{friction},k,k+1} = v_{T,k+\theta}^{\top} p_{T,k,k+1}, \quad (122) $$

where we have used the approximations given in (73). Now, provided $\frac{1}{2} \le \theta \le 1$, we have an estimate for the variation of the energy given by

$$ \Delta K_{k,k+1} + \Delta U_{k,k+1} + \Delta G_{k,k+1} - \Delta W_{\mathrm{ext},k,k+1} \le \Delta W_{\mathrm{impact},k,k+1} + \Delta W_{\mathrm{friction},k,k+1}. \quad (123) $$

It has been shown (Acary, 2016) that our chosen discretisation of the impact law leads to

$$ \Delta W_{\mathrm{impact},k,k+1} \le 0 \quad \text{if} \quad \theta \le \frac{1}{1+e}. \quad (124) $$

Similarly, it has also been shown (Acary and Collins-Craft, 2025) that our chosen discretisation of the friction law gives

$$ \Delta W_{\mathrm{friction},k,k+1} \le 0 \quad \text{if} \quad \frac{1}{2} \le \theta \le \frac{1}{1+e} \le 1. \quad (125) $$

We are thus able to obtain the result

$$ \Delta K_{k,k+1} + \Delta U_{k,k+1} + \Delta G_{k,k+1} - \Delta W_{\mathrm{ext},k,k+1} \le \Delta W_{\mathrm{impact},k,k+1} + \Delta W_{\mathrm{friction},k,k+1} \le 0. \quad (126) $$

For $\theta = \frac{1}{2}$, we can sharpen this to

$$ \Delta K_{k,k+1} + \Delta U_{k,k+1} + \Delta G_{k,k+1} - \Delta W_{\mathrm{ext},k,k+1} = \Delta W_{\mathrm{impact},k,k+1} + \Delta W_{\mathrm{friction},k,k+1} \le 0. \quad (127) $$

We can compare this with the continuous-in-space-and-time discrete energy balance in (47), and see that the discretised scheme does not spuriously add energy into the system, and is in fact dissipative. In the special case of no friction or impact, the integration scheme is energy-preserving.

## 4.5 Quasi-static scheme

Often in laboratory tests displacements are applied that allow the rate of crack propagation to be confined to the quasi-static regime, and as such it is useful to have a numerical scheme that is adapted to the case of slow crack growth over long periods of time, with the caveat that such a scheme can only work if the material is stiff and the decohesion process releases a significant amount of energy. Such a caveat is not unique to our model, as it is known that the combination of linear elasticity with a linear cohesive zone model as we have presented here is not always well-posed in quasi-statics (Bouchitté et al., 1995; Rodella et al., 2025). We once again use the $\theta$-method to discretise the system in time, and work with the contact and friction forces, rather than percussions. $D$ and $\mathbb{1}$ retain their meaning as before. In quasi-statics, we solve for normal relative displacement directly, and so we have no need to perform a contact detection step or distinguish between $H$ and $\bar{H}$. Consequently there is no distinction between the Delassus matrices $U, V, W = HK^{-1}H$ (other than those due to the selection of tangential and normal degrees of freedom), so we will denote all of them as $U$. The development of the equations for the quasi-static case is
detailed in Appendix B.2, and here we simply present the LCP($L_{stat}$, $q_{stat}$) that we solve:

$$w_{stat} = \begin{bmatrix} h\theta\lambda_{k+\theta} \\ \xi_{k+\theta} \\ u_{N,k+\theta} \\ 1h\theta\zeta_{k+\theta} + D^{\top}h\theta v_{T,k+\theta} \\ \mu Sr_{N,k+\theta}^{c} - 1^{\top}\hat{r}_{T,k+\theta}^{ir} \\ 1\chi_{k+\theta} + D^{\top}u_{T,k+\theta} \\ \sigma_c\gamma S\beta_{k+\theta} - 1^{\top}\hat{r}_{T,k+\theta}^{r} \end{bmatrix}, \quad z_{stat} = \begin{bmatrix} SA^{r}_{k+\theta} \\ S\beta_{k+\theta} \\ Sr_{N,k+\theta}^{c} \\ \hat{r}_{T,k+\theta}^{ir} \\ h\theta\zeta_{k+\theta} \\ \hat{r}_{T,k+\theta}^{r} \\ \chi_{k+\theta} \end{bmatrix},$$

$$L_{stat} = \begin{bmatrix} \mathbf{0}^{m \times m} & -S^{-1} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ S^{-1} & \sigma_c (\delta_{c,N}S^{-1} - \sigma_c U_{NN}) & \sigma_c U_{NN} & \sigma_c U_{NT} D & \mathbf{0}^{m \times m} & \sigma_c U_{NT} D & \sigma_c \gamma I \\ \mathbf{0}^{m \times m} & -\sigma_c U_{NN} & U_{NN} & U_{NT} D & \mathbf{0}^{m \times m} & U_{NT} D & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & -\sigma_c D^{\top} U_{TN} & D^{\top} U_{TN} & D^{\top} U_{TT} D & 1 & D^{\top} U_{TT} D & \mathbf{0}^{2m \times m} \\ \mathbf{0}^{m \times m} & \mathbf{0}^{m \times m} & \mu I & -1^{\top} & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} \\ \mathbf{0}^{2m \times m} & -\sigma_c D^{\top} U_{TN} & D^{\top} U_{TN} & D^{\top} U_{TT} D & \mathbf{0}^{2m \times m} & D^{\top} U_{TT} D & 1 \\ \mathbf{0}^{m \times m} & \sigma_c \gamma I & \mathbf{0}^{m \times m} & \mathbf{0}^{m \times 2m} & \mathbf{0}^{m \times m} & -1^{\top} & \mathbf{0}^{m \times m} \end{bmatrix},$$

$$q_{stat} = \begin{bmatrix} \beta_k \\ \sigma_c (H_N K^{-1} F_{k+\theta} + b_{N,k+\theta} - \delta_{c,N} \mathbf{1}^m) \\ H_N K^{-1} F_{k+\theta} + b_{N,k+\theta} \\ D^{\top} H_T K^{-1} F_{k+\theta} - D^{\top} H_T u_k \\ \mathbf{0}^m \\ D^{\top} H_T K^{-1} F_{k+\theta} + D^{\top} b_{T,k+\theta} \\ \eta_q \mathbf{1}^m \end{bmatrix}. \qquad (128)$$

We do not develop any proof of the existence of solutions to this LCP, although such a proof would follow straightforwardly the proof of the dynamic scheme in §4.3 as the LCP has the same structure. However, the existence of solutions in this case depends on the material parameters $\sigma_c$, $\delta_{c,N}$, E and $\nu$ taking the appropriate values to ensure that $\sigma_c (\delta_{c,N}S^{-1} - \sigma_c U_{NN})$ is positive definite. As such, materials which are insufficiently stiff cannot be guaranteed to have a solution with this system, and must be simulated with the dynamic LCP regardless of the rate of loading. We also note that even where the stiffness matrix K has been modified appropriately to include Dirichlet boundary conditions, it will only be positive definite so long as all elements remain connected. If any part of the system is able to translate as a rigid body (as may be the case when a crack runs through a system from one side to the other), the matrix will no longer be invertible and the system cannot be solved. Once again, in this case simulation with the dynamic LCP resolves the difficulty.

# 5 Numerical simulations

We implement our discrete system in Python 3.12.1, and resolve the LCP using Siconos 4.5.0 (Acary et al., 2024). The LCP (92) can be solved reliably by Lemke's algorithm (Lemke and Howson, Jr., 1964). Specifically, Lemke's algorithm provides a guarantee that if solutions of the LCP exist (as we have proved above), then the algorithm will find one. The nonsmooth Newton approach would be the most promising for large-scale and parallelised simulations, but for this work where we want to demonstrate the proof-of-concept, we prefer the greater accuracy and mathematical guarantees of Lemke's algorithm.

As in Collins-Craft et al. (2022), in implementation we remove the entries corresponding to the cohesive zone when $\beta$ is less than a certain threshold $c_{cohesion}$, namely $1 \times 10^{-3}$, compared to $10^{-12}$ in Collins-Craft et al. (2022). This larger threshold is chosen for numerical efficiency, as the tangential constraints are two-sided i.e. $\hat{r}_{T,k+1}^r \ge 0$ and $1^{\top}\hat{r}_{T,k+1}^r \le \sigma_c \gamma S \beta_{k+1}$, which becomes difficult for the solver to satisfy at small values of $\beta$ without the use of very small time-steps. While in principle this can be addressed by the use of an enumerative solver, this becomes unworkable for systems of even moderate size, as the scaling properties of enumerative solvers are extremely unfavourable. Hence, we favour the less elegant, but simpler, solution of simply removing these lines from the problem when their resolution becomes overly onerous. We similarly remove those entries corresponding to the contact problem when $u_{N,k} \ge 0$. The effect of both of these removal procedures is to make the LCP smaller, and thus increase the speed of solution.
## 5.1 Quasi-static scalar case with elastic springs

In the interests of demonstrating the benefits of working with the dynamic formulation in (92), we demonstrate a pathology that can arise in quasi-static systems. We consider the case of an elastic spring bound to a rigid substrate at one end, but otherwise free to move, in the spirit of Acary and Monerie (2006) and Chaboche et al. (2001).

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

Figure 3: We model the elastic spring as two nodes connected by a linear elastic material. Node 1 is bound to the rigid substrate and has a normal confining force applied to it while node 2 has a controlled displacement or force applied to it. We add two springs, one in the vertical direction and one running diagonally, and forbid any displacements at the far end.

We are obliged to add a vertical spring to the system, as in the quasi-static case an absence of stiffness in the normal direction causes the system to be insensitive to the applied confining force and changes in the elasticity of the spring. We add a diagonal spring to complete the structure and avoid pathological singularities in the stiffness matrix.

We give the system the material parameters $G_c = 0.25 \text{ N/mm}$, $\sigma_c = 0.5 \text{ MPa}$, $\gamma = 1$, $\mu = 0.5$, a surface area $S = 1 \text{ mm}^2$, a vertical and horizontal spring length of 1 mm, with a corresponding length of $\sqrt{2} \text{ mm}$ for the diagonal spring. We set the regularisation parameter $\eta_q = 1 \times 10^{-8}$ and the time-stepping parameter $\theta = 1$. We apply a driving displacement of $u = 0.5t$ to the far end of the horizontal spring, apply a confining force of $F = -5 \text{ N}$ in the vertical direction at the cohesive zone node, and then vary the stiffness of the bars in the system.

### 5.1.1 The soft spring system

First, setting the bar stiffness to $E = 0.45 \text{ MPa}$, simulating for 18 ms with 360 time steps ($h = 0.05 \text{ ms}$), we have the behaviour shown in Figure 4.
<figure><img src="image_4.png" /><figcaption>Figure 4: The system with $E = 0.45$ MPa. (a) The nodal displacements $u_1$ and $u_3$ as a function of time $t$. (b) The cohesion $\beta$ as a function of time $t$. (c) The cohesion $\beta$ as a function of crack sliding displacement $u_T$, normalised by the critical sliding distance $\delta_{c,T}$. (d) The tangential cohesive for $r_T^r$ as a function of time $t$. (e) The tangential frictional force $r_T^{fr}$ as a function of time $t$. (f) The sum of the cohesive and frictional tangential forces $r_T^r + r_T^{fr}$ as a function of time $t$.</figcaption></figure>

We observe behaviour from our system consistent with it being ill-posed. The tangential displacement of the node to which the cohesive zone is attached jumps rapidly, and the sum of the tangential cohesive and frictional forces likewise experiences a sudden jump. Of particular note, the system decoheres entirely in one time-step, and the graph of the cohesion against the normalised crack sliding displacement demonstrates dramatic “overshoot” (that is to say, if the system behaved correctly, $\beta$ would be zero when $u_T/\delta_{c,T} = 1$, but in this system $\beta$ goes to zero when $u_T/\delta_{c,T} = 1.6694$). As a consequence, an amount of energy substantially more than $G_c$ will be released by the decohesion process, and thus the numerical results are no longer coherent with the analytical model.

### 5.1.2 The stiff spring system

By contrast to the system in §5.1.1, if we instead set the bar stiffness to $E = 4.5$ MPa and simulate the system for 6 ms with 120 time steps, we observe the behaviour shown in Figure 5.
<figure><img src="image_5.png" /><figcaption>Figure 8: The edge-cracked block geometry, with tractions applied to the top and bottom surfaces. The crack will propagate rightwards from the pre-existing crack tip, and at a certain point begin to bifurcate.</figcaption></figure>

While the idealised pre-cut crack geometry is infinitely sharp, in practice specifying this choice causes Gmsh to make the body continuous across the crack. We thus set the lower and upper crack lips to make contact with the left-hand side of the geometry at 19.999 mm and 20.001 mm respectively, which is sufficiently sharp for our purposes. We further treat the block as having a depth of 1 mm, and consequently of being in the plane stress condition. In meshing the block, we instruct Gmsh to produce a mesh with a characteristic length of (at most) 0.5 mm along the centre line from the initial crack tip, as well as at the top right and bottom right corners. The characteristic length at the left-hand side of the block is 5 mm. The combined effect of these instructions is to produce a mesh that is coarser on the left-hand (pre-cracked) side, with a more densely meshed zone that "fans out" from the centre line towards the corners. The mesh is displayed in Figure 9.

<figure><img src="image_6.png" /><figcaption>Figure 9: The mesh generated by Gmsh which we use for the simulation. The visualisation is conducted in ParaView 5.11.2 (Ahrens et al., 2005; *ParaView* 2023).</figcaption></figure>

We use the parameters detailed in Table 1, which are generally representative of a brittle polymer such as poly(methyl methacrylate) (PMMA) or Homalite-100, which are classical model materials in experimental fracture mechanics, particularly for this geometry (Ramulu and Kobayashi, 1985).
<table><caption>Table 1: Values of the simulation parameters. Of these parameters, the mechanical parameters, $\theta$, $h_{crack}$, $\eta_L$, $\eta_q$ and a choice of the LCP solver are inherent to our method and so must be specified. Specifying $c_{cohesion}$ is substantially required in practice as despite our proof of solution in the general discrete setting, the concrete implementation in floating point arithmetic can occasionally pose difficulties for the feasibility of the problem. Setting a value of $c_{threshold}$ gives our method power to simulate the rigid branch of the cohesive law, if we so choose. The other parameter values are merely details of implementation, and it is perfectly possible to use the method without adaptive time-stepping and accepting the default parameter values of the chosen LCP solver.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td colspan="3">Mechanical Parameters</td></tr><tr><td>Young's modulus <i>E</i></td><td>32 &times; 10<sup>3</sup></td><td>MPa</td></tr><tr><td>Poisson's ratio <i>&nu;</i></td><td>0.2</td><td></td></tr><tr><td>Density <i>&rho;</i></td><td>0.00245</td><td>g/mm<sup>3</sup></td></tr><tr><td>Critical traction <i>&sigma;</i><sub>c</sub></td><td>12</td><td>MPa</td></tr><tr><td>Critical fracture energy <i>G</i><sub>c</sub></td><td>0.003</td><td>N/mm</td></tr><tr><td>Critical traction ratio <i>&gamma;</i></td><td>1</td><td></td></tr><tr><td>Newton's coefficient of restitution <i>e</i></td><td>0</td><td></td></tr><tr><td>Coulomb friction coefficient <i>&mu;</i></td><td>0 / 0.4</td><td></td></tr><tr><td colspan="3">Time step parameters</td></tr><tr><td><i>&theta;</i>-method parameter</td><td>0.5</td><td></td></tr><tr><td>Maximum time step before cracking <i>h</i><sub>max</sub></td><td>1 &times; 10<sup>-4</sup></td><td>ms</td></tr><tr><td>Maximum time step after cracking <i>h</i><sub>crack</sub></td><td>1 &times; 10<sup>-5</sup></td><td>ms</td></tr><tr><td colspan="3">Adaptive time-stepping parameters</td></tr><tr><td>Minimum time step threshold <i>c</i><sub>min</sub></td><td>1 &times; 10<sup>-4</sup></td><td></td></tr><tr><td>Division factor for unsuccessful resolution <i>c</i><sub>divide</sub></td><td>2</td><td></td></tr><tr><td>Multiplication factor for successful resolution <i>c</i><sub>multiply</sub></td><td>1.2</td><td></td></tr><tr><td>Number of successful resolutions to trigger multiplication <i>n</i><sub>success</sub></td><td>5</td><td></td></tr><tr><td colspan="3">Numerical parameters</td></tr><tr><td>Matrix regularisation factor <i>&eta;</i><sub>L</sub></td><td>1 &times; 10<sup>-10</sup></td><td></td></tr><tr><td>Vector regularisation factor <i>&eta;</i><sub>q</sub></td><td>1 &times; 10<sup>-8</sup></td><td></td></tr><tr><td>Cohesion threshold <i>c</i><sub>cohesion</sub></td><td>0.001</td><td></td></tr><tr><td>Insertion threshold <i>c</i><sub>threshold</sub></td><td>1</td><td></td></tr><tr><td>Excess decohesion threshold <i>h</i>&lambda;<sub>threshold</sub></td><td>0.1</td><td></td></tr><tr><td>LCP solver</td><td>SICONOS_LCP_LEMKE</td><td></td></tr><tr><td>LCP tolerance</td><td>1 &times; 10<sup>-12</sup></td><td></td></tr><tr><td>Maximum number of LCP iterations</td><td>1 &times; 10<sup>5</sup></td><td></td></tr><tr><td>LCP tolerance multiplication factor <i>c</i><sub>LCP</sub></td><td>1</td><td></td></tr></tbody></table>

The mechanical parameters are set to be identical to those in Nguyen (2014) (where those parameters exist in our model). The parameter $e$ exists in our model but not that used by Nguyen, and so we set $e = 0$, as this is known to be the value most suitable for finite element simulations (non-zero values will cause the elements to “jiggle” continuously in this case). Nguyen’s system is frictionless, so we conduct one simulation with $\mu = 0$. We are interested in demonstrating that our system can successfully resolve a frictional system, so we also conduct one simulation with $\mu = 0.4$, taking the value as being representative of the results for PMMA in Bouissou et al. (1998). We choose $\theta = 0.5$ to have the best possible energy conservation properties, and $h_{max} = 1 \times 10^{-4}$ ms to balance detecting the crack initiation accurately in time with numerical efficiency. $h_{crack} = 1 \times 10^{-5}$ ms is set to obtain very fine local resolution of the elasto-dynamics near the crack tip, although our method allows substantially larger crack steps while still satisfying the proof of the existence of solutions. Our minimum time step threshold $c_{min}$, division factor $c_{divide}$, multiplication factor $c_{multiply}$ and successful resolution number $n_{success}$ are chosen to balance an accurate resolution of the system at numerically challenging points (particularly when the system first bifurcates) with a reasonable time of simulation (noting that the system simulation time is very sensitive to $c_{min}$ in particular). The regularisation factors $\eta_L$ and $\eta_q$ are chosen to be as small as possible to minimise the deviation from the LCP obtained via algebraic manipulation in Appendix B.1. The cohesion threshold $c_{cohesion}$ is chosen to be small enough that for typical rates of decohesion and typical time step sizes in the simulation, the point would decohere in the next time step, while also being large enough that the solver does not have difficulty in respecting the double-sided constraint. The insertion threshold $c_{threshold}$ is set to 1 for numerical reasons, as a smaller value entails inserting more cohesive zones, which increases the time required to solve the LCP at each time step. The
excess decohesion threshold $h\lambda_{threshold}$ is set to allow a reasonable amount of decohesion to occur per time step (for larger time steps), while still ensuring that the decohesion process is finely resolved in time. The SICONOS\_LCP\_LEMKE solver is chosen due to the advantageous theoretical property of Lemke's algorithm (Lemke and Howson, Jr., 1964), namely the guarantee of finding a solution if one exists. This choice of algorithm also renders the values of LCP tolerance and tolerance multiplication factors moot, but they can be used for the other feasible algorithm choices (SICONOS\_LCP\_NEWTON\_FB\_FBLSA and SICONOS\_LCP\_AVI\_CAOFERRIS), so we still pass the listed values as dummy arguments to the simulation. In the case of the Lemke algorithm, the maximum iteration number only functions as a threshold to indicate an unsuccessful solve and does not affect the accuracy, but for the other algorithms it can determine the accuracy of the solution.

We set $\sigma_{applied} = 1.5$ MPa, applied as a step load (that is, the system is not allowed to first find the corresponding static equilibrium, but starts with all $u = 0$ and $v = 0$). We simulate the system for 0.06 ms, and write to .pvtu and .vtk files every $1 \times 10^{-5}$ ms. The end states of both the frictionless and frictional systems are shown in Figure 10.

<figure><img src="image_7.png" /><figcaption>Figure 10: The state of the system at the end of the simulation, with each mesh shown using a $1 \times$ warp factor (that is, the mesh in the image is displaced to match the displacements it possesses in the simulation). (a) The frictionless system, with the mesh coloured by total displacement and the points coloured by value of cohesion. (b) The frictional system, with the mesh coloured by total displacement and the points coloured by value of cohesion.</figcaption></figure>

We see in Figure 10 that in both cases the system bifurcates into two separate cracks, traceable by following the path of zero cohesion points in dark purple. The fact of this bifurcation is reasonably significant, as previous studies have shown an apparent inability for extrinsic cohesive zone models to support a crack bifurcation (Falk et al., 2001). Given the subsequent success of several authors using extrinsic cohesive zone models (Nguyen, 2014; Park et al., 2012; Rabczuk and Belytschko, 2004; Song et al., 2008) in demonstrating crack bifurcation with exactly the geometry and loading conditions detailed above, we infer that while intrinsic cohesive zone models may demonstrate crack bifurcation on coarser meshes, extrinsic models will demonstrate crack bifurcation provided the mesh is adequately refined. We also note several pertinent features of the crack pattern, in particular that one of the branches of the initial bifurcation eventually terminates. This initial bifurcation occurs at
about 0.028 ms into the simulation, matching the moment of bifurcation obtained by Nguyen (2014). Numerically, this point in time is also associated with very small time steps as the system resolves the spatial instability. It is also soon after this point that in spite of the symmetric mesh, the simulation becomes slightly asymmetric due to numerical error. Similar deviations from symmetry are observed in symmetric test specimens (such as those shown in Eid et al. (2023)), which highlights the sensitivity of global crack trajectories to very small perturbations in both models and the real world. This points to the well-known issue of mesh-sensitivity of the simulation, as the crack path is constrained to run along the edges of the elements. This issue can be alleviated somewhat by setting each side to have a different value of $\sigma_c$, using a modified Weibull distribution to enforce larger elements as being weaker (Zhou and Molinari, 2004). While for the sake of simplicity we have used a constant $\sigma_c$, our method straightforwardly allows for this strategy by treating all of the $\sigma_c$ and $\delta_{c,N}$ terms in (92) as vectors. However, for a solution that more closely approaches the ideal of mesh-independence, we must turn to more spatially sophisticated approaches such as the virtual element method (Marfia et al., 2022; Nguyen-Thanh et al., 2018), discontinuous Galerkin approaches (Nguyen, 2014; Versino et al., 2015), the extended finite element method (Moës and Belytschko, 2002) or pinwheel-based meshes (Papoulia et al., 2006). While we do not anticipate any particular difficulties in combining our approach with these frameworks, any such implementation remains beyond the scope of this article. Finally, we also observe that the frictional system is less prone to further crack branching, presumably as the additional apparent strength induced by the frictional resistance makes triggering the crack insertion condition less likely.

<figure><img src="image_8.png" /><figcaption>Figure 11: The energetic quantities of the frictionless and frictional edge-cracked block systems. (a) The total energy, and cumulated impact, friction and external work in the frictionless system. (b) The net difference of the work input with the stored and dissipated energies over the time step in the frictionless system. (c) The total energy, and cumulated impact, friction and external work in the frictional system. (d) The net difference of the work input with the stored and dissipated energies over the time step in the frictional system.</figcaption></figure>

In Figure 11, we calculate the total energy as given in §4.4 and cumulate the work terms in the standard way to obtain subfigures (a) and (c). We obtain the change in the energy balance over the time step by moving all terms to the right hand side of the equality in (127) and evaluating to obtain subfigures (b) and (d). In observing these systems, we see that the numerical scheme accurately accounts for the energetic input without introducing any numerical dissipation, hence the size of the net change in energy at each step is on the order of the floating point error. Further, at no point is the scheme anything other than energy-preserving (in the absence of impact and friction), or dissipative (in the presence of impact and friction), meaning that the practical results in implementation correspond exactly with the theoretical proofs obtained in §4.4. We note also that the frictionless system dissipates more energy than the frictional system, despite having fewer dissipative mechanisms. This occurs as the greater apparent strength of the frictional system causes fewer crack surfaces to be created, as it is the normal impact on these surfaces that is the main means of dissipation. In the frictional system, we can see that the frictional sliding does very little energy dissipation, with the normal impact being the main dissipative mechanism. These results highlight that
even in geometries and loading conditions that favour mode I cracking without sliding, inclusion of the friction is important
to accurately model the dissipation, and doing so is more conservative for design purposes, as on net less energy is dissipated.

6 Conclusions

In this paper, we used the principal of virtual power to establish the equilibrium and boundary conditions for a body featuring
a cohesive zone. We then postulated a free energy potential for the surface that allowed us to obtain nonsmooth state laws
for the system. These state laws describe a family of extrinsic cohesive zone models that do not feature the unload–reload
elasticity characteristic of those models that have the shifted intrinsic model structure. This mathematical feature guarantees
appropriate physical behaviour under complex, nonmonotonic loading conditions, where traditional extrinsic cohesive zone
models exhibit the same spurious behaviours as intrinsic models. We then specified a nonsmooth dissipative pseudo-potential
that allowed us to obtain a complete generalised family of constitutive models for the cohesive zone system for both normal and
tangential crack opening, and showed that in the case of frictional contact the model is dissipative (otherwise, it is irreversible
but nondissipative). By specifying a particular energetic potential, we obtained a model that gives a linear evolution of the
cohesion variable with increasing displacement jump across the crack faces.

We then spatially discretised our system and re-framed certain parts using differential measures. This allowed us to include
the nonsmooth contact dynamics formulation within our model. By using a θ-method to discretise the system in time, in com-
bination with the Frémond impact law for friction, we obtained a time-stepping scheme that we were able to write as a linear
complementarity problem. We then showed that given a sufficiently small time-step, the solution of the problem is guaranteed
to exist. The time-stepping scheme was then shown to be dissipative in the case of contact and friction, and energy-preserving
in their absence.

We then implemented the model numerically by solving the linear complementarity problem at each time-step, followed by a
global updating of the degrees of freedom not involved in the monolithic complementarity problem. By using the finite element
method, we were able to simulate a system of interest drawn from the literature, and show that our numerical method is able to
replicate the crack bifurcation behaviour that is observed experimentally and using other numerical methods. Empirically, the
method demonstrates excellent energetic conservation properties and is able to balance accuracy of solution with numerical
efficiency.

This work has demonstrated that the applications of techniques from convex analysis lead to physically correct and numeri-
cally efficient models that can encompass fracture, contact and friction. In formulating this work, we decided to remain in the
quasi-two-dimensional and small displacement framework so that we could continue to benefit from the linear complemen-
tarity formulation. The natural extension of the work is to pass to a fully three-dimensional model, however this will require
formulating the problem as a second-order cone complementarity problem, and is a substantial leap in complexity. Other
possible directions of research would be to explore other (nonlinear) forms of the cohesion law, extend the model to model
interface phenomena such as “rate-and-state” friction, or combine crack propagation with other physical phenomena such as
phase change.

Acknowledgements

The authors acknowledge the support of the Marie Skłodowska-Curie Actions program under the Horizon Europe research and
innovation framework program (Grant agreement ID 101064805 LEMMA). Views and opinions expressed are however those of
the authors only and do not necessarily reflect those of the European Union or Horizon Europe. Neither the European Union
nor the granting authority can be held responsible for them. The authors also express their gratitude to Thibault Ghesquière-
Diérickx, Guillaume Anciaux and Jean-François Molinari for their aid in integrating our model with the Akantu finite element
software. Any errors in doing so are our own. The authors also thank Franck Bourrier for insightful comments during the
development of the model.

Code and data availability

The codes required to run the simulations described in this paper are available in a [GitLab](https://gitlab.com/Sklo-dowska-Curie) repository archived on [Software Heritage](https://software-heritage.org), or alternatively on [Zenodo](). The data outputs of the codes are also available on [Zenodo]().

A Reformulation of a bi-dimensional Coulomb-like inclusion into a comple-
mentarity problem

Inspired by Stewart and Trinkle (1996), we give the results concerning the reformulation of Coulomb-like friction problem into
complementarity relations.
**Lemma 2.** The solution $y_N, y_T, x_N, x_T$ of the following inclusion

$$ \begin{cases} 0 \le y_N \perp x_N \ge 0, \\ -y_T \in y_N \operatorname{sgn}(x_T) = y_N \partial |x_T|, \end{cases} \qquad (129) $$

is contained in the solution of the following complementarity system

$$ \begin{cases} 0 \le y_N \perp x_N \ge 0, \\ 0 \le \hat{y}_T \perp 1 \lambda + D^\top x_T \ge 0, \\ 0 \le \lambda \perp y_N - 1^\top \hat{y}_T \ge 0, \end{cases} \qquad (130) $$

with $y_T = D\hat{y}_T$ and $D = [1, -1]$. Furthermore, we have $y_N|x_T| = y_N\lambda$.

**Proof** Let us consider the cases associated with

$$ \begin{cases} 0 \le \hat{y}_T \perp 1 \lambda + D^\top x_T \ge 0, \\ 0 \le \lambda \perp y_N - 1^\top \hat{y}_T \ge 0. \end{cases} \qquad (131) $$

1. $\lambda = 0$. In that case, the complementarity condition implies

$$ x_T = 0, \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N, \qquad (132) $$

and then $\lambda = |x_T|$. Since $\hat{y}_{T,1} \ge 0$, $\hat{y}_{T,2} \ge 0$, we also have

$$ \begin{cases} y_T = \hat{y}_{T,1} - \hat{y}_{T,2} = \hat{y}_{T,1} + \hat{y}_{T,2} - 2\hat{y}_{T,2} \le \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N, \\ -y_T = \hat{y}_{T,2} - \hat{y}_{T,1} = \hat{y}_{T,2} + \hat{y}_{T,1} - 2\hat{y}_{T,1} \le \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N, \end{cases} \qquad (133) $$

and we conclude $|y_T| \le y_N$.

2. $\lambda > 0, \hat{y}_{T,1} > 0, \hat{y}_{T,2} = 0$. In that case, the complementarity condition implies

$$ \lambda = -x_T > 0, x_T < 0, y_{T,1} = y_N, \qquad (134) $$

and then $y_T = y_N, x_T < 0, \lambda = |x_T|$.

3. $\lambda > 0, \hat{y}_{T,1} = 0, \hat{y}_{T,2} > 0$. In that case, the complementarity condition implies

$$ \lambda = x_T > 0, x_T > 0, y_{T,2} = y_N, \qquad (135) $$

and then $y_T = y_N, x_T > 0, \lambda = |x_T|$.

4. $\lambda > 0, \hat{y}_{T,1} > 0, \hat{y}_{T,2} > 0$. In that case, the complementarity condition implies

$$ \lambda = 0. \qquad (136) $$

This case is not possible.

5. $\lambda > 0, \hat{y}_{T,1} = 0, \hat{y}_{T,2} = 0$. In that case, the complementarity condition implies

$$ y_N = 0, \lambda \ge -x_T, \lambda \ge x_T. \qquad (137) $$

and then $y_N\lambda = y_N|x_T|$.

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

As mentioned in §4.2, the last item of the proof may lead to the addition of unwanted solutions. To remedy this, we propose adding a regularisation variable to the system that will very slightly increase the threshold of the sign function and eliminate artificial solutions.

**Lemma 3.** Let $\epsilon > 0$. The solution $y_N, y_T, x_N, x_T$ of the following inclusion

$$ \begin{cases} 0 \le y_N \perp x_N \ge 0, \\ -y_T \in (y_N + \epsilon) \operatorname{sgn}(x_T) = (y_N + \epsilon) \partial |x_T|, \end{cases} \qquad (138) $$

is given by solving the following complementarity system

$$ \begin{cases} 0 \le y_N \perp x_N \ge 0, \\ 0 \le \hat{y}_T \perp 1 \lambda + D^\top x_T \ge 0, \\ 0 \le \lambda \perp (y_N + \epsilon) 1^\top \hat{y}_T \ge 0, \end{cases} \qquad (139) $$

with $y_T = D\hat{y}_T$ and $D = [1, -1]$. Furthermore, we have $|x_T| = \lambda$.
**Proof** Let us consider the cases associated with

$$ \begin{cases} 0 \le \hat{y}_T \perp 1\lambda + D^\top x_T \ge 0, \\ 0 \le \lambda \perp (y_N + \epsilon) - 1^\top \hat{y}_T \ge 0. \end{cases} \qquad (140) $$

1. $\lambda = 0$. In that case, the complementarity condition implies

$$ x_T = 0, \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N + \epsilon, \qquad (141) $$

and then $\lambda = |x_T|$. Since $\hat{y}_{T,1} \ge 0$, $\hat{y}_{T,2} \ge 0$, we also have

$$ \begin{cases} y_T = \hat{y}_{T,1} - \hat{y}_{T,2} = \hat{y}_{T,1} + \hat{y}_{T,2} - 2\hat{y}_{T,2} \le \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N + \epsilon, \\ -y_T = \hat{y}_{T,2} - \hat{y}_{T,1} = \hat{y}_{T,2} + \hat{y}_{T,1} - 2\hat{y}_{T,1} \le \hat{y}_{T,1} + \hat{y}_{T,2} \le y_N + \epsilon, \end{cases} \qquad (142) $$

and we conclude $|y_T| \le y_N + \epsilon$.

2. $\lambda > 0, \hat{y}_{T,1} > 0, \hat{y}_{T,2} = 0$. In that case, the complementarity condition implies

$$ \lambda = -x_T > 0, x_T < 0, y_{T,1} = y_N + \epsilon, \qquad (143) $$

and then $y_T = y_N + \epsilon, x_T < 0, \lambda = |x_T|$.

3. $\lambda > 0, \hat{y}_{T,1} = 0, \hat{y}_{T,2} > 0$. In that case, the complementarity condition implies

$$ \lambda = x_T > 0, x_T > 0, y_{T,2} = y_N + \epsilon, \qquad (144) $$

and then $y_T = y_N + \epsilon, x_T > 0, \lambda = |x_T|$.

4. $\lambda > 0, \hat{y}_{T,1} > 0, \hat{y}_{T,2} > 0$. In that case, the complementarity condition implies

$$ \lambda = 0. \qquad (145) $$

This case is not possible.

5. $\lambda > 0, \hat{y}_{T,1} = 0, \hat{y}_{T,2} = 0$. In that case, the complementarity condition implies

$$ y_N = -\epsilon. \qquad (146) $$

This case is not possible. □

# B Linear complementarity problem development

Here, we detail the manipulations necessary to arrive at the LCPs detailed in this paper.

## B.1 Dynamic linear complementarity problem development

First, we need to express each components of the vector $w$ in terms of the components of the vector $z$ using the time-discretisation (86). Having chosen the LCP variables in this way, the expression for $w_1$ is trivial:

$$ w_1 = S^{-1} (S\beta_k - S\beta_{k+\theta}) = -S^{-1}z_2 + \beta_k. \qquad (147) $$

Likewise, the expression for $w_7$:

$$ w_7 = \sigma_c \gamma z_2 - 1^\top z_6, \qquad (148) $$

and for $w_5$:

$$ w_5 = \mu z_3 - 1^\top z_4. \qquad (149) $$

To continue to express the components of $w$, we will use the following relations

$$ v_{k+\theta} - v_k = \theta(v_{k+1} - v_k) \text{ and } u_{k+\theta} = u_k + h\theta v_{k+\theta}. \qquad (150) $$

We start the development of the expressions for the more involved terms by expanding the first line of (86) multiplied by $\theta$. We arrive at

$$ \begin{aligned} M(v_{k+\theta} - v_k) + h\theta K (u_k + h\theta v_{k+\theta}) = & h\theta F_{k+\theta} - h\theta \sigma_c H_N^\top S \beta_{k+\theta} + \theta \bar{H}_N^\top p_{N,k,k+1} \\ & + \theta \bar{H}_T^\top D \hat{p}_{T,k,k+1} + h\theta H_T^\top D \hat{r}_{T,k+\theta}^r. \end{aligned} \qquad (151) $$
We denote the augmented mass matrix as $\hat{M} = M + h^2\theta^2 K$ and the free impulse (without the contribution of the cohesive zone model) $\hat{i}_{k,k+1} = Mv_k - h\theta Ku_k + h\theta F_{k+\theta}$. When necessary, we modify the augmented mass matrix $\hat{M}$ and the free impulse $\hat{i}$ to take into account Dirichlet boundary conditions. Thus, the velocities $v_{k+\theta}$ can be determined by

$$v_{k+\theta} = \hat{M}^{-1} \left( \hat{i}_{k,k+1} - h\theta \sigma_c H_N^\top S \beta_{k+\theta} + \theta \bar{H}_N^\top p_{N,k,k+1} + \theta \bar{H}_T^\top D \hat{p}_{T,k,k+1} + h\theta H_T^\top D \hat{r}_{T,k+\theta} \right). \quad (152)$$

In terms of the LCP variables, we obtain

$$v_{k+\theta} = \hat{M}^{-1} \left( \hat{i}_{k,k+1} - h\theta \sigma_c H_N^\top z_2 + \theta \bar{H}_N^\top z_3 + \theta \bar{H}_T^\top D z_4 + h\theta H_T^\top D z_6 \right). \quad (153)$$

We obtain the normal relative velocities by applying $\bar{H}_N$ to (153):

$$\begin{aligned} v_{N,k+\theta} &= \bar{H}_N v_{k+\theta}, \\ &= -h\theta \sigma_c \bar{H}_N \hat{M}^{-1} H_N^\top z_2 + \theta \bar{H}_N \hat{M}^{-1} \bar{H}_N^\top z_3 + \theta \bar{H}_N \hat{M}^{-1} \bar{H}_T^\top D z_4 + h\theta \bar{H}_N \hat{M}^{-1} H_T^\top D z_6 + \bar{H}_N \hat{M}^{-1} \hat{i}_{k,k+1}, \\ &= -h\theta \sigma_c V_{NN} z_2 + \theta W_{NN} z_3 + \theta W_{NT} D z_4 + h\theta V_{NT} D z_6 + \bar{H}_N \hat{M}^{-1} \hat{i}_{k,k+1}, \end{aligned} \quad (154)$$

and the tangential relative velocities by applying $\bar{H}_T$ to (153):

$$\begin{aligned} v_{T,k+\theta} &= \bar{H}_T v_{k+\theta}, \\ &= -h\theta \sigma_c \bar{H}_T \hat{M}^{-1} H_N^\top z_2 + \theta \bar{H}_T \hat{M}^{-1} \bar{H}_N^\top z_3 + \theta \bar{H}_T \hat{M}^{-1} \bar{H}_T^\top D z_4 + h\theta \bar{H}_T \hat{M}^{-1} H_T^\top D z_6 + \bar{H}_T \hat{M}^{-1} \hat{i}_{k,k+1}, \\ &= -h\theta \sigma_c V_{TN} z_2 + \theta W_{TN} z_3 + \theta W_{TT} D z_4 + h\theta V_{TT} D z_6 + \bar{H}_T \hat{M}^{-1} \hat{i}_{k,k+1}, \end{aligned} \quad (155)$$

where $W_{NN} = \bar{H}_N \hat{M}^{-1} \bar{H}_N^\top$ and $W_{TT} = \bar{H}_T \hat{M}^{-1} \bar{H}_T^\top$ are the Delassus matrices for the degrees of freedom involved in the purely normal and purely tangential contact problems, while $W_{NT} = \bar{H}_N \hat{M}^{-1} \bar{H}_N^\top$ and $W_{TN} = \bar{H}_T \hat{M}^{-1} \bar{H}_T^\top$ are the Delassus matrices for the interaction terms between the normal and tangential parts of the contact problem. The matrices $V_{NN} = \bar{H}_N \hat{M}^{-1} H_N^\top$, $V_{TT} = \bar{H}_T \hat{M}^{-1} H_T^\top$, $V_{NT} = \bar{H}_N \hat{M}^{-1} H_N^\top$ and $V_{TN} = \bar{H}_T \hat{M}^{-1} H_T^\top$ capture the effect of the cohesive terms on the contact and friction degrees of freedom. From (154) and (155) we can straightforwardly get first $w_3$:

$$\begin{aligned} w_3 &= \theta v_{N,k+\theta} + \theta(\theta(1+e)-1)v_{N,k}, \\ &= -h\theta^2 \sigma_c V_{NN} z_2 + \theta^2 W_{NN} z_3 + \theta^2 W_{NT} D z_4 + h\theta^2 V_{NT} D z_6 + \theta \bar{H}_N \hat{M}^{-1} \hat{i}_{k,k+1} + \theta(\theta(1+e)-1) \bar{H}_N v_k, \end{aligned} \quad (156)$$

and then $w_4$:

$$\begin{aligned} w_4 &= \mathbb{1}\theta\zeta_{k+\theta} + \theta D^\top v_{T,k+\theta}, \\ &= -h\theta^2 \sigma_c D^\top V_{TN} z_2 + \theta^2 D^\top W_{TN} z_3 + \theta^2 D^\top W_{TT} D z_4 + \mathbb{1}z_5 + h\theta^2 D^\top V_{TT} D z_6 + \theta D^\top \bar{H}_T \hat{M}^{-1} \hat{i}_{k,k+1}. \end{aligned} \quad (157)$$

In the same way as for the velocity, we can expand the expressions for the normal and tangential relative displacements. Firstly in the normal direction:

$$\begin{aligned} u_{N,k+\theta} &= H_N u_{k+\theta} + b_{N,k+\theta}, \\ &= H_N (u_k + h\theta v_{k+\theta}) + b_{N,k+\theta}, \\ &= -h^2\theta^2 \sigma_c H_N \hat{M}^{-1} H_N^\top z_2 + h\theta^2 H_N \hat{M}^{-1} \bar{H}_N^\top z_3 + h\theta^2 H_N \hat{M}^{-1} \bar{H}_T^\top D z_4 + h^2\theta^2 H_N \hat{M}^{-1} H_T^\top D z_6 + q_{u_N}, \\ &= -h^2\theta^2 \sigma_c U_{NN} z_2 + h\theta^2 V_{NN}^\top z_3 + h\theta^2 V_{TN}^\top D z_4 + h^2\theta^2 U_{NT} D z_6 + q_{u_N}, \end{aligned} \quad (158)$$

and then in the tangential direction:

$$\begin{aligned} u_{T,k+\theta} &= H_T u_{k+\theta} + b_{T,k+\theta}, \\ &= H_T (u_k + h\theta v_{k+\theta}) + b_{T,k+\theta}, \\ &= -h^2\theta^2 \sigma_c H_T \hat{M}^{-1} H_N^\top z_2 + h\theta^2 H_T \hat{M}^{-1} \bar{H}_N^\top z_3 + h\theta^2 H_T \hat{M}^{-1} \bar{H}_T^\top D z_4 + h^2\theta^2 H_T \hat{M}^{-1} H_T^\top D z_6 + q_{u_T}, \\ &= -h^2\theta^2 \sigma_c U_{TN} z_2 + h\theta^2 V_{NT}^\top z_3 + h\theta^2 V_{TT}^\top D z_4 + h^2\theta^2 U_{TT} D z_6 + q_{u_T}, \end{aligned} \quad (159)$$

where for compactness we have introduced the terms $q_{u_N}$ and $q_{u_T}$ given by

$$q_{u_N} = h\theta H_N \hat{M}^{-1} \hat{i}_{k,k+1} + H_N u_k + b_{N,k+\theta}, \quad (160)$$

$$q_{u_T} = h\theta H_T \hat{M}^{-1} \hat{i}_{k,k+1} + H_T u_k + b_{T,k+\theta}, \quad (161)$$

and where $U_{NN} = H_N \hat{M}^{-1} H_N^\top$ and $U_{TT} = H_T \hat{M}^{-1} H_T^\top$ are the complete Delassus matrices for the normal and tangential degrees of freedom, while $U_{NT} = H_N \hat{M}^{-1} H_N^\top$ and $U_{TN} = H_T \hat{M}^{-1} H_T^\top$ are the Delassus matrices for the interaction terms between the two.
Equation (159) leads directly to one for $w_6$:

$$ 
\begin{aligned}
w_6 &= \mathbf{1}\chi_{k+\theta} + D^\top u_{\mathbf{T},k+\theta} \\
&= -h^2\theta^2\sigma_c D^\top U_{\mathbf{T}\mathbf{N}} z_2 + h\theta^2 D^\top V_{\mathbf{N}\mathbf{T}}^\top z_3 + h\theta^2 D^\top V_{\mathbf{T}\mathbf{T}}^\top D z_4 + h^2\theta^2 D^\top U_{\mathbf{T}\mathbf{T}} D z_6 + \mathbf{1}z_7 + D^\top q_{u_{\mathbf{T}}}.
\end{aligned} 
\quad (162) $$

Finally, we consider our expression for the slack variable $\xi_{k+\theta} = A^r_{k+\theta} + \sigma_c \delta_{c,N}(\beta_{k+\theta} - \mathbf{1}) + \sigma_c u_{N,k+\theta} + \sigma_c \gamma |u_{\mathbf{T},k+\theta}|$ where $\mathbf{1}$ represents a vector of ones, using the expressions of $u_{N,k+\theta}$ and $u_{\mathbf{T},k+\theta}$ for all cohesive points. As detailed in §4.2, we make the choice to replace $|u_{\mathbf{T},k+\theta}|$ with $\chi_{k+\theta}$ in order to obtain an LCP, at the cost of introducing non-uniqueness into the system at $u_{\mathbf{T}} = 0$. Having made this choice, we obtain an expression for $w_2$:

$$ 
\begin{aligned}
w_2 &= \sigma_c \delta_{c,N}(\beta_{k+\theta} - \mathbf{1})) + \sigma_c \gamma \chi_{k+\theta} + A^r_{k+\theta} + \sigma_c u_{N,k+\theta} \\
&= \sigma_c \delta_{c,N}(S^{-1}z_2 - \mathbf{1})) + \sigma_c \gamma z_7 + S^{-1}z_1 + \sigma_c [-h^2\theta^2\sigma_c U_{NN}z_2 + h\theta^2 V_{NN}^\top z_3 + h\theta^2 V_{\mathbf{T}\mathbf{T}}^\top D z_4 + h^2\theta^2 U_{\mathbf{T}\mathbf{T}} D z_6 + q_{u_N}] \\
&= S^{-1}z_1 + \sigma_c (\delta_{c,N}S^{-1} - h^2\theta^2\sigma_c U_{NN})z_2 + h\theta^2\sigma_c V_{NN}^\top z_3 + h\theta^2\sigma_c V_{\mathbf{T}\mathbf{T}}^\top D z_4 + h^2\theta^2\sigma_c U_{\mathbf{T}\mathbf{T}} D z_6 + \sigma_c \gamma z_7 + \sigma_c q_{u_N} - \sigma_c \delta_{c,N}\mathbf{1}
\end{aligned} 
\quad (163) $$

At this point we have obtained the complete set of expressions necessary to construct the LCP. We then add a regularisation term $\eta_L$ along all the main diagonal entries of the L matrix, as well as another regularisation term $\eta_q$ to the last entry of the q vector, which together have the effect of facilitating the operation of the pivot solver and suppressing the non-physical solution at $u_{\mathbf{T}} = 0$, and giving us the LCP as presented in (92).

## B.2 Quasi-static linear complementarity problem development

As noted above, in the quasi-static system we will not make any distinction between $H$ and $\bar{H}$, and hence all the matrices will be $U$. The global elastic-cohesive-frictional-contact problem is thus:

$$ 
\left\{
\begin{array}{l}
K u_{k+\theta} = F_{k+\theta} - \sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r, \\
S r_{T,k+\theta}^r = D \hat{r}_{T,k+\theta}^r, \\
S r_{T,k+\theta}^{ir} = D \hat{r}_{T,k+\theta}^{ir}, \\
u_{k+1} = u_k + h v_{k+\theta}, \\
u_{N,k+\theta} = H_N u_{k+\theta} + b_{N,k+\theta}, \quad u_{\mathbf{T},k+\theta} = H_{\mathbf{T}} u_{k+\theta} + b_{\mathbf{T},k+\theta}, \quad v_{\mathbf{T},k+\theta} = H_{\mathbf{T}} v_{k+\theta}, \\
\beta_{k+1} = \beta_k - h \lambda_{k+\theta}, \\
\sigma_c \delta_{c,N} (\beta_{k+\theta} - \mathbf{1}) + \sigma_c u_{N,k+\theta} + \sigma_c \gamma |u_{\mathbf{T},k+\theta}| + A^r_{k+\theta} = \xi_{k+\theta}, \\
r_{N,k+\theta}^c = r_{N,k+\theta} + \sigma_c \beta_{k+\theta}, \\
0 \le S A^r_{k+\theta} \perp \lambda_{k+\theta} \ge 0, \\
0 \le S \beta_{k+\theta} \perp \xi_{k+\theta} \ge 0, \\
0 \le S r_{N,k+\theta}^c \perp u_{N,k+\theta} \ge 0, \\
0 \le \hat{r}_{T,k+\theta}^{ir} \perp \mathbf{1} \zeta_{k+\theta} + D^\top v_{\mathbf{T},k+\theta} \ge 0, \\
0 \le \zeta_{k+\theta} \perp \mu S r_{N,k+\theta} - \mathbf{1}^\top \hat{r}_{T,k+\theta}^{ir} \ge 0, \\
0 \le \hat{r}_{T,k+\theta}^r \perp \mathbf{1} \chi_{k+\theta} + D^\top u_{\mathbf{T},k+\theta} \ge 0, \\
0 \le \chi_{k+\theta} \perp S \beta_{k+\theta} \sigma_c \gamma - \mathbf{1}^\top \hat{r}_{T,k+\theta}^r \ge 0.
\end{array}
\right.
\quad (164) $$

We establish our complementarity vectors $w$ and $z$ as

$$ w_{\text{stat}} = \begin{bmatrix} h\theta\lambda_{k+\theta} \\ \xi_{k+\theta} \\ u_{N,k+\theta} \\ \mathbf{1}h\theta\zeta_{k+\theta} + D^\top h\theta v_{\mathbf{T},k+\theta} \\ \mu S r_{N,k+\theta}^c - \mathbf{1}^\top \hat{r}_{T,k+\theta}^{ir} \\ \mathbf{1}\chi_{k+\theta} + D^\top u_{\mathbf{T},k+\theta} \\ \sigma_c \gamma S \beta_{k+\theta} - \mathbf{1}^\top \hat{r}_{T,k+\theta}^r \end{bmatrix}, \quad z_{\text{stat}} = \begin{bmatrix} S A^r_{k+\theta} \\ S \beta_{k+\theta} \\ S r_{N,k+\theta}^c \\ \hat{r}_{T,k+\theta}^{ir} \\ h\theta\zeta_{k+\theta} \\ \hat{r}_{T,k+\theta}^r \\ \chi_{k+\theta} \end{bmatrix}, \quad (165) $$

where we consider $h\theta\lambda_{k+\theta}$, $h\theta v_{\mathbf{T},k+\theta}$ and $h\theta\zeta_{k+\theta}$ to avoid an ill-conditioned matrix as $h \to 0$. We can begin rearranging our equations to form the LCP proper. Once again, we obtain trivially $w_1 = S^{-1}(S\beta_k - S\beta_{k+\theta}) = -S^{-1}z_2 + \beta_k$, $w_5 = \mu z_3 - \mathbf{1}^\top z_4$
and $w_7 = \sigma_c \gamma z_2 - 1^\top z_6$. Then, assuming a modification of K and F to take into account the boundary conditions, we have

$$K u_{k+\theta} = F_{k+\theta} - \sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r,$$

$$K(u_k + h\theta v_{k+\theta}) = F_{k+\theta} - \sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r,$$

$$h\theta v_{k+\theta} = K^{-1} \left( -\sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r + F_{k+\theta} \right) - u_k. \quad (166)$$

Now, applying the tangential selection matrix $H_T$ to $h\theta v_{k+\theta}$ we obtain

$$\begin{aligned}
h\theta v_{T,k+\theta} &= -\sigma_c H_T K^{-1} H_N^\top S \beta_{k+\theta} + H_T K^{-1} H_N^\top S r_{N,k+\theta}^c + H_T K^{-1} H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T K^{-1} H_T^\top D \hat{r}_{T,k+\theta}^r \\
&\quad + H_T K^{-1} F_{k+\theta} - H_T u_k, \\
&= -\sigma_c U_{TN} S \beta_{k+\theta} + U_{TN} S r_{N,k+\theta}^c + U_{TT} D \hat{r}_{T,k+\theta}^\text{ir} + U_{TT} D \hat{r}_{T,k+\theta}^r + H_T K^{-1} F_{k+\theta} - H_T u_k, \\
&= -\sigma_c U_{TN} z_2 + U_{TN} z_3 + U_{TT} D z_4 + U_{TT} D z_6 + H_T K^{-1} F_{k+\theta} - H_T u_k.
\end{aligned} \quad (167)$$

This directly gives us $w_4$:

$$w_4 = -\sigma_c D^\top U_{TN} z_2 + D^\top U_{TN} z_3 + D^\top U_{TT} D z_4 + 1 z_5 + D^\top U_{TT} D z_6 + D^\top H_T K^{-1} F_{k+\theta} - D^\top H_T u_k. \quad (168)$$

Now, the normal relative displacement is given by

$$\begin{aligned}
u_{N,k+\theta} &= H_N u_{k+\theta} + b_{N,k+\theta}, \\
&= H_N h\theta v_{k+\theta} + H_N u_k + b_{N,k+\theta}, \\
&= H_N K^{-1} \left( -\sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r + F_{k+\theta} \right) - H_N u_k + H_N u_k + b_{N,k+\theta}, \\
&= -\sigma_c U_{NN} S \beta_{k+\theta} + U_{NN} S r_{N,k+\theta}^c + U_{NT} D \hat{r}_{T,k+\theta}^\text{ir} + U_{NT} D \hat{r}_{T,k+\theta}^r + H_N K^{-1} F_{k+\theta} + b_{N,k+\theta}, \\
&= -\sigma_c U_{NN} z_2 + U_{NN} z_3 + U_{NT} D z_4 + U_{NT} D z_6 + H_N K^{-1} F_{k+\theta} + b_{N,k+\theta},
\end{aligned} \quad (169)$$

which gives us $w_3$ directly. Similarly, the tangential relative displacement is given by

$$\begin{aligned}
u_{T,k+\theta} &= H_T u_{k+\theta} + b_{T,k+\theta}, \\
&= H_T h\theta v_{k+\theta} + H_T u_k + b_{T,k+\theta}, \\
&= H_T K^{-1} \left( -\sigma_c H_N^\top S \beta_{k+\theta} + H_N^\top S r_{N,k+\theta}^c + H_T^\top D \hat{r}_{T,k+\theta}^\text{ir} + H_T^\top D \hat{r}_{T,k+\theta}^r + F_{k+\theta} \right) - H_T u_k + H_T u_k + b_{T,k+\theta}, \\
&= -\sigma_c U_{TN} S \beta_{k+\theta} + U_{TN} S r_{N,k+\theta}^c + U_{TT} D \hat{r}_{T,k+\theta}^\text{ir} + U_{TT} D \hat{r}_{T,k+\theta}^r + H_T K^{-1} F_{k+\theta} + b_{T,k+\theta}, \\
&= -\sigma_c U_{TN} z_2 + U_{TN} z_3 + U_{TT} D z_4 + U_{TT} D z_6 + H_T K^{-1} F_{k+\theta} + b_{T,k+\theta}.
\end{aligned} \quad (170)$$

From, this we easily obtain $w_6$:

$$\begin{aligned}
w_6 &= 1 \chi_{k+\theta} + D^\top u_{T,k+\theta}, \\
&= -\sigma_c D^\top U_{TN} S \beta_{k+\theta} + D^\top U_{TN} S r_{N,k+\theta}^c + D^\top U_{TT} D \hat{r}_{T,k+\theta}^\text{ir} + D^\top U_{TT} D \hat{r}_{T,k+\theta}^r + 1 \chi_{k+\theta} + D^\top H_T K^{-1} F_{k+\theta} + D^\top b_{T,k+\theta}, \\
&= -\sigma_c D^\top U_{TN} z_2 + D^\top U_{TN} z_3 + D^\top U_{TT} D z_4 + D^\top U_{TT} D z_6 + 1 z_7 + D^\top H_T K^{-1} F_{k+\theta} + D^\top b_{T,k+\theta}.
\end{aligned} \quad (171)$$

Finally, we require an expression for $\xi_{k+\theta}$. We make the same replacement ($|u_{T,k+\theta}| \to \chi_{k+\theta}$) as in the dynamic case, and obtain

$$\begin{aligned}
\xi_{k+\theta} &= A^T_{k+\theta} + \sigma_c \delta_{c,N} (\beta_{k+\theta} - \mathbf{1}) + \sigma_c u_{N,k+\theta} + \sigma_c \gamma \chi_{k+\theta}, \\
&= S^{-1} S A^T_{k+\theta} + \sigma_c \delta_{c,N} S^{-1} S \beta_{k+\theta} + \sigma_c \left( -\sigma_c U_{NN} S \beta_{k+\theta} + U_{NN} S r_{N,k+\theta}^c + U_{NT} D \hat{r}_{T,k+\theta}^\text{ir} + U_{NT} D \hat{r}_{T,k+\theta}^r + H_N K^{-1} F_{k+\theta} + b_{N,k+\theta} \right. \\
&\qquad \left. + \sigma_c \gamma \chi_{k+\theta} - \sigma_c \delta_{c,N} \mathbf{1} \right), \\
&= S^{-1} z_1 + \sigma_c \left( \delta_{c,N} S^{-1} - \sigma_c U_{NN} \right) z_2 + \sigma_c U_{NN} z_3 + \sigma_c U_{NT} D z_4 + \sigma_c U_{NT} D z_6 + \sigma_c \gamma z_7 + \sigma_c \left( H_N K^{-1} F_{k+\theta} + b_{N,k+\theta} - \delta_{c,N} \mathbf{1} \right),
\end{aligned} \quad (172)$$

which gives us $w_2$ directly. With all of the equations specified, we can write the LCP presented in (128).

## C Finite element implementation details

After having read in the finite element mesh and generated the appropriate structural stiffness matrix K and consistent mass matrix M, we set an initial time step size that is chosen for efficiency and does not necessarily resolve the very fast dynamics of the system (that is to say that elastic waves can travel over more than one element per time step). We then solve the system without any cohesive zones inserted by obtaining the velocities by (152) (using SciPy's sparse matrix library (Virtanen et al.,
2020) to solve the linear system with an LU decomposition) and then updating the velocities and displacements at the end of the time step following the principles of the θ-method. The nodal variables v, u, F and the derived quantities at the Gauss points (ε and σ) are written into a .ptvu file.

At the end of each time step, Akantu checks whether a cohesive zone needs to be inserted between two elements. This is done by calculating the stress at each Gauss point, obtaining the interpolation functions on the basis of the “Gauss element” formed by the Gauss points, and performing an extrapolation to the quadrature points (Felippa, 2004; Hinton and Campbell, 1974). It should be noted that for this purpose Akantu always uses at least two interpolation points along the edge of the element, even for linear triangular elements. In this work, we use this type of element, and so the stress is calculated at a single Gauss point and should be constant within the element and identical at the two quadrature points lying on the interface (at ±1/√3 in the internal isoparametric coordinate running along the interface). Nevertheless, the method described is generic, and can work straightforwardly with higher order elements such as the quadratic triangle (T6).

Once the stresses are obtained at the quadrature points, we calculate the traction τ at each point using the stress value calculated from each element and the normal at the quadrature point (that is to say in our case of two elements with two quadrature points, we calculate four sets of tractions). For each set of tractions we calculate the normal traction τᵢ = σ · n · n, the compressive component of the normal traction τcomp = |max(0, −τᵢ)| (so for a tensile stress, this value will be zero), and from this the frictional resistance is calculated by τfriction = μτcomp. Then, the tangential traction τt = σ · n · t is calculated and if

$$|\tau_{\text{t}}| \ge \begin{cases} -\gamma (\tau_i - c_{\text{threshold}}\sigma) & \text{if } \tau_{\text{comp}} = 0, \text{ or} \\ \mu\tau_{\text{comp}} + c_{\text{threshold}}\sigma_c\gamma & \text{if } \tau_{\text{comp}} > 0, \end{cases} \quad (173)$$

is true, the quadrature point indicates that it fulfils the cohesive zone insertion criterion. The parameter cthreshold ∈ [0, 1] is a numerical parameter that controls at what fraction of the critical traction a cohesive zone should be inserted.

**Remark 8.** Conventionally in extrinsic cohesive zone implementations, this parameter exists only implicitly as it must be equal to one. While we also set the value of the parameter to 1 for reasons of numerical efficiency, in principle our method is perfectly able to set cthreshold < 1 as our formulation is able to simulate the rigid branch of the cohesive law, which other implementations are not capable of.

If all of the quadrature points indicate that they fulfil the cohesive zone insertion criterion a cohesive zone is inserted along the element interface, and one or more of the corresponding nodes are duplicated.

**Remark 9.** Our choice to use a condition that all the quadrature points on the interface must fulfil the insertion criterion separately using the stresses from both elements is merely an implementation detail, and other choices are possible. For instance, we could also use a criterion that if any of the quadrature points fulfil the insertion criterion using the stresses from either element then a cohesive zone should be inserted, or we could first calculate the average stress using both elements, and then calculate the corresponding tractions to check the insertion condition against (which is the most common approach). Our model and methodology is perfectly coherent with any of these approaches, and once again the choice is made for reasons of numerical efficiency.

At this point, the code indicates that it is passing to the crack simulation stage. Akantu duplicates the nodes and rebuilds the K and M matrices, taking this duplication into account, but remaining in the total Lagrangian finite element framework (that is, the matrices are calculated with respect to the nodes’ initial positions, not their current positions). A new maximum time step size hcrack is chosen, where we select one that is sufficiently small to fully resolve the very fast elastodynamics of the system. For practical purposes, this is smaller than what is required to fulfil Assumption 1, although we nevertheless check this condition each time we build the corresponding matrices. Then, we reconstruct M⁻¹, reimposing our boundary conditions, as well as the selection and interpolation matrices required for the cohesive zone model. In constructing the H and H matrices, we allow the user to choose whether to diverge very slightly from the otherwise linear finite element framework, and calculate the normal and tangential vectors based on the current positions of the corresponding nodes at the instant of creation, rather than their original positions. This acts as a small regularisation of the model and enables much larger time steps to be used at points at which a crack bifurcates. Then, the LCP matrix L is calculated.

At this point, the program enters a new time-stepping loop and at each step, the LCP vector q is calculated. The vectors wwarm and zwarm which contain the values of the appropriate variables from the previous time step are created to give the LCP a “warm start” to aid convergence. Then, the cohesive zones are characterised to determine whether they are intact and/or in contact. For those zones that are no longer intact, their entries corresponding to the first, second, sixth and seventh lines of (91) are removed from consideration, while for those cohesive zones that are not in contact, their entries corresponding to the third, fourth and fifth lines of (91) are removed from consideration. Then, w and z are created by selecting from wwarm and zwarm only those entries that are still under consideration, and a corresponding reduced L and q are created. This operation is performed with the intention of obtaining the smallest possible LCP to solve, which aids numerical performance. Then, the reduced LCP is solved with Siconos.

Depending on the solution of the LCP, three paths are possible. In the first path, the LCP is successfully solved (that is to say, Siconos does not indicate that the LCP solution algorithm failed to converge). We then check a condition on the size of
$h\lambda_{k+1}$, ensuring that it is less than some value $h\lambda_{threshold}$ to ensure that the system does not decohere too rapidly. Providing that this condition is fulfilled, the solution is post-processed to obtain the values of the state variables at step $k + 1$, which are then stored and written to file as appropriate. Then, the cohesive zone insertion condition is checked once again at the end of the step. The time is advanced, and in the event that a new cohesive zone needs to be inserted, the steps to reconstruct the appropriate matrices are repeated. In both the case of no new cohesive zones or new cohesive zones being inserted, the system is solved once again as described above. In the event that the current time step is smaller than $h_{crack}$, and the system has been successfully solved a certain number of times in succession (with the exact number $n_{success}$ to be specified by the user), the time step size is multiplied by a user-specified factor $c_{multiply}$ at the end of the solution step (unless the resultant $h$ would be greater than $h_{crack}$, in which case it is set to that value), and at the start of the next time step, the corresponding matrices are reconstructed using the new time step size, and using the updated node positions in the case of the $H$ and $\bar{H}$ matrices (if that option is selected by the user).

In the second path, the LCP indicates a successful solution, but the condition on the size of $h\lambda_{k+1}$ fails. In this case, we do not conduct any post-processing, write any variables or advance the time. We instead divide the time step by a user-determined factor $c_{divide}$, and reconstruct the appropriate matrices and attempt the solution again. In the event that the time step size falls below a user-specified threshold $c_{min} \times h_{crack}$, rather than further dividing the time step, it is instead reset to $h_{crack}$ and the matrices are once again reconstructed and a new solution is attempted.

In the third path, the LCP indicates an unsuccessful solution, namely that Siconos could not converge. In this case, we follow the same procedure as for the second path, performing no post-processing or system updating, adapting the time step, reconstructing the relevant matrices and attempting the solution again.

The solution procedure continues until the maximum simulation time has been reached, at which point the code terminates.

## References

Acary, V. (May 2016). "Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact". In: *Journal of Applied Mathematics and Mechanics* / *Zeitschrift für Angewandte Mathematik und Mechanik* 96.5, pp. 585–603. DOI: 10.1002/zamm.201400231. URL: https://hal.inria.fr/hal-01235240.

Acary, V. and B. Brogliato (2008). *Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics*. Ed. by F. Pfeiffer and P. Wriggers. Vol. 35. Lecture Notes in Applied and Computational Mechanics. Berlin, Heidelberg: Springer. ISBN: 978-3-540-75391-9 978-3-540-75392-6. DOI: 10.1007/978-3-540-75392-6.

Acary, V. and N.A. Collins-Craft (June 30, 2025). "On the Moreau-Jean scheme with the Frémond impact law: energy conservation and dissipation properties for elastodynamics with contact, impact and friction". In: *Journal of Theoretical, Computational and Applied Mechanics*, pp. 1–30. DOI: 10.46298/jtcam.13480. URL: https://jtcam.episciences.org/13480.

Acary, V. and Y. Monerie (2006). *Nonsmooth fracture dynamics using a cohesive zone approach*. English. Research Report RR-6032. Inria, p. 56. URL: http://hal.inria.fr/inria-00110560/en/.

[SW] Acary, V. et al., *Siconos* version 4.5.0, July 24, 2024. Inria, Humboldt-Universität zu Berlin, and CNRS. LIC: Apache-2.0. VCS: https://gricad gitlab.univ-grenoble-alpes.fr/nonsmooth/siconos, SWHID: <swh:1:snp:44fe2b21cc16c6e6976e41c9a8efd57fd94bcc7f; origin= https://github.com/siconos/siconos>.

Ahrens, J., B. Geveci, and C. Law (2005). "36 - ParaView: An End-User Tool for Large-Data Visualization". In: *Visualization Handbook*. Ed. by C.D. Hansen and C.R. Johnson. Burlington: Butterworth-Heinemann, pp. 717–731. ISBN: 978-0-12-387582-2. DOI: https://doi.org/10.1016/B978-012387582-2/50038-1. URL: https://www.sciencedirect.com/science/article/pii/B9780123875822500381.

Armero, F. and C. Linder (Dec. 1, 2009). "Numerical Simulation of Dynamic Fracture Using Finite Elements with Embedded Discontinuities". In: *International Journal of Fracture* 160.2, pp. 119–141. ISSN: 1573-2673. DOI: 10.1007/s10704-009-9413-9. URL: https://doi.org/10.1007/s10704-009-9413-9.

Auth, K.L., J. Brouzoulis, and M. Ekh (July 1, 2022). "A Fully Coupled Chemo-Mechanical Cohesive Zone Model for Oxygen Embrittlement of Nickel-Based Superalloys". In: *Journal of the Mechanics and Physics of Solids* 164, p. 104880. ISSN: 00225096. DOI: 10.1016/j.jmps.2022.104880. URL: https://linkinghub.elsevier.com/retrieve/pii/S0022509622000898.

Azab, M., G. Parry, and R. Estevez (Mar. 2020). "An analytical model for DCB/wedge tests based on Timoshenko beam kinematics for accurate determination of cohesive zone lengths". In: *International Journal of Fracture* 222.1-2, pp. 137–153. ISSN: 15732673. DOI: 10.1007/s10704-020-00438-2. URL: https://doi.org/10.1007/s10704-020-00438-2.

Barenblatt, G.I. (Jan. 1, 1962). "The Mathematical Theory of Equilibrium Cracks in Brittle Fracture". In: *Advances in Applied Mechanics*. Ed. by H. L. Dryden et al. Vol. 7. Elsevier, pp. 55–129. DOI: 10.1016/S0065-2156(08)70121-2. URL: https://www.sciencedirect.com/science/article/pii/S0065215608701212.
Bergfeld, B. et al. (Dec. 2021). "Crack propagation speeds in weak snowpack layers". In: *Journal of Glaciology*, pp. 1-14. ISSN: 0022-1430. DOI: 10.1017/JOG.2021.118. URL: https://www.cambridge.org/core/journals/journal-of-glaciology/article/crack-propagation-speeds-in-weak-snowpack-layers/C28905FA3198DF87B0A18AB135FC4766.

Berman, N., G. Cohen, and J. Fineberg (2020). "Dynamics and properties of the cohesive zone in rapid fracture and friction". In: *Physical Review Letters* 125.12, p. 125503. ISSN: 10797114. DOI: 10.1103/PhysRevLett.125.125503. URL: https://doi.org/10.1103/PhysRevLett.125.125503.

Bleyer, J., C. Roux-Langlois, and J.-F. Molinari (Mar. 1, 2017). "Dynamic Crack Propagation with a Variational Phase-Field Model: Limiting Speed, Crack Branching and Velocity-Toughening Mechanisms". In: *International Journal of Fracture* 204.1, pp. 79-100. ISSN: 1573-2673. DOI: 10.1007/s10704-016-0163-1. URL: https://doi.org/10.1007/s10704-016-0163-1.

Bolte, J. and E. Pauwels (July 2021). "Conservative Set Valued Fields, Automatic Differentiation, Stochastic Gradient Methods and Deep Learning". In: *Mathematical Programming* 188.1, pp. 19-51. ISSN: 1436-4646. DOI: 10.1007/s10107-020-01501-5.

Bouchitté, G., G. Buttazzo, and A. Braides (Jan. 1995). "Relaxation results for some free discontinuity problems." In: *Journal für die reine und angewandte Mathematik* 1995.458, pp. 1-18. ISSN: 1435-5345. DOI: 10.1515/crl1.1995.458.1.

Bouissou, S., J.P. Petit, and M. Barquins (Feb. 1998). "Normal Load, Slip Rate and Roughness Influence on the Polymethyl-methacrylate Dynamics of Sliding 1. Stable Sliding to Stick-Slip Transition". In: *Wear* 214.2, pp. 156-164. ISSN: 00431648. DOI: 10.1016/S0043-1648(97)00242-1.

Bybordiani, M. and D. Dias-da-Costa (Apr. 2021). "A consistent finite element approach for dynamic crack propagation with explicit time integration". In: *Computer Methods in Applied Mechanics and Engineering* 376, p. 113652. ISSN: 00457825. DOI: 10.1016/j.cma.2020.113652. URL: https://doi.org/10.1016/j.cma.2020.113652%20https://linkinghub.elsevier.com/retrieve/pii/S0045782520308379.

Camacho, G.T. and M. Ortiz (1996). "Computational modelling of impact damage in brittle materials". In: *International Journal of Solids and Structures* 33.20, pp. 2899-2938. ISSN: 0020-7683. DOI: 10.1016/0020-7683(95)00255-3. URL: http://www.sciencedirect.com/science/article/pii/0020768395002553.

Cаманho, P.P., C.G. Davila, and M.F. De Moura (2003). "Numerical Simulation of Mixed-Mode Progressive Delamination in Composite Materials". In: *Journal of Composite Materials* 37.16, pp. 1415-1438. DOI: 10.1177/002199803034505. URL: https://www.researchgate.net/publication/37650396.

Cazes, F., M. Coret, and A. Combescure (2013). "A two-field modified lagrangian formulation for robust simulations of extrinsic cohesive zone models". In: *Computational Mechanics* 51, pp. 865-884. URL: https://hal.archives-ouvertes.fr/hal-00938515.

Célarié, F. et al. (2003). "Glass Breaks like Metal, but at the Nanometer Scale". In: *Physical Review Letters* 90.7, pp. 075504-1-075504-4. ISSN: 10797114. DOI: 10.1103/PhysRevLett.90.075504.

Chaboche, J.L., F. Feyel, and Y. Monerie (2001). "Interface debonding models: A viscous regularization with a limited rate dependency". In: *International Journal of Solids and Structures* 38.18, pp. 3127-3160. ISSN: 00207683. DOI: 10.1016/S0020-7683(00)00053-6.

Collins-Craft, N.A., F. Bourrier, and V. Acary (Oct. 1, 2022). "On the Formulation and Implementation of Extrinsic Cohesive Zone Models with Contact". In: *Computer Methods in Applied Mechanics and Engineering* 400, p. 115545. ISSN: 0045-7825. DOI: 10.1016/j.cma.2022.115545. URL: https://www.sciencedirect.com/science/article/pii/S0045782522005369.

Cottle, R.W., J.-S. Pang, and R.E. Stone (2009). *The Linear Complementarity Problem*. Ed. by R.E. O'Malley. Second. Philadelphia: Society for Industrial and Applied Mathematics, p. 781. ISBN: 978-0-89871-686-3. DOI: 10.1137/1.9780898719000. URL: http://epubs.siam.org/doi/book/10.1137/1.9780898719000.

Dugdale, D.S. (May 1960). "Yielding of Steel Sheets Containing Slits". In: *Journal of the Mechanics and Physics of Solids* 8.2, pp. 100-104. ISSN: 00225096. DOI: 10.1016/0022-5096(60)90013-2. URL: https://linkinghub.elsevier.com/retrieve/pii/0022509660900132.

Dureisseix, D. et al. (Sept. 18, 2024). "An Explicit Dynamics Framework Suited to Highly Non-Smooth Interface Behaviors". In: *Journal of Theoretical, Computational and Applied Mechanics*. ISSN: 2726-6141. DOI: 10.46298/jtcam.11611. URL: https://jtcam.episciences.org/14330.

Eid, E., R. Seghir, and J. Réthoré (Aug. 18, 2023). "Crack Branching at Low Tip Speeds: Spilling the T". In: *Journal of Theoretical, Computational and Applied Mechanics*. ISSN: 2726-6141. DOI: 10.46298/jtcam.10172. URL: https://jtcam.episciences.org/12218.

Falk, M.L., A. Needleman, and J.R. Rice (2001). "A critical evaluation of dynamic fracture simulations using cohesive surfaces". In: *Journal de Physique IV*. Vol. 11, Pr5-43-Pr5-50. DOI: 10.1051/jp4:2001506.

Felippa, C.A. (2004). *Introduction to Finite Element Methods*. Boulder, Colorado, United States of America: University of Colorado. 794 pp. URL: https://quickfem.com/finite-element-analysis/.

Ferry, R. et al. (Feb. 2025). "Depth Dependence of Coseismic Off-Fault Damage and Its Effects on Rupture Dynamics". In: *Journal of Geophysical Research: Solid Earth* 130.2, e2024JB029787. ISSN: 2169-9313, 2169-9356. DOI: 10.1029/2024JB029787.

Foulk, J.W. (Jan. 2010). "An examination of stability in cohesive zone modeling". In: *Computer Methods in Applied Mechanics and Engineering* 199.9-12, pp. 465-470. ISSN: 0045-7825. DOI: 10.1016/J.CMA.2009.08.025.
Frémond, M. (1988). "Contact with Adhesion". In: *Nonsmooth Mechanics and Applications*. Vienna: Springer Vienna, pp. 177-221. DOI: 10.1007/978-3-7091-2624-0_3. URL: http://link.springer.com/10.1007/978-3-7091-2624-0_3.

- (2002). *Non-Smooth Thermomechanics*. Berlin, Heidelberg: Springer Berlin Heidelberg. ISBN: 978-3-642-08578-9. DOI: 10.1007/978-3-662-04800-9. URL: http://link.springer.com/10.1007/978-3-662-04800-9.

- (2012a). "Contact with Adhesion". In: *Phase Change in Mechanics*. Vol. 13. Lecture Notes of the Unione Matematica Italiana 4. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 151-156. ISBN: 978-3-642-24608-1. DOI: 10.1007/978-3-642-24609-8. URL: http://link.springer.com/10.1007/978-3-642-24609-8.

- (2012b). "Damage of Solids Glued on One Another: Coupling of Volume and Surface Damages". In: *Phase Change in Mechanics*. Lecture Notes of the Unione Matematica Italiana. Berlin, Heidelberg: Springer Berlin Heidelberg, pp. 115-130. ISBN: 978-3-642-24608-1. DOI: 10.1007/978-3-642-24609-8. URL: http://link.springer.com/10.1007/978-3-642-24609-8.

- (2017). *Collisions Engineering: Theory and Applications*. Vol. 6. Springer Series in Solid and Structural Mechanics. Berlin, Heidelberg: Springer. ISBN: 978-3-662-52694-1 978-3-662-52696-5. DOI: 10.1007/978-3-662-52696-5. URL: http://link.springer.com/10.1007/978-3-662-52696-5.

Geuzaine, C. and J.-F. Remacle (Sept. 2009). "Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities". In: *International Journal for Numerical Methods in Engineering* 79.11, pp. 1309-1331. ISSN: 00295981. DOI: 10.102/nme.2579. URL: http://doi.wiley.com/10.1002/nme.2579.

[SW] Geuzaine, C. and J.-F. Remacle, Gmsh 2022. LIC: Modified GPL-2.0. URL: https://gitlab.onelab.info/gmsh/gmsh/-/tree/gmsh_4_11_1?ref_type=tags, SWHID: (swh:1:rev:4608e03c213fb6d19fb9561dfe47d7261c20e872; origin=https://gitlab.onelab.info/gmsh/gmsh; visit=swh:1:sn p:8363430f244b683c0d9a0d1ce531b49edd9b7f60).

Griffith, A.A. (Jan. 1921). "The phenomena of rupture and flow in solids". In: *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 221.582-593, pp. 163-198. ISSN: 0264-3952. DOI: 10.1098/rsta.1921.0006. URL: https://royalsocietypublishing.org/doi/10.1098/rsta.1921.0006.

Guilloteau, E., H. Charrue, and F. Creuzet (1996). "The direct observation of the core region of a propagating fracture crack in glass". In: *Europhysics Letters (EPL)* 34.7, pp. 549-554. ISSN: 0295-5075. DOI: 10.1209/epl/i1996-00493-3.

Gvirtzman, S. and J. Fineberg (Aug. 2, 2021). "Nucleation Fronts Ignite the Interface Rupture That Initiates Frictional Motion". In: *Nature Physics* 2021, pp. 1-6. ISSN: 1745-2481. DOI: 10.1038/s41567-021-01299-9. URL: https://www.nature.com/articles/s41567-021-01299-9.

Halphen, B. and Q.S. Nguyen (Jan. 1, 1975). "Sur les Matériaux Standard Généralisés". In: *Journal de Mécanique* 14, pp. 39-63. URL: https://hal.science/hal-03600755/document.

Hinton, E. and J.S. Campbell (1974). "Local and Global Smoothing of Discontinuous Finite Element Functions Using a Least Squares Method". In: *International Journal for Numerical Methods in Engineering* 8.3, pp. 461-480. ISSN: 1097-0207. DOI: 10.1002/nme.1620080303. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620080303.

Houlsby, G.T. (2019). "Frictional Plasticity in a Convex Analytical Setting". In: *Open Geomechanics* 1.3, pp. 1-10. DOI: 10.5802/ogeo.2. URL: https://opengeomechanics-centre-mersenne.org/item/OGEO_2019_1_A3_0.

Jean, M. (July 1999). "The non-smooth contact dynamics method". In: *Computer Methods in Applied Mechanics and Engineering* 177.3-4, pp. 235-257. ISSN: 00457825. DOI: 10.1016/S0045-7825(98)00383-1. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782598003831.

Jean, M., V. Acary, and Y. Monerie (Dec. 2001). "Non-Smooth Contact Dynamics Approach of Cohesive Materials". In: *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* 359.1789, pp. 2497-2518. ISSN: 1364503X. DOI: 10.1098/rsta.2001.0906.

Jean, M. and J.J. Moreau (1992). "Unilaterality and dry friction in the dynamics of rigid body collections". In: *1st Contact Mechanics International Symposium*. Ed. by A. Curnier. Lausanne, pp. 31-48. URL: https://hal.archives-ouvertes.fr/hal-01863710.

Jiang, S. et al. (Sept. 2021). "Characterisation of fracture evolution of a single cemented brittle grain using in-situ X-ray computed tomography". In: *International Journal of Rock Mechanics and Mining Sciences* 145, p. 104835. ISSN: 1365-1609. DOI: 10.1016/j.ijrmms.2021.104835.

Kubair, D.V. and P.H. Geubelle (2003). "Comparative analysis of extrinsic and intrinsic cohesive models of dynamic fracture". In: *International Journal of Solids and Structures* 40.15, pp. 3853-3868. ISSN: 0020-7683. DOI: https://doi.org/10.1016/S0020-7683(03)00171-9. URL: http://www.sciencedirect.com/science/article/pii/S0020768303001719.

Lebihain, M., T. Roch, and J.-F. Molinari (Nov. 2022). "Quasi-Static Crack Front Deformations in Cohesive Materials". In: *Journal of the Mechanics and Physics of Solids* 168, p. 105025. ISSN: 0022-5096. DOI: 10.1016/j.jmps.2022.105025.

Lemke, C.E. and J.T. Howson, Jr. (June 1964). "Equilibrium Points of Bimatrix Games". In: *Journal of the Society for Industrial and Applied Mathematics* 12.2, pp. 413-423. ISSN: 0368-4245. DOI: 10.1137/0112033. URL: http://epubs.siam.org/doi/10.1137/0112033.
Lorentz, E. (2008). "A mixed interface finite element for cohesive zone models". In: *Computer Methods in Applied Mechanics and Engineering* 198.2, pp. 302-317. ISSN: 0045-7825. DOI: https://doi.org/10.1016/j.cma.2008.08.006.
URL: https://www.sciencedirect.com/science/article/pii/S0045782508002880.

Marfia, S., E. Monaldo, and E. Sacco (June 2022). "Cohesive Fracture Evolution within Virtual Element Method". In: *Engineering Fracture Mechanics* 269, p. 108464. ISSN: 00137944. DOI: 10.1016/j.engfracmech.2022.108464. URL: https://linkinghub.elsevier.com/retrieve/pii/S0013794422002107.

Marigo, J.J. (1981). "Formulation d'une loi d'endommagement d'un materiau élastique". In: *Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre* 292.May, pp. 1309-1312. URL: http://gallica.bnf.fr/ark:/12148/bpt6k36331471h/f337.item.

Moës, N. and T. Belytschko (2002). "Extended Finite Element Method for Cohesive Crack Growth". In: *Engineering Fracture Mechanics* 69.7, pp. 813-833. ISSN: 00137944. DOI: 10.1016/S0013-7944(01)00128-X.

Monerie, Y. and V. Acary (2001). "Formulation dynamique d'un modèle de zone cohésive tridimensionnel couplant endommagement et frottement". In: *Revue Européenne des Éléments Finis* 10.2-4, pp. 489-503. DOI: 10.1080/12506559.2001.11869264. eprint: https://doi.org/10.1080/12506559.2001.11869264. URL: https://doi.org/10.1080/12506559.2001.11869264.

Moreau, J.J. (1970). "Sur les lois de frottement, de plasticité et de viscosité". In: *Comptes rendus de l'Académie des sciences. Série A - Sciences mathématiques* 271, pp. 608-611. URL: https://hal.science/hal-01868140/document.

- (1974). "On Unilateral Constraints, Friction and Plasticity". In: *New Variational Techniques in Mathematical Physics*. Springer Berlin Heidelberg, pp. 171-322. DOI: 10.1007/978-3-642-10960-7_7. URL: https://link.springer.com/chapter/10.1007/978-3-642-10960-7_7.

- (1986). "Une formulation du contact à frottement sec; application au calcul numérique". In: *Comptes rendus de l'Académie des sciences. Série 2, Mécanique, Physique, Chimie, Sciences de l'univers, Sciences de la Terre* 302.13, pp. 799-801. ISSN: 0764-4450. URL: https://gallica.bnf.fr/ark:/12148/bpt6k5675350v/f805.item.

- (July 1999). "Numerical aspects of the sweeping process". In: *Computer Methods in Applied Mechanics and Engineering* 177.3-4, pp. 329-349. ISSN: 0045-7825. DOI: 10.1016/S0045-7825(98)00387-9.

Nguyen, V.P. (2014). "Discontinuous Galerkin/extrinsic cohesive zone modeling: Implementation caveats and applications in computational fracture mechanics". In: *Engineering Fracture Mechanics* 128, pp. 37-68. ISSN: 0013-7944. DOI: 10.1016/j.engfracmech.2014.07.003. URL: http://www.sciencedirect.com/science/article/pii/S0013794414002136.

Nguyen-Thanh, V.M. et al. (Oct. 2018). "A Virtual Element Method for 2D Linear Elastic Fracture Analysis". In: *Computer Methods in Applied Mechanics and Engineering* 340, pp. 366-395. ISSN: 00457825. DOI: 10.1016/j.cma.2018.05.021. URL: https://linkinghub.elsevier.com/retrieve/pii/S0045782518302664.

Nkoumbou Kaptchouang, N.B. et al. (2021). "Cohesive GTN model for ductile fracture simulation". In: *Engineering Fracture Mechanics* 242, p. 107437. ISSN: 0013-7944. DOI: 10.1016/j.engfracmech.2020.107437. URL: https://www.sciencedirect.com/science/article/pii/S0013794420310080.

Okubo, K. et al. (Nov. 2019). "Dynamics, Radiation, and Overall Energy Budget of Earthquake Rupture With Coseismic Off-Fault Damage". In: *Journal of Geophysical Research: Solid Earth* 124.11, pp. 11771-11801. ISSN: 2169-9313. DOI: 10.1029/2019JB017304. eprint: 1901.01771. URL: https://onlinelibrary.wiley.com/doi/10.1029/2019JB017304.

Papoulia, K.D., C.-H. Sam, and S.A. Vavasis (Oct. 2003). "Time continuity in cohesive finite element modeling". In: *International Journal for Numerical Methods in Engineering* 58.5, pp. 679-701. ISSN: 0029-5981. DOI: 10.1002/nme.778. URL: http://doi.wiley.com/10.1002/nme.778.

Papoulia, K.D., S.A. Vavasis, and P. Ganguly (2006). "Spatial Convergence of Crack Nucleation Using a Cohesive Finite-Element Model on a Pinwheel-Based Mesh". In: *International Journal for Numerical Methods in Engineering* 67.1, pp. 1-16. ISSN: 1097-0207. DOI: 10.1002/nme.1598. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1598.

[SW], *ParaView* version 5.11.2, Sept. 21, 2023. Kitware. LIC: CCPL-3.0. URL: https://archive.softwareheritage.org/swh:1:rev:9d280162844f78897446c95a7dac3c977f41d746;origin=https://gitlab.kitware.com/paraview/paraview.git;visit=swh:1:snp:767607ff936f76b248e40a84a1b087994e7d6335, vcs: https://gitlab.kitware.com/paraview/paraview.

Park, K. et al. (2012). "Adaptive Mesh Refinement and Coarsening for Cohesive Zone Modeling of Dynamic Fracture". In: *International Journal for Numerical Methods in Engineering* 92.1, pp. 1-35. ISSN: 1097-0207. DOI: 10.1002/nme.3163. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.3163.

Parrinello, F. (Sept. 2020). "Hybrid equilibrium element with interelement interface for the analysis of delamination and crack propagation problems". In: *International Journal for Numerical Methods in Engineering* August. ISSN: 0029-5981. DOI: 10.1002/nme.6531. URL: https://onlinelibrary.wiley.com/doi/10.1002/nme.6531.

Parrinello, F. and G. Borino (2020). "Cohesive-frictional interface in an equilibrium based finite element formulation". In: *Proceedings of XXIV AIMETA Conference 2019*. Ed. by A. Carcaterra, A. Paolone, and G. Graziani. Lecture Notes in Mechanical Engineering September. Cham: Springer International Publishing, pp. 419-426. ISBN: 978-3-030-41056-8. DOI: 10.1007/978-3-030-41057-5. URL: http://link.springer.com/10.1007/978-3-030-41057-5.
Perales, F. et al. (2010). "A NonSmooth Contact Dynamics-based multi-domain solver". In: *European Journal of Computational Mechanics* 19.4, pp. 389-417. DOI: 10.3166/ejcm.19.389-417. eprint: https://doi.org/10.3166/ejcm.19.389-417. URL: https://doi.org/10.3166/ejcm.19.389-417.

Rabczuk, T. and T. Belytschko (2004). "Cracking Particles: A Simplified Meshfree Method for Arbitrary Evolving Cracks". In: *International Journal for Numerical Methods in Engineering* 61.13, pp. 2316-2343. ISSN: 1097-0207. DOI: 10.1002/nme.1151. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1151.

Ramulu, M. and A.S. Kobayashi (Mar. 1, 1985). "Mechanics of Crack Curving and Branching – a Dynamic Fracture Analysis". In: *International Journal of Fracture* 27.3, pp. 187-201. ISSN: 1573-2673. DOI: 10.1007/BF00017967. URL: https://doi.org/10.1007/BF00017967.

Raous, M., L. Cangémi, and M. Cocou (1999). "A consistent model coupling adhesion, friction, and unilateral contact". In: *Computer Methods in Applied Mechanics and Engineering*. DOI: 10.1016/S0045-7825(98)00389-2. URL: https://hal.archives-ouvertes.fr/hal-03178187.

Réthoré, J. and R. Estevez (2013). "Identification of a cohesive zone model from digital images at the micron-scale". In: *Journal of the Mechanics and Physics of Solids* 61.6, pp. 1407-1420. ISSN: 00225096. DOI: 10.1016/j.jmps.2013.01.011.

[SW] Richart, N. et al., *Akantu 5.0.7-post1* version 5.0.7-post1, Sept. 20, 2024. École Polytechnique Fédérale de Lausanne. LIC: LGPL-3.0. URL: https://archive.softwareheritage.org/swh:1:rev:23af7633761c4b0d93dd34d71080eee8d8b385d4;origin=https://gitlab.com/akantu/akantu;visit=swh:1:snp:804a1588aa902e61a5dacdaca3c3dc88b35b66e1, vcs: https://gitlab.com/akantu/akantu, SWHID: <swh:1:rev:23af7633761c4b0d93dd34d71080eee8d8b385d4;origin=https://gitlab.com/akantu/akantu;visit=swh:1:snp:804a1588aa902e61a5dacdaca3c3dc88b35b66e1>.

Richart, N. et al. (Feb. 22, 2024b). "Akantu: An HPC Finite-Element Library for Contact and Dynamic Fracture Simulations". In: *Journal of Open Source Software* 9.94, p. 5253. ISSN: 2475-9066. DOI: 10.21105/joss.05253. URL: https://joss.theoj.org/papers/10.21105/joss.05253.

Roch, T., M. Lebihain, and J.-F. Molinari (Aug. 2023). "Dynamic Crack-Front Deformations in Cohesive Materials". In: *Physical Review Letters* 131.9, p. 096101. DOI: 10.1103/PhysRevLett.131.096101.

Rodella, A. et al. (Mar. 2025). *Sharp-Interface Cohesive Fracture Models with Relaxed Bulk Energies: Numerical Investigations*. URL: https://hal.sorbonne-universite.fr/hal-05008658/.

Sam, C.-H., K.D. Papoulia, and S.A. Vavasis (Sept. 2005). "Obtaining initially rigid cohesive finite element models that are temporally convergent". In: *Engineering Fracture Mechanics* 72.14, pp. 2247-2267. ISSN: 00137944. DOI: 10.1016/j.engfracmech.2004.12.008. URL: https://linkinghub.elsevier.com/retrieve/pii/S0013794405000895.

Samimi, M., J.A.W. van Dommelen, and M.G.D. Geers (Dec. 2011). "A three-dimensional self-adaptive cohesive zone model for interfacial delamination". In: *Computer Methods in Applied Mechanics and Engineering* 200.49-52, pp. 3540-3553. ISSN: 0045-7825. DOI: 10.1016/J.CMA.2011.08.021.

[SW] Schlömer, N., *meshio: Tools for mesh files* 2024. LIC: MIT. DOI: 10.5281/zenodo.1288334, URL: https://git hub.com/nschloe/meshio/tree/v5.3.5, SWHID: <swh:1:dir:6e1cd14ba43f09c1e729f41f4b6e2ad5fb0d6c24;origin=github.com/nschloe/meshio;visit=swh:1:snp:9d55efa5b1f7ff0238d8557c5aee1673b55e3d29;anchor=swh:1:rev:b2ee99842e119901349fdeee06b5bf61e01f450a>.

Seagraves, A. and R. Radovitzky (2009). "Advances in Cohesive Zone Modeling of Dynamic Fracture". In: *Dynamic Failure of Materials and Structures*. Boston, MA: Springer US, pp. 349-405. ISBN: 978-1-4419-0445-4. DOI: 10.1007/978-1-4419-0446-1_12. URL: https://link.springer.com/chapter/10.1007/978-1-4419-0446-1_12.

Snozzi, L. and J.-F. Molinari (Feb. 2013). "A Cohesive Element Model for Mixed Mode Loading with Frictional Contact Capability". In: *International Journal for Numerical Methods in Engineering* 93.5, pp. 510-526. ISSN: 00295981. DOI: 10.1002/nme.4398.

Song, J.-H., H. Wang, and T. Belytschko (July 1, 2008). "A Comparative Study on Finite Element Methods for Dynamic Fracture". In: *Computational Mechanics* 42.2, pp. 239-250. ISSN: 1432-0924. DOI: 10.1007/s00466-007-0210-x. URL: https://doi.org/10.1007/s00466-007-0210-x.

Stewart, D.E. and J.C. Trinkle (Aug. 15, 1996). "An Implicit Time-Stepping Scheme for Rigid Body Dynamics with Inelastic Collisions and Coulomb Friction". In: *International Journal for Numerical Methods in Engineering* 39.15, pp. 2673-2691. ISSN: 0029-5981. DOI: 10.1002/(SICI)1097-0207(19960815)39:15<2673::AID-NME972>3.0.CO;2-I. URL: https://onlinelibrary.wiley.com/doi/10.1002/(SICI)1097-0207(19960815)39:15%3C2673::AID-NME972%3E3.0.CO;2-I.

Talon, C. and A. Curnier (2003). "A model of adhesion coupled to contact and friction". In: *European Journal of Mechanics - A/Solids* 22.4, pp. 545-565. ISSN: 0997-7538. DOI: https://doi.org/10.1016/S0997-7538(03)00046-9. URL: https://www.sciencedirect.com/science/article/pii/S0997753803000469.

Vargas, R. et al. (Dec. 2020). "On the identification of cohesive zone model for curved crack in mortar". In: *Strain* 56.6. ISSN: 0039-2103. DOI: 10.1111/str.12364. URL: https://doi.org/10.1111/str.12364%20https://onlinelibrary.wiley.com/doi/10.1111/str.12364.

Venzal, V. et al. (Aug. 2020). "Frictional cohesive zone model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional behaviors". In: *International Journal of Solids and Structures* 198, pp. 17-30. ISSN: 00207683. DOI:
10.1016/j.ijsolstr.2020.04.023. URL: https://linkinghub.elsevier.com/retrieve/pii/S0020768320301402.

Versino, D. et al. (2015). "A thermodynamically consistent discontinuous Galerkin formulation for interface separation". In: *Composite Structures* 133, pp. 595-606. ISSN: 0263-8223. DOI: 10.1016/j.compstruct.2015.07.080. URL: http://www.sciencedirect.com/science/article/pii/S0263822315006261.

Virtanen, P. et al. (Mar. 2020). "SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python". In: *Nature Methods* 17.3, pp. 261-272. ISSN: 1548-7105. DOI: 10.1038/s41592-019-0686-2. URL: https://www.nature.com/articles/s41592-019-0686-2.

Xu, X.-P. and A. Needleman (Sept. 1994). "Numerical simulations of fast crack growth in brittle solids". In: *Journal of the Mechanics and Physics of Solids* 42.9, pp. 1397-1434. ISSN: 00225096. DOI: 10.1016/0022-5096(94)90003-5. URL: https://linkinghub.elsevier.com/retrieve/pii/0022509694900035.

Zhou, F. and J.-F. Molinari (Jan. 2004). "Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency". In: *International Journal for Numerical Methods in Engineering* 59.1, pp. 1-24. ISSN: 0029-5981. DOI: 10.1002/nme.857. URL: http://doi.wiley.com/10.1002/nme.857.