# Fragmentation analysis of a bar with the Lip-field approach

Nicolas Moës, Benoît Lé, Andrew Stershic

► To cite this version:

Nicolas Moës, Benoît Lé, Andrew Stershic. Fragmentation analysis of a bar with the Lip-field approach. Mechanics of Materials, 2022, 172, pp.104365. 10.1016/j.mechmat.2022.104365. hal-04585728

HAL Id: hal-04585728

https://hal.science/hal-04585728

Submitted on 23 May 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# FRAGMENTATION ANALYSIS OF A BAR WITH THE LIP-FIELD APPROACH

A PREPRINT

Nicolas Moës

École Centrale de Nantes
GeM Institute, UMR CNRS 6183
Institut Universitaire de France (IUF)
1 rue de la Noë, 44321 Nantes, France
nicolas.moes@ec-nantes.fr

Benoît Lé *

École Centrale de Nantes
GeM Institute, UMR CNRS 6183
1 rue de la Noë, 44321 Nantes, France
benoit.le@ec-nantes.fr

Andrew Stershic

Sandia National Laboratories, California
7011 East Avenue, Livermore, CA 94550, USA
ajsters@sandia.gov

January 2022

## ABSTRACT

The Lip-field approach was introduced in [1] as a new way to regularize softening material models. It was tested in 1D quasistatic in [1] and 2D quasistatic in [2]: this paper extends it to 1D dynamics, on the challenging problem of dynamic fragmentation. The Lip-field approach formulates the mechanical problem to be solved as an optimization problem, where the incremental potential to be minimized is the non-regularized one. Spurious localization is prevented by imposing a Lipschitz constraint on the damage field. The displacement and damage field at each time step are obtained by a staggered algorithm, that is the displacement field is computed for a fixed damage field, then the damage field is computed for a fixed displacement field. Indeed, these two problems are convex, which is not the case of the global problem where the displacement and damage fields are sought at the same time. The incremental potential is obtained by equivalence with a cohesive zone model, which makes material parameters calibration simple. A non-regularized local damage equivalent to a cohesive zone model is also proposed. It is used as a reference for the Lip-field approach, without the need to implement displacement jumps. These approaches are applied to the brittle fragmentation of a 1D bar with randomly perturbed material properties to accelerate spatial convergence. Both explicit and implicit dynamic implementations are compared. Favorable comparison to several analytical, numerical and experimental references serves to validate the modeling approach.

**Keywords** Fragmentation, Damage, Fracture, Dynamic, Lipschitz, Lip-field

## 1 Introduction

Fragmentation is a challenging application for failure modeling approaches, as failure is pervasive through the domain. It is indeed a multiscale problem, balancing the localization processes that lead to fracture surfaces with a significant interactivity between the multiple failure points. Accordingly, these competing physics lead to a material response varies which significantly with loading rate; with quasistatic and low-velocity loadings, a small number of dominant cracks is sufficient to relieve the structure; in contrast, strongly dynamic loadings lead to rate-dependent multiple fragmentation [3].

Much study on fragmentation has been motivated by military application, with foundational studies performed by Mott [4], supplemented by a number of experimental studies (e.g. [5]), analytical models (e.g. [6, 7]), and numerical
FRAGMENTATION ANALYSIS OF A BAR WITH THE LIP-FIELD
APPROACH

A PREPRINT

Nicolas Moës

École Centrale de Nantes
GeM Institute, UMR CNRS 6183
Institut Universitaire de France (IUF)
1 rue de la Noë, 44321 Nantes, France
nicolas.moes@ec-nantes.fr

Benoît Lé *

École Centrale de Nantes
GeM Institute, UMR CNRS 6183
1 rue de la Noë, 44321 Nantes, France
benoit.le@ec-nantes.fr

Andrew Stershic

Sandia National Laboratories, California
7011 East Avenue, Livermore, CA 94550, USA
ajsters@sandia.gov

January 2022

ABSTRACT

The Lip-field approach was introduced in [1] as a new way to regularize softening material models. It was tested in 1D quasistatic in [1] and 2D quasistatic in [2]: this paper extends it to 1D dynamics, on the challenging problem of dynamic fragmentation. The Lip-field approach formulates the mechanical problem to be solved as an optimization problem, where the incremental potential to be minimized is the non-regularized one. Spurious localization is prevented by imposing a Lipschitz constraint on the damage field. The displacement and damage field at each time step are obtained by a staggered algorithm, that is the displacement field is computed for a fixed damage field, then the damage field is computed for a fixed displacement field. Indeed, these two problems are convex, which is not the case of the global problem where the displacement and damage fields are sought at the same time. The incremental potential is obtained by equivalence with a cohesive zone model, which makes material parameters calibration simple. A non-regularized local damage equivalent to a cohesive zone model is also proposed. It is used as a reference for the Lip-field approach, without the need to implement displacement jumps. These approaches are applied to the brittle fragmentation of a 1D bar with randomly perturbed material properties to accelerate spatial convergence. Both explicit and implicit dynamic implementations are compared. Favorable comparison to several analytical, numerical and experimental references serves to validate the modeling approach.

**Keywords** Fragmentation, Damage, Fracture, Dynamic, Lipschitz, Lip-field

# 1 Introduction

Fragmentation is a challenging application for failure modeling approaches, as failure is pervasive through the domain. It is indeed a multiscale problem, balancing the localization processes that lead to fracture surfaces with a significant interactivity between the multiple failure points. Accordingly, these competing physics lead to a material response varies which significantly with loading rate; with quasistatic and low-velocity loadings, a small number of dominant cracks is sufficient to relieve the structure; in contrast, strongly dynamic loadings lead to rate-dependent multiple fragmentation [3].

Much study on fragmentation has been motivated by military application, with foundational studies performed by Mott [4], supplemented by a number of experimental studies (e.g. [5]), analytical models (e.g. [6, 7]), and numerical
models (e.g. [8, 9, 10, 11]). The original motivation persists (e.g. [12, 13]), and the topic has additional relevance in such diverse fields such as ice-sheet morphology (e.g. [14, 15]), geophysics (e.g. [16, 17]), and nuclear reactor fuels (e.g. [18]).

Fracture of softening materials is physically characterized by a process zone ahead of the cracks tips where micro-
cracking nucleate and coalesce to give macro-cracks. Taking into account the characteristic length of this process zone
is a key point in numerical modelling of softening materials.

A first and very popular model to consider this characteristic length is the cohesive zone model (CZM) [19, 20, 21, 22]. It consists in modeling the macro-cracks by displacement jumps, with cohesive forces acting on the crack lips over a certain distance before the crack tips. These cohesive zones are usually introduced at the boundaries between finite elements, which make them very dependent on the mesh orientation. Despite this limitation, this cohesive element approach has been advanced for fragmentation studies, especially for quasi-one-dimensional problems [8, 9, 10, 11, 23]. Some approaches like the eXtended finite element method (X-FEM) [24] alleviate the mesh dependency by allowing the introduction of cohesive cracks inside the finite elements, however modelling complex crack patterns like branching or coalescence remains challenging from a geometry perspective.

Other means of discrete fracture representation that have been used to model fragmentation include various meshfree and particle methods, such as: smooth particle hydrodynamics [25, 26], material point method [27], discrete element method [28], reproducing kernel particle method (RKPM) [29], and peridynamics [30]. There is an impressive diversity in these approaches, and an exhaustive review is beyond the scope of this work. However, it bears mention that several meshfree methods have been proposed that integrate aspects of cohesive modeling in order to represent Griffith failure. Exemplars include the virtual internal bond method [31], RKPM with cohesive tractions [32, 33], cracking particles method [34], cohesive molecular dynamics methods [35], bond-based peridynamics with energetic criteria [36], multiscale cohesive atomistic methods [37], and the optimal transport method coupled with eigenerosion [27].

An alternative class of numerical failure modeling is continuum damage models. These models represent the influence of micro-cracks by means of some internal damage variable [38, 39]. With this approach complex crack morphology (e.g. branching, networks, fragmentation) is handled natively, even if the exact position of macro-cracks is not available anymore. The failure of the material is represented by a softening constitutive model; however, local damage models suffer from the well-known problem of spurious mesh dependency [40], that is a dissipated energy which tends to zero when refining finite element meshes.

These models need a length scale $l_c$ to recover well-posedness. This problem has been widely studied, see for instance higher-order gradient models [41, 42] or regularization of internal variables [43, 44, 45, 46, 47]. Alternatively, a local damage model can be regularized by introducing a characteristic time instead of a characteristic length [48, 49, 50, 51]. Delay damage models and nonlocal integral models are compared in [52]; it seems that even if both models can regularize local damage models, their validity domain depends on the strain rates. A recent assessment of the delay damage model in limit cases may also be found in [53].

Another approach of recent interest is the phase-field fracture model [54, 55], which derives from a smoothing of fracture mechanics. Based on the variational approach to fracture [56], it uses a smoothed representation of the macro-cracks to get a process zone with a finite thickness and can deal with complex crack topologies. It was initially developed for brittle fracture, but more recent works have extended it to cohesive fracture [57, 58], enabling application to ductile fracture [59]. The phase-field approach has also been applied to dynamic brittle fracture [60, 61, 62, 63].

The Thick Level Set (TLS) approach introduced by [64, 65] defines the damage variable as a user-defined function of a level-set function corresponding to the distance to the damage front (boundary between the damaged and undamaged material). In this framework, the iso-$l_c$ of this level-set function gives the crack position, which can be enriched to introduce displacement jumps [66]. Complex crack patterns can be captured using specific element cutting algorithms [67]. The TLS has been applied to dynamic fracture [68] and fragmentation [69].

A new regularization method, called the Lip-field approach, was introduced in [1]. It is based on formulating the mechanical problem to be solved as an optimization problem, where the unknowns are the displacement field and the damage variable (with eventually other internal variables like plastic deformation for instance). The expression of the potential to be minimized is exactly the same as the non-regularized problem, however a Lipschitz constraint is applied on the damage variable to keep its gradient bounded. The idea is similar to the graded damage model from [70, 71], where some constraints are applied to the damage variable and its gradient through some Lagrange multipliers. However, the implementation is different; the Lip-field approach does not need the gradient of the damage variable, which can be kept at the integration points. This makes the implementation of Lip-field in a standard finite element code easier. Also, [1] demonstrated the existence of bounds which can be used to decrease the size of the domain where the damage variable has to be computed, reducing computing times. Note that the idea of keeping the gradient of some problem variable that Lip-field and graded damage have in common is not new, see for instance slope
models (e.g. [8, 9, 10, 11]). The original motivation persists (e.g. [12, 13]), and the topic has additional relevance in such diverse fields such as ice-sheet morphology (e.g. [14, 15]), geophysics (e.g. [16, 17]), and nuclear reactor fuels (e.g. [18]).

Fracture of softening materials is physically characterized by a process zone ahead of the cracks tips where micro-
cracking nucleate and coalesce to give macro-cracks. Taking into account the characteristic length of this process zone
is a key point in numerical modelling of softening materials.

A first and very popular model to consider this characteristic length is the cohesive zone model (CZM) [19, 20, 21, 22]. It consists in modeling the macro-cracks by displacement jumps, with cohesive forces acting on the crack lips over a certain distance before the crack tips. These cohesive zones are usually introduced at the boundaries between finite elements, which make them very dependent on the mesh orientation. Despite this limitation, this cohesive element approach has been advanced for fragmentation studies, especially for quasi-one-dimensional problems [8, 9, 10, 11, 23]. Some approaches like the eXtended finite element method (X-FEM) [24] alleviate the mesh dependency by allowing the introduction of cohesive cracks inside the finite elements, however modelling complex crack patterns like branching or coalescence remains challenging from a geometry perspective.

Other means of discrete fracture representation that have been used to model fragmentation include various meshfree and particle methods, such as: smooth particle hydrodynamics [25, 26], material point method [27], discrete element method [28], reproducing kernel particle method (RKPM) [29], and peridynamics [30]. There is an impressive diversity in these approaches, and an exhaustive review is beyond the scope of this work. However, it bears mention that several meshfree methods have been proposed that integrate aspects of cohesive modeling in order to represent Griffith failure. Exemplars include the virtual internal bond method [31], RKPM with cohesive tractions [32, 33], cracking particles method [34], cohesive molecular dynamics methods [35], bond-based peridynamics with energetic criteria [36], multiscale cohesive atomistic methods [37], and the optimal transport method coupled with eigenerosion [27].

An alternative class of numerical failure modeling is continuum damage models. These models represent the influence of micro-cracks by means of some internal damage variable [38, 39]. With this approach complex crack morphology (e.g. branching, networks, fragmentation) is handled natively, even if the exact position of macro-cracks is not available anymore. The failure of the material is represented by a softening constitutive model; however, local damage models suffer from the well-known problem of spurious mesh dependency [40], that is a dissipated energy which tends to zero when refining finite element meshes.

These models need a length scale $l_c$ to recover well-posedness. This problem has been widely studied, see for instance higher-order gradient models [41, 42] or regularization of internal variables [43, 44, 45, 46, 47]. Alternatively, a local damage model can be regularized by introducing a characteristic time instead of a characteristic length [48, 49, 50, 51]. Delay damage models and nonlocal integral models are compared in [52]; it seems that even if both models can regularize local damage models, their validity domain depends on the strain rates. A recent assessment of the delay damage model in limit cases may also be found in [53].

Another approach of recent interest is the phase-field fracture model [54, 55], which derives from a smoothing of fracture mechanics. Based on the variational approach to fracture [56], it uses a smoothed representation of the macro-cracks to get a process zone with a finite thickness and can deal with complex crack topologies. It was initially developed for brittle fracture, but more recent works have extended it to cohesive fracture [57, 58], enabling application to ductile fracture [59]. The phase-field approach has also been applied to dynamic brittle fracture [60, 61, 62, 63].

The Thick Level Set (TLS) approach introduced by [64, 65] defines the damage variable as a user-defined function of a level-set function corresponding to the distance to the damage front (boundary between the damaged and undamaged material). In this framework, the iso-$l_c$ of this level-set function gives the crack position, which can be enriched to introduce displacement jumps [66]. Complex crack patterns can be captured using specific element cutting algorithms [67]. The TLS has been applied to dynamic fracture [68] and fragmentation [69].

A new regularization method, called the Lip-field approach, was introduced in [1]. It is based on formulating the mechanical problem to be solved as an optimization problem, where the unknowns are the displacement field and the damage variable (with eventually other internal variables like plastic deformation for instance). The expression of the potential to be minimized is exactly the same as the non-regularized problem, however a Lipschitz constraint is applied on the damage variable to keep its gradient bounded. The idea is similar to the graded damage model from [70, 71], where some constraints are applied to the damage variable and its gradient through some Lagrange multipliers. However, the implementation is different; the Lip-field approach does not need the gradient of the damage variable, which can be kept at the integration points. This makes the implementation of Lip-field in a standard finite element code easier. Also, [1] demonstrated the existence of bounds which can be used to decrease the size of the domain where the damage variable has to be computed, reducing computing times. Note that the idea of keeping the gradient of some problem variable that Lip-field and graded damage have in common is not new, see for instance slope
constrained topology optimization [72]. This approach was tested in 1D and 2D on quasistatic simulations [1, 2], with damage and softening plasticity (1D case) models. Compared to the TLS approach, the Lip-field approach does not require the level-set technology and yields a convex optimization on the damage variable. Also, as for the phase-field approach the Lip-field solution proceeds with a staggered scheme in which each step is rather straightforward to implement in the finite element framework (a Python Lip-field implementation is provided in [2]).

In this paper, the Lip-field is extended to 1D dynamics on the aforementioned fragmentation problem. The paper is organized as follows. Section 2 describes the mechanical formulation of the 1D fragmentation problem. Section 3 focuses on the Lip-field regularization. Section 4 explains how to get a Lip-field model equivalent to a linear CZM. It also gives a simplified implementation of a CZM which can be used as comparison. Details about the numerical resolution are given in Section 5. Some numerical results are presented in Section 6. These results are discussed in Section 7. Section 8 concludes the paper.

## 2 The 1D fragmentation problem

We consider the deformation with respect to the time variable $t$ of a 1D bar $\Omega = [0; L]$, represented by its elongation $u(x, t)$ for $(x, t) \in \Omega \times [0; T_{\max}]$. In the framework of small-strains assumption, the uniaxial strain $\epsilon$ writes

$$ \epsilon = u_{,x} \qquad (1) $$

We look for the bar displacement $u \in \mathcal{U}$ such that

$$ \int_{0}^{L} \sigma \epsilon(u^*) \, dx = \int_{0}^{L} \rho \ddot{u} u^* \, dx, \quad \forall u^* \in \mathcal{U}^* \qquad (2) $$

where $\sigma$ is the uniaxial stress and $\rho$ the density. $\mathcal{U}(t)$ is the space of the kinematically admissible solutions

$$ \mathcal{U}(t) = \{u \in H^1(\Omega) : u(0, t) = 0, u(L, t) = \dot{\epsilon}_0 Lt\} \qquad (3) $$

and $\mathcal{U}^*$ the test functions space

$$ \mathcal{U}^* = \{u \in H^1(\Omega) : u^*(0) = 0, u^*(L) = 0\} \qquad (4) $$

$\dot{\epsilon}_0$ is the initial, uniform strain rate that is imposed on the bar as initial condition:

$$ u(x, 0) = 0, \quad \forall x \in \Omega \qquad (5) $$

$$ \dot{u}(x, 0) = \dot{\epsilon}_0 x, \quad \forall x \in \Omega \qquad (6) $$

Imposing experimentally such an initial velocity condition in $\Omega$ can be done, among others, by pulsed magnetic fields, see for instance the work of [73, 74] who used such an experimental setup to study the fragmentation of metal rings. The bar is made of a brittle material characterized by a free energy

$$ \varphi(\epsilon, d) = \frac{1}{2} g(d) E \epsilon^2 \qquad (7) $$

where $d$ is a damage variable and $E$ the Young modulus. $g(d)$ is a decreasing function such that $g(0) = 1$ (undamaged material) and $g(1) = 0$ (fully damaged material). The state equations are

$$ \sigma = \frac{\partial \varphi}{\partial \epsilon} = Eg(d)\epsilon, \quad Y = -\frac{\partial \varphi}{\partial d} = -\frac{1}{2}g'(d)E\epsilon^2 \qquad (8) $$

where the prime indicates a derivative with respect to $d$. The evolution laws are chosen as

$$ \dot{d} \ge 0, \quad Y - Y_c H(d) \le 0, \quad (Y - Y_c H(d)) \dot{d} = 0 \qquad (9) $$

where $Y_c$ is the critical energy release rate and $H$ is an increasing function such that $H(0) = 1$, called the “softening function” in what follows. Note that the energy release rate is symmetrical in tension/compression. Asymmetry did not seem necessary to obtain good results as it will be shown in Section 6, but could be easily implemented if needed, see [2].

Solving the fragmentation problem described by the above equations requires to determine, at each time $t_n$, the displacement field $u_n$ and the damage field $d_n$. In what follows, we will assume that $(u_n, d_n)$ are known, and we seek to compute $(u_{n+1}, d_{n+1})$ at time $t_{n+1} = t_n + \Delta t$ (the “$n+1$” indexes are hereafter dropped for simplicity). It can be formulated as an optimization problem:

$$ (u, d) = \arg \min_{\tilde{u} \in U, \tilde{d} \in D} F(\tilde{u}, \tilde{d}) \qquad (10a) $$

$$ F(u, d) = \int_{0}^{L} \frac{1}{2} \rho \frac{(u - u_{\text{p}})^2}{\beta \Delta t^2} + \varphi(\epsilon(u), d) + Y_{c}h(d) \, dx \qquad (10b) $$
constrained topology optimization [72]. This approach was tested in 1D and 2D on quasistatic simulations [1, 2], with damage and softening plasticity (1D case) models. Compared to the TLS approach, the Lip-field approach does not require the level-set technology and yields a convex optimization on the damage variable. Also, as for the phase-field approach the Lip-field solution proceeds with a staggered scheme in which each step is rather straightforward to implement in the finite element framework (a Python Lip-field implementation is provided in [2]).

In this paper, the Lip-field is extended to 1D dynamics on the aforementioned fragmentation problem. The paper is organized as follows. Section 2 describes the mechanical formulation of the 1D fragmentation problem. Section 3 focuses on the Lip-field regularization. Section 4 explains how to get a Lip-field model equivalent to a linear CZM. It also gives a simplified implementation of a CZM which can be used as comparison. Details about the numerical resolution are given in Section 5. Some numerical results are presented in Section 6. These results are discussed in Section 7. Section 8 concludes the paper.

## 2 The 1D fragmentation problem

We consider the deformation with respect to the time variable $t$ of a 1D bar $\Omega = [0; L]$, represented by its elongation $u(x, t)$ for $(x, t) \in \Omega \times [0; T_{\max}]$. In the framework of small-strains assumption, the uniaxial strain $\epsilon$ writes

$$ \epsilon = u_{,x} \qquad (1) $$

We look for the bar displacement $u \in \mathcal{U}$ such that

$$ \int_{0}^{L} \sigma \epsilon(u^*) \, dx = \int_{0}^{L} \rho \ddot{u} u^* \, dx, \quad \forall u^* \in \mathcal{U}^* \qquad (2) $$

where $\sigma$ is the uniaxial stress and $\rho$ the density. $\mathcal{U}(t)$ is the space of the kinematically admissible solutions

$$ \mathcal{U}(t) = \{u \in H^1(\Omega) : u(0, t) = 0, u(L, t) = \dot{\epsilon}_0 Lt\} \qquad (3) $$

and $\mathcal{U}^*$ the test functions space

$$ \mathcal{U}^* = \{u \in H^1(\Omega) : u^*(0) = 0, u^*(L) = 0\} \qquad (4) $$

$\dot{\epsilon}_0$ is the initial, uniform strain rate that is imposed on the bar as initial condition:

$$ u(x, 0) = 0, \quad \forall x \in \Omega \qquad (5) $$

$$ \dot{u}(x, 0) = \dot{\epsilon}_0 x, \quad \forall x \in \Omega \qquad (6) $$

Imposing experimentally such an initial velocity condition in $\Omega$ can be done, among others, by pulsed magnetic fields, see for instance the work of [73, 74] who used such an experimental setup to study the fragmentation of metal rings. The bar is made of a brittle material characterized by a free energy

$$ \varphi(\epsilon, d) = \frac{1}{2} g(d) E \epsilon^2 \qquad (7) $$

where $d$ is a damage variable and $E$ the Young modulus. $g(d)$ is a decreasing function such that $g(0) = 1$ (undamaged material) and $g(1) = 0$ (fully damaged material). The state equations are

$$ \sigma = \frac{\partial \varphi}{\partial \epsilon} = Eg(d)\epsilon, \quad Y = -\frac{\partial \varphi}{\partial d} = -\frac{1}{2}g'(d)E\epsilon^2 \qquad (8) $$

where the prime indicates a derivative with respect to $d$. The evolution laws are chosen as

$$ \dot{d} \ge 0, \quad Y - Y_c H(d) \le 0, \quad (Y - Y_c H(d)) \dot{d} = 0 \qquad (9) $$

where $Y_c$ is the critical energy release rate and $H$ is an increasing function such that $H(0) = 1$, called the “softening function” in what follows. Note that the energy release rate is symmetrical in tension/compression. Asymmetry did not seem necessary to obtain good results as it will be shown in Section 6, but could be easily implemented if needed, see [2].

Solving the fragmentation problem described by the above equations requires to determine, at each time $t_n$, the displacement field $u_n$ and the damage field $d_n$. In what follows, we will assume that $(u_n, d_n)$ are known, and we seek to compute $(u_{n+1}, d_{n+1})$ at time $t_{n+1} = t_n + \Delta t$ (the “$n+1$” indexes are hereafter dropped for simplicity). It can be formulated as an optimization problem:

$$ (u, d) = \arg \min_{\tilde{u} \in \mathcal{U}, \tilde{d} \in \mathcal{D}} F(\tilde{u}, \tilde{d}) \qquad (10a) $$

$$ F(u, d) = \int_{0}^{L} \frac{1}{2} \rho \frac{(u - u_{\text{p}})^2}{\beta \Delta t^2} + \varphi(\epsilon(u), d) + Y_{c}h(d) \, dx \qquad (10b) $$
where

$$h(d) = \int_{0}^{d} H(\tilde{d}) \, d\tilde{d} \qquad (11)$$

and

$$U = \{u \in H^1(\Omega) : u(0) = 0, u(L) = \dot{\epsilon}_0 L t_{n+1}\} \qquad (12)$$

$$D = \{d \in L^{\infty}(\Omega) : d_n \le d \le 1\} \qquad (13)$$

The way the Dirichlet boundary conditions that appear in (12) are applied is detailed in Section 5. Note that the definition of D automatically ensures the irreversibility of the evolution of d. The predictor $u_p$ is defined by

$$u_p = u_n + \Delta t v_n + \frac{\Delta t^2}{2} (1 - 2\beta) a_n \qquad (14)$$

$$v_p = v_n + (1 - \gamma) \Delta t a_n \qquad (15)$$

with $v$ the velocity and $a$ the acceleration. The numerical parameters $\beta$ and $\gamma$ determine the type of numerical integration (explicit or implicit). Then, $a$ is computed (see detail in Section 5) and gives the final values for $u$ and $v$

$$u = u_p + \beta \Delta t^2 a \qquad (16)$$

$$v = v_p + \gamma \Delta t a \qquad (17)$$

The choice $\beta = 0$ corresponds to the explicit dynamics (central difference method). The new displacement is found by $u = u_p$ and thus will not depend on the coming value of damage. The case $\beta = 1/4$, $\gamma = 1/2$ corresponds to the implicit Newmark scheme (undamped trapezoidal rule). Note that for the explicit scheme, it is common to consider a lumped mass matrix, whereas for the implicit scheme, an exact mass matrix integration is performed. Also, taking the variation of (10b) with respect to $u$ gives the time discretization of the weak form of the equilibrium (2), where the acceleration is obtained by equation (16) :

$$\int_{0}^{L} \sigma \epsilon(u^*) \, dx = \int_{0}^{L} \rho \frac{(u - u_p)}{\beta \Delta t^2} u^* \, dx, \quad \forall u^* \in \mathcal{U}^* \qquad (18)$$

while the variation of (10b) with respect to $d$ and the constrain $d \in D$ gives equations (9).

Solving (10) in space $D$ is an ill-posed problem, which leads to solutions exhibiting spurious mesh dependency. Therefore, the Lip-field approach [1] introduces an alternative problem

$$(u, d) = \arg \min_{\tilde{u} \in U, \tilde{d} \in D \cup L} F(\tilde{u}, \tilde{d}) \qquad (19)$$

where

$$L = \{d \in L^{\infty}([0, L]) : |d(x) - d(y)| \le |x - y|/\ell, \quad \forall x, y \in [0, L]\} \qquad (20)$$

The only difference between problems (10) and (19) is the space where the solution $d$ is computed: by enforcing that $d$ must belong to $L$, we prevent the gradient of $d$ to go to infinity, thus avoiding spurious localization. In practice, this problem is not convex with respect to $(u, d)$, therefore, it is solved with a staggered iterative process²:

$$u^{k+1} = \arg \min_{\tilde{u} \in U} F(\tilde{u}, d^k) \qquad (21a)$$

$$d^{k+1} = \arg \min_{\tilde{d} \in D \cup L} F(u^{k+1}, \tilde{d}) \qquad (21b)$$

Problem (21a) is trivially convex with respect to $u$ for a fixed $d$, whereas problem (21b) is convex provided that $h(d)$ is convex and $g(d)$ convex, both bounded from below, so that it has a unique solution for a given $u$.

### 3 Lip-field regularization

In this section, we focus on the minimization problem with respect to $d$ given by equation (21b), finding a solution $d$ which is in space $L$. Detailed explanation and demonstration can be found in [1]. In 1D, the domain is discretized
where

$$h(d) = \int_{0}^{d} H(\tilde{d}) \, d\tilde{d} \qquad (11)$$

and

$$U = \{u \in H^1(\Omega) : u(0) = 0, u(L) = \dot{\epsilon}_0 L t_{n+1}\} \qquad (12)$$

$$D = \{d \in L^{\infty}(\Omega) : d_n \le d \le 1\} \qquad (13)$$

The way the Dirichlet boundary conditions that appear in (12) are applied is detailed in Section 5. Note that the definition of D automatically ensures the irreversibility of the evolution of d. The predictor $u_p$ is defined by

$$u_p = u_n + \Delta t v_n + \frac{\Delta t^2}{2} (1 - 2\beta) a_n \qquad (14)$$

$$v_p = v_n + (1 - \gamma) \Delta t a_n \qquad (15)$$

with $v$ the velocity and $a$ the acceleration. The numerical parameters $\beta$ and $\gamma$ determine the type of numerical integration (explicit or implicit). Then, $a$ is computed (see detail in Section 5) and gives the final values for $u$ and $v$

$$u = u_p + \beta \Delta t^2 a \qquad (16)$$

$$v = v_p + \gamma \Delta t a \qquad (17)$$

The choice $\beta = 0$ corresponds to the explicit dynamics (central difference method). The new displacement is found by $u = u_p$ and thus will not depend on the coming value of damage. The case $\beta = 1/4$, $\gamma = 1/2$ corresponds to the implicit Newmark scheme (undamped trapezoidal rule). Note that for the explicit scheme, it is common to consider a lumped mass matrix, whereas for the implicit scheme, an exact mass matrix integration is performed. Also, taking the variation of (10b) with respect to $u$ gives the time discretization of the weak form of the equilibrium (2), where the acceleration is obtained by equation (16) :

$$\int_{0}^{L} \sigma \epsilon(u^*) \, dx = \int_{0}^{L} \rho \frac{(u - u_p)}{\beta \Delta t^2} u^* \, dx, \quad \forall u^* \in \mathcal{U}^* \qquad (18)$$

while the variation of (10b) with respect to $d$ and the constrain $d \in D$ gives equations (9).

Solving (10) in space $D$ is an ill-posed problem, which leads to solutions exhibiting spurious mesh dependency. Therefore, the Lip-field approach [1] introduces an alternative problem

$$(u, d) = \arg \min_{\tilde{u} \in U, \tilde{d} \in D \cup L} F(\tilde{u}, \tilde{d}) \qquad (19)$$

where

$$L = \{d \in L^{\infty}([0, L]) : |d(x) - d(y)| \le |x - y|/\ell, \quad \forall x, y \in [0, L]\} \qquad (20)$$

The only difference between problems (10) and (19) is the space where the solution $d$ is computed: by enforcing that $d$ must belong to $L$, we prevent the gradient of $d$ to go to infinity, thus avoiding spurious localization. In practice, this problem is not convex with respect to $(u, d)$, therefore, it is solved with a staggered iterative process²:

$$u^{k+1} = \arg \min_{\tilde{u} \in U} F(\tilde{u}, d^k) \qquad (21a)$$

$$d^{k+1} = \arg \min_{\tilde{d} \in D \cup L} F(u^{k+1}, \tilde{d}) \qquad (21b)$$

Problem (21a) is trivially convex with respect to $u$ for a fixed $d$, whereas problem (21b) is convex provided that $h(d)$ is convex and $g(d)$ convex, both bounded from below, so that it has a unique solution for a given $u$.

### 3 Lip-field regularization

In this section, we focus on the minimization problem with respect to $d$ given by equation (21b), finding a solution $d$ which is in space $L$. Detailed explanation and demonstration can be found in [1]. In 1D, the domain is discretized
<figure><img src="image_1.png" /><figcaption>Figure 1: Discretization of a 1D bar. Displacements are stored at nodes, and damage variables at the integration points, located at center of each element.</figcaption></figure>

with $N_e$ linear finite elements of equal size $h_e$, where the displacement $u$ is stored at each finite element node, while the damage $d$ is stored at the element integration points (element centroids), as illustrated on Figure 1.

In this case, the Lipschitz constraint on the damage variable simply writes

$$d_i - d_{i+1} - h_e/\ell \le 0, \quad i = 1, \dots, N_e - 1 \qquad (22a)$$

$$d_i - d_{i-1} - h_e/\ell \le 0, \quad i = 2, \dots, N_e \qquad (22b)$$

The above condition can be provided to any optimizer to solve (21b). This problem can be solved in the entire computation domain, however [1] proposed an efficient way to reduce the size of the domain where problem (21b) needs to be solved, which is used in this paper. The first step is to solve (21b) without the Lipschitz constraint:

$$\bar{d} = \arg \min_{\bar{d} \in D} F(u^{k+1}, \tilde{d}) \qquad (23)$$

This provides a prediction $\bar{d}$ of the damage field. If it satisfies the Lipschitz constraint, then $d^{k+1} = \bar{d}$ (see Figure 2). Otherwise, it was shown by [1] that $d^{k+1}$ can be bounded by a lower projection $\pi^\ell$ and an upper projection $\pi^u$

$$\pi^\ell \bar{d}(x) = \min_{y \in \Omega} \left( \bar{d}(y) + \frac{1}{\ell} |x - y| \right) \qquad (24a)$$

$$\pi^u \bar{d}(x) = \min_{y \in \Omega} \left( \bar{d}(y) - \frac{1}{\ell} |x - y| \right) \qquad (24b)$$

such that

$$d_n \le \pi^\ell \bar{d} \le \bar{d} \le \pi^u \bar{d} \le 1 \qquad (25)$$

and

$$\pi^\ell \bar{d} \le d^{k+1} \le \pi^u \bar{d} \qquad (26)$$

An important consequence is that

$$\pi^\ell \bar{d}(x) = \pi^u \bar{d}(x) \Rightarrow d^{k+1}(x) = \bar{d}(x) \qquad (27)$$

We can make use of this property to solve (21b) only where $\pi^\ell \bar{d}(x) \ne \pi^u \bar{d}(x)$, thus reducing the size of the problem to be solved.

## 4 Equivalence with a cohesive zone model (CZM)

In this section, we discuss the expression of the softening function, and in particular its primitive $h$, introduced in equation (11). In [75], it was shown on a 1D case that for a particular choice of $h$, it was possible to mimic the behavior of a linear cohesive zone model (CZM) (see Figure 4). The same reasoning is used in this paper.

Considering for $g$ the following expression:

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

The expression of $h$ which allows to mimic a linear CZM is

$$h(d) = \frac{2d - d^2}{(1 - d + \lambda d^2)^2} \qquad (29)$$

where

$$\lambda = \frac{2Y_c \ell}{G_c} \qquad (30)$$

<footnote n="2">²Here, we distinguish between timesteps and iterations, with the former annotated as subscripts (e.g. $d_n$) and the latter as superscripts ($d^k$).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Discretization of a 1D bar. Displacements are stored at nodes, and damage variables at the integration points, located at center of each element.</figcaption></figure>

with $N_e$ linear finite elements of equal size $h_e$, where the displacement $u$ is stored at each finite element node, while the damage $d$ is stored at the element integration points (element centroids), as illustrated on Figure 1.

In this case, the Lipschitz constraint on the damage variable simply writes

$$d_i - d_{i+1} - h_e/\ell \le 0, \quad i = 1, \dots, N_e - 1 \qquad (22a)$$

$$d_i - d_{i-1} - h_e/\ell \le 0, \quad i = 2, \dots, N_e \qquad (22b)$$

The above condition can be provided to any optimizer to solve (21b). This problem can be solved in the entire computation domain, however [1] proposed an efficient way to reduce the size of the domain where problem (21b) needs to be solved, which is used in this paper. The first step is to solve (21b) without the Lipschitz constraint:

$$\bar{d} = \arg \min_{\bar{d} \in D} F(u^{k+1}, \tilde{d}) \qquad (23)$$

This provides a prediction $\bar{d}$ of the damage field. If it satisfies the Lipschitz constraint, then $d^{k+1} = \bar{d}$ (see Figure 2). Otherwise, it was shown by [1] that $d^{k+1}$ can be bounded by a lower projection $\pi^\ell$ and an upper projection $\pi^u$

$$\pi^\ell \bar{d}(x) = \min_{y \in \Omega} \left( \bar{d}(y) + \frac{1}{\ell} |x - y| \right) \qquad (24a)$$

$$\pi^u \bar{d}(x) = \min_{y \in \Omega} \left( \bar{d}(y) - \frac{1}{\ell} |x - y| \right) \qquad (24b)$$

such that

$$d_n \le \pi^\ell \bar{d} \le \bar{d} \le \pi^u \bar{d} \le 1 \qquad (25)$$

and

$$\pi^\ell \bar{d} \le d^{k+1} \le \pi^u \bar{d} \qquad (26)$$

An important consequence is that

$$\pi^\ell \bar{d}(x) = \pi^u \bar{d}(x) \Rightarrow d^{k+1}(x) = \bar{d}(x) \qquad (27)$$

We can make use of this property to solve (21b) only where $\pi^\ell \bar{d}(x) \ne \pi^u \bar{d}(x)$, thus reducing the size of the problem to be solved.

# 4 Equivalence with a cohesive zone model (CZM)

In this section, we discuss the expression of the softening function, and in particular its primitive $h$, introduced in equation (11). In [75], it was shown on a 1D case that for a particular choice of $h$, it was possible to mimic the behavior of a linear cohesive zone model (CZM) (see Figure 4). The same reasoning is used in this paper.

Considering for $g$ the following expression:

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

The expression of $h$ which allows to mimic a linear CZM is

$$h(d) = \frac{2d - d^2}{(1 - d + \lambda d^2)^2} \qquad (29)$$

where

$$\lambda = \frac{2Y_c \ell}{G_c} \qquad (30)$$

<footnote n="2">²Here, we distinguish between timesteps and iterations, with the former annotated as subscripts (e.g. $d_n$) and the latter as superscripts ($d^k$).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: Sketch of computation of local update $d^{k+1}$ from $d_n$ [1]. If $\bar{d}$ is in L, then $d^{k+1} = \bar{d}$. Otherwise, $\pi^l \bar{d}$ and $\pi^l d$ give bounds for the projection of $\bar{d}$ on L.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Illustration of the properties of the upper and lower projections [1]. Outside the interval $[a, b]$, the projections are equal and thus $d^{k+1} = \bar{d}$. Inside $[a, b]$, we have $\pi^l \bar{d} \le d^{k+1} \le \pi^u \bar{d}$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 4: Linear traction-separation law for cohesive zone</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2: Sketch of computation of local update $d^{k+1}$ from $d_n$ [1]. If $\bar{d}$ is in L, then $d^{k+1} = \bar{d}$. Otherwise, $\pi^l \bar{d}$ and $\pi^l d$ give bounds for the projection of $\bar{d}$ on L.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the properties of the upper and lower projections [1]. Outside the interval $[a, b]$, the projections are equal and thus $d^{k+1} = \bar{d}$. Inside $[a, b]$, we have $\pi^l \bar{d} \le d^{k+1} \le \pi^u \bar{d}$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Linear traction-separation law for cohesive zone</figcaption></figure>
where $G_c$ is the material toughness. Note that $h$ is convex provided that $\lambda \le 1/2$. Also, in the framework of the CZM equivalence, the material must start breaking when the stress is equal to the critical stress $\sigma_c$, which imposes

$$Y_c = \frac{\sigma_c^2}{2E} \qquad (31)$$

In what follows, we will simply call “Lip-field” the Lip-field model associated with the above $g(d)$ and $h(d)$ functions.

We have explained how to get a Lip-field model equivalent to a linear CZM<fnref n="3" />, so it would be logical to compare the results obtained with a true linear CZM. We implemented and alternate version of the CZM, inspired by the so-called crack-band model from [76, 77]. It consists in using a damage model depending on the finite element size $h_e$. We consider a purely local damage model, so that $(u, d)$ is solution of the following problem

$$(u, d) = \arg \min_{\tilde{u} \in U, \tilde{d} \in D} F_{\text{CZM}}(\tilde{u}, \tilde{d}) \qquad (32a)$$

$$F_{\text{CZM}}(u, d) = \int_{0}^{L} \frac{1}{2} \rho \frac{(u - u_{\text{p}})^2}{\beta \Delta t^2} + \varphi(\epsilon(u), d) + Y_c h_{\text{CZM}}(d) \, dx \qquad (32b)$$

where

$$h_{\text{CZM}}(d) = \frac{1}{(1 - \lambda_c)} \left( \frac{1}{(1 - \lambda_c)g(d) + \lambda_c} - 1 \right) \qquad (33)$$

with

$$\lambda_c = \frac{\sigma_c h_e}{E w_c} \qquad (34)$$

The first difference with the Lip-field problem described by equations (19) is that $d$ does not need to belong to L anymore. Thus, the resolution algorithm used to solve (19) can also be used to solve (32), just by deactivating the Lipschitz constraint on $d$. The second difference is the softening function $h_{\text{CZM}}$, which is computed so that the behavior of a single finite element of size $h_e$ has the same behavior as a finite element modeled with a CZM. The details of the derivation of $h_{\text{CZM}}$ can be found in the appendix. Note that a 1D bar finite element described with the equations above has analytically the same “stress-displacement at the ends of the bar” behavior as a 1D finite element modeled using a classical CZM. The only difference is the absence of displacement jump on the finite element, which is in 1D not compulsory if the only quantity of interested considered is the dissipated energy, which is the case in this paper. This alleviates the need to introduce new displacement degrees of freedom as cracks appear, and makes the implementation more straightforward.

Note that in problem (32), $d$ does not belong to L, so its gradient may eventually go to infinity. However, because of the dependency of $h_{\text{CZM}}(d)$ on the finite element size $h_e$, the dissipated energy does not go to zero when the mesh is refined.

In what follows, we will simply call “CZM” the local damage model associated with the above $h_{\text{CZM}}(d)$ function.

# 5 Numerical resolution

The displacement, speed, and acceleration fields are discretized using classical linear interpolation functions:

$$u(x, t) = \sum_{i \in N} u_i(t) N_i(x) \qquad (35)$$

$$v(x, t) = \sum_{i \in N} v_i(t) N_i(x) \qquad (36)$$

$$a(x, t) = \sum_{i \in N} a_i(t) N_i(x) \qquad (37)$$

where $N$ is the set of the finite element nodes, as illustrated on Figure 1, $u_i$, $v_i$, $a_i$ are the degrees of freedom associated to node $i$ and $N_i$ the corresponding interpolation function. In what follows, we will note $(U_n, V_n, A_n)$ the vectors

<footnote n="3">³We note that the Lip-field model derived herein is equivalent to a linear CZM in the sense that it features the same parameters as a linear CZM ({G_c, σ_c}) and produces the same crack-opening behavior in 1D quasistatics for a single flaw. This is analogous to the derivation for TLS-CZM equivalence derived in [75]. For general application, such as the present focus on 1D dynamics with multiple opening flaws, the behavior is expected to be similar, but not to correspond absolutely.</footnote>
where $G_c$ is the material toughness. Note that $h$ is convex provided that $\lambda \le 1/2$. Also, in the framework of the CZM equivalence, the material must start breaking when the stress is equal to the critical stress $\sigma_c$, which imposes

$$Y_c = \frac{\sigma_c^2}{2E} \qquad (31)$$

In what follows, we will simply call “Lip-field” the Lip-field model associated with the above $g(d)$ and $h(d)$ functions.

We have explained how to get a Lip-field model equivalent to a linear CZM<fnref n="3" />, so it would be logical to compare the results obtained with a true linear CZM. We implemented and alternate version of the CZM, inspired by the so-called crack-band model from [76, 77]. It consists in using a damage model depending on the finite element size $h_e$. We consider a purely local damage model, so that $(u, d)$ is solution of the following problem

$$(u, d) = \arg \min_{\tilde{u} \in U, \tilde{d} \in D} F_{\text{CZM}}(\tilde{u}, \tilde{d}) \qquad (32a)$$

$$F_{\text{CZM}}(u, d) = \int_{0}^{L} \frac{1}{2} \rho \frac{(u - u_{\text{p}})^2}{\beta \Delta t^2} + \varphi(\epsilon(u), d) + Y_c h_{\text{CZM}}(d) \, dx \qquad (32b)$$

where

$$h_{\text{CZM}}(d) = \frac{1}{(1 - \lambda_c)} \left( \frac{1}{(1 - \lambda_c)g(d) + \lambda_c} - 1 \right) \qquad (33)$$

with

$$\lambda_c = \frac{\sigma_c h_e}{E w_c} \qquad (34)$$

The first difference with the Lip-field problem described by equations (19) is that $d$ does not need to belong to L anymore. Thus, the resolution algorithm used to solve (19) can also be used to solve (32), just by deactivating the Lipschitz constraint on $d$. The second difference is the softening function $h_{\text{CZM}}$, which is computed so that the behavior of a single finite element of size $h_e$ has the same behavior as a finite element modeled with a CZM. The details of the derivation of $h_{\text{CZM}}$ can be found in the appendix. Note that a 1D bar finite element described with the equations above has analytically the same “stress-displacement at the ends of the bar” behavior as a 1D finite element modeled using a classical CZM. The only difference is the absence of displacement jump on the finite element, which is in 1D not compulsory if the only quantity of interested considered is the dissipated energy, which is the case in this paper. This alleviates the need to introduce new displacement degrees of freedom as cracks appear, and makes the implementation more straightforward.

Note that in problem (32), $d$ does not belong to L, so its gradient may eventually go to infinity. However, because of the dependency of $h_{\text{CZM}}(d)$ on the finite element size $h_e$, the dissipated energy does not go to zero when the mesh is refined.

In what follows, we will simply call “CZM” the local damage model associated with the above $h_{\text{CZM}}(d)$ function.

# 5 Numerical resolution

The displacement, speed, and acceleration fields are discretized using classical linear interpolation functions:

$$u(x, t) = \sum_{i \in N} u_i(t) N_i(x) \qquad (35)$$

$$v(x, t) = \sum_{i \in N} v_i(t) N_i(x) \qquad (36)$$

$$a(x, t) = \sum_{i \in N} a_i(t) N_i(x) \qquad (37)$$

where $N$ is the set of the finite element nodes, as illustrated on Figure 1, $u_i$, $v_i$, $a_i$ are the degrees of freedom associated to node $i$ and $N_i$ the corresponding interpolation function. In what follows, we will note $(U_n, V_n, A_n)$ the vectors

<footnote n="3">³We note that the Lip-field model derived herein is equivalent to a linear CZM in the sense that it features the same parameters as a linear CZM ({G_c, σ_c}) and produces the same crack-opening behavior in 1D quasistatics for a single flaw. This is analogous to the derivation for TLS-CZM equivalence derived in [75]. For general application, such as the present focus on 1D dynamics with multiple opening flaws, the behavior is expected to be similar, but not to correspond absolutely.</footnote>
containing the degrees of freedom of $(u_n, v_n, a_n)$. They follow the same naming convention (predicted values, iterate, values at $t$, etc.). By injecting this discretization in equation (2), we get the following matrix system:

$$[\mathbf{M}]{\{\mathbf{A}_n\} + [\mathbf{K}(d_n)]}{\{\mathbf{U}_n\} = \{\mathbf{0}\}} \qquad (38)$$

where

$$M_{ij} = \int_{0}^{L} \rho N_i N_j \, dx \quad \forall (i, j) \in N \times N \qquad (39)$$

$$K_{ij}(d_n) = \int_{0}^{L} E g(d_n) \rho N_{i,x} N_{j,x} \, dx \quad \forall (i, j) \in N \times N \qquad (40)$$

Now we discuss how to solve (38). As stated in Section 2, the computation of $u$ depends on the values of the parameters $\beta$ and $\gamma$.

**Explicit dynamic :** $\beta = 0$ and $\gamma = 1/2$. If $\beta = 0$, one can see from equation (10b) that the only solution of the problem is

$$\{\mathbf{U}\} = \{\mathbf{U}_p\} = \{\mathbf{U}_n\} + \Delta t \{\mathbf{V}_n\} + \frac{\Delta t^2}{2} \{\mathbf{A}_n\} \qquad (41)$$

The predicted nodal velocity vector is

$$\{\mathbf{V}_p\} = \{\mathbf{V}_n\} + \frac{\Delta t}{2} \{\mathbf{A}_n\} \qquad (42)$$

Knowing $\{\mathbf{U}\}$, problem (21b) can be solved to compute $d$. Then, the nodal acceleration vector $\{\mathbf{A}\}$ is computed by solving

$$[\tilde{\mathbf{M}}]{\{\mathbf{A}\} = -[\mathbf{K}(d)]\{\mathbf{U}\}} \qquad (43)$$

where $[\tilde{\mathbf{M}}]$ is the lumped version of operator $[\mathbf{M}]$. To impose the Dirichlet boundary conditions at $x = 0$ and $x = L$, the first and last degrees of freedom of $\{\mathbf{A}_n\}$ are set to 0 (which gives the expected results, since their velocity is imposed at the beginning of the simulation from the initial conditions, and the velocity at both tips of the bar is constant.). Finally, nodal velocity vector $\{\mathbf{V}\}$ is obtained by equation (17):

$$\{\mathbf{V}\} = \{\mathbf{V}_p\} + \frac{\Delta t}{2} \{\mathbf{A}\} \qquad (44)$$

The explicit dynamics resolution is summarized in Algorithm 1.

<caption>**Algorithm 1:** Explicit dynamic resolution</caption>

1 Initialization: $n = 0$, $\{\mathbf{U}_0\}$, $\{\mathbf{V}_0\}$, $\{\mathbf{A}_0\}$, $d_0$;
2 while $n < n_{max}$ do
3 Compute nodal displacement $\{\mathbf{U}\}$ using (41) and predicted nodal velocity $\{\mathbf{V}_p\}$ using (42);
4 Solve (23) to compute damage prediction $\bar{d}$;
5 Compute projections $\pi^u \bar{d}$ and $\pi^l \bar{d}$ using (24);
6 Solve (21b) where $\pi^u \bar{d} \neq \pi^l \bar{d}$ to find $d$;
7 Compute nodal acceleration $\{\mathbf{A}\}$ by solving (43);
8 Set first and last values of $\{\mathbf{A}\}$ to 0;
9 Compute nodal velocity using $\{\mathbf{V}\}$ (44);
10 $n \leftarrow n + 1$;

**Implicit dynamic :** $\beta = 1/2$ and $\gamma = 1/4$. The first step is to compute the vector of nodal displacement prediction $\{\mathbf{U}_p\}$ using (14)

$$\{\mathbf{U}_p\} = \{\mathbf{U}_n\} + \Delta t \{\mathbf{V}_n\} + \frac{\Delta t^2}{2}(1 - 2\beta) \{\mathbf{A}_n\} \qquad (45)$$

Then, unlike the explicit dynamics case, iterations are necessary to compute $\{\mathbf{U}\}$ and $d$. By injecting the expression of $u_p$ in (21a), we get the linear system that needs to be solved to find $\{\mathbf{U}^{k+1}\}$:

$$(\beta \Delta t^2 [\mathbf{K}(d^k)] + [\mathbf{M}]) \{\mathbf{U}^{k+1}\} = [\mathbf{M}] \{\mathbf{U}_p\} \qquad (46)$$
containing the degrees of freedom of $(u_n, v_n, a_n)$. They follow the same naming convention (predicted values, iterate, values at $t$, etc.). By injecting this discretization in equation (2), we get the following matrix system:

$$[\mathbf{M}]{\{\mathbf{A}_n\} + [\mathbf{K}(d_n)]}{\{\mathbf{U}_n\} = \{\mathbf{0}\}} \qquad (38)$$

where

$$M_{ij} = \int_{0}^{L} \rho N_i N_j \, dx \quad \forall (i, j) \in N \times N \qquad (39)$$

$$K_{ij}(d_n) = \int_{0}^{L} E g(d_n) \rho N_{i,x} N_{j,x} \, dx \quad \forall (i, j) \in N \times N \qquad (40)$$

Now we discuss how to solve (38). As stated in Section 2, the computation of $u$ depends on the values of the parameters $\beta$ and $\gamma$.

**Explicit dynamic :** $\beta = 0$ and $\gamma = 1/2$. If $\beta = 0$, one can see from equation (10b) that the only solution of the problem is

$$\{\mathbf{U}\} = \{\mathbf{U}_p\} = \{\mathbf{U}_n\} + \Delta t \{\mathbf{V}_n\} + \frac{\Delta t^2}{2} \{\mathbf{A}_n\} \qquad (41)$$

The predicted nodal velocity vector is

$$\{\mathbf{V}_p\} = \{\mathbf{V}_n\} + \frac{\Delta t}{2} \{\mathbf{A}_n\} \qquad (42)$$

Knowing $\{\mathbf{U}\}$, problem (21b) can be solved to compute $d$. Then, the nodal acceleration vector $\{\mathbf{A}\}$ is computed by solving

$$[\tilde{\mathbf{M}}]{\{\mathbf{A}\}} = -[\mathbf{K}(d)]{\{\mathbf{U}\}} \qquad (43)$$

where $[\tilde{\mathbf{M}}]$ is the lumped version of operator $[\mathbf{M}]$. To impose the Dirichlet boundary conditions at $x = 0$ and $x = L$, the first and last degrees of freedom of $\{\mathbf{A}_n\}$ are set to 0 (which gives the expected results, since their velocity is imposed at the beginning of the simulation from the initial conditions, and the velocity at both tips of the bar is constant.). Finally, nodal velocity vector $\{\mathbf{V}\}$ is obtained by equation (17):

$$\{\mathbf{V}\} = \{\mathbf{V}_p\} + \frac{\Delta t}{2} \{\mathbf{A}\} \qquad (44)$$

The explicit dynamics resolution is summarized in Algorithm 1.

<caption>**Algorithm 1:** Explicit dynamic resolution</caption>

1 Initialization: $n = 0$, $\{\mathbf{U}_0\}$, $\{\mathbf{V}_0\}$, $\{\mathbf{A}_0\}$, $d_0$;
2 **while** $n < n_{max}$ **do**
3 &nbsp;&nbsp;&nbsp;&nbsp;Compute nodal displacement $\{\mathbf{U}\}$ using (41) and predicted nodal velocity $\{\mathbf{V}_p\}$ using (42);
4 &nbsp;&nbsp;&nbsp;&nbsp;Solve (23) to compute damage prediction $\bar{d}$;
5 &nbsp;&nbsp;&nbsp;&nbsp;Compute projections $\pi^u \bar{d}$ and $\pi^l \bar{d}$ using (24);
6 &nbsp;&nbsp;&nbsp;&nbsp;Solve (21b) where $\pi^u \bar{d} \neq \pi^l \bar{d}$ to find $d$;
7 &nbsp;&nbsp;&nbsp;&nbsp;Compute nodal acceleration $\{\mathbf{A}\}$ by solving (43);
8 &nbsp;&nbsp;&nbsp;&nbsp;Set first and last values of $\{\mathbf{A}\}$ to 0;
9 &nbsp;&nbsp;&nbsp;&nbsp;Compute nodal velocity using $\{\mathbf{V}\}$ (44);
10 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$n \leftarrow n + 1$;

**Implicit dynamic :** $\beta = 1/2$ and $\gamma = 1/4$. The first step is to compute the vector of nodal displacement prediction $\{\mathbf{U}_p\}$ using (14)

$$\{\mathbf{U}_p\} = \{\mathbf{U}_n\} + \Delta t \{\mathbf{V}_n\} + \frac{\Delta t^2}{2} (1 - 2\beta) \{\mathbf{A}_n\} \qquad (45)$$

Then, unlike the explicit dynamics case, iterations are necessary to compute $\{\mathbf{U}\}$ and $d$. By injecting the expression of $u_p$ in (21a), we get the linear system that needs to be solved to find $\{\mathbf{U}^{k+1}\}$:

$$(\beta \Delta t^2 [\mathbf{K}(d^k)] + [\mathbf{M}]) \{\mathbf{U}^{k+1}\} = [\mathbf{M}] \{\mathbf{U}_p\} \qquad (46)$$
Here the boundary conditions are classically applied on {$U^{k+1}$}, by proper modifications of the matrix and right-hand side of (46). Problem (21b) is solved to compute $d^k$. Alternate resolution of (46) and (21b) are performed until convergence. Finally, the nodal acceleration vector is obtained by replacing in (16) the expression of $u_p$ from equation (14):

$$ \{\mathbf{A}\} = \frac{1}{\beta\Delta t^2} (\{\mathbf{U}\} - \{\mathbf{U}_p\}) - \frac{1}{\beta\Delta t} \{\mathbf{V}_n\} - \left(\frac{1}{2\beta} - 1\right) \{\mathbf{A}_n\} \quad (47) $$

and the nodal velocity vector by equation (17)

$$ \{\mathbf{V}\} = \{\mathbf{V}_n\} + \frac{1}{2}(1-\gamma)\Delta t\{\mathbf{A}_n\} + \gamma\Delta t\{\mathbf{A}\} \quad (48) $$

The implicit dynamic resolution is summarized in Algorithm 2.

<caption>**Algorithm 2: Implicit dynamic resolution**</caption>

1 Initialization: $n = 0$, $\{U_0\}$, $\{V_0\}$, $\{A_0\}$, $d_0$;
2 **while** $n < n_{max}$ **do**
3 Compute nodal displacement prediction $\{U_p\}$ using (45);
4 $k = 0$;
5 $d^0 = d_n$;
6 **while** $k < k_{max}$ **and** $err_u > tol_u$ **and** $err_d > tol_d$ **do**
7 Assemble matrix and right-hand side of (46);
8 Modify matrix and right-hand side of (46) to apply boundary conditions;
9 Compute nodal displacement $\{U^{k+1}\}$ by solving (46);
10 Solve (23) to compute damage prediction $\bar{d}^{k+1}$;
11 Compute projections $\pi^u \bar{d}^{k+1}$ and $\pi^l \bar{d}^{k+1}$ using (24);
12 Solve (21b) where $\pi^u \bar{d}^{k+1} \neq \pi^l \bar{d}^{k+1}$ to find $d^{k+1}$;
13 Compute errors:
$$ err_u = \left| \frac{\{U^{k+1}\} - \{U^k\}}{\{U^{k+1}\} - \{U_n\}} \right|, \quad err_d = \left| \frac{d^{k+1} - d^k}{d^{k+1} - d_n} \right| \quad (49) $$
$$ k \leftarrow k + 1; $$
14 Compute nodal acceleration $\{A\}$ using (47) and nodal velocity $\{V\}$ using (48);
15 $n \leftarrow n + 1$;

# 6 Numerical examples

In this section, we present the results obtained with the Lip-field and the CZM. The considered example is a 1D bar of length $L = 2 \cdot 10^{-3}$ m and cross-sectional area $A = 2 \cdot 10^{-7}$ m² from [9, 7] made of dense alumina (Al₂O₃). The numerical values of the material parameters are given in Table 1. Following the reasoning of [69], the Lip-field regularization length is set to $\ell = 2.21 \cdot 10^{-6}$ m, which represents half the size of the expected fragment size from Drugan's analysis [7] for the highest strain rate of $\dot{\epsilon}_0 = 7.5 \cdot 10^6$ s⁻¹ tested in [69]. Both will be tested with explicit and implicit dynamic resolution, with a time step of $\Delta t = 0.99\Delta t_c$, where $\Delta t_c = h_e/c$ with $c = \sqrt{E/\rho}$ the elastic wave speed in the same material. Note that the critical time step increases when damage increases, so using a time step smaller than $\Delta t_c$ ensures the stability of explicit dynamics resolution. The values of the numerical parameters for the implicit dynamic resolution are $k_{max} = 10$, $tol_u = 10^{-6}$, $tol_d = 10^{-6}$.

A stochastic approach is employed in order to hasten mesh convergence of global quantities in the otherwise-uniform fragmentation simulations, following the observations of Molinari et al. [78] (random mesh spacing) and Stershic et al. [69] (random material properties). Here, the Young's modulus follows a Weibull distribution, which can be obtained as explained in [69] by:

$$ E(r) = E_0(-\ln(r))^{1/m} + E_{min} \quad (50) $$

<footnote n="4">⁴Reproducing the approach of [69], an effective Young's modulus is employed in the one-dimensional finite element analysis to capture the plane-strain conditions from [9]. Note that [69] used an uniaxial-strain modulus, whereas [9] presents a plane-strain geometry.</footnote>
Here the boundary conditions are classically applied on {$U^{k+1}$}, by proper modifications of the matrix and right-hand side of (46). Problem (21b) is solved to compute $d^k$. Alternate resolution of (46) and (21b) are performed until convergence. Finally, the nodal acceleration vector is obtained by replacing in (16) the expression of $u_p$ from equation (14):

$$ \{\mathbf{A}\} = \frac{1}{\beta\Delta t^2} (\{\mathbf{U}\} - \{\mathbf{U}_p\}) - \frac{1}{\beta\Delta t} \{\mathbf{V}_n\} - \left(\frac{1}{2\beta} - 1\right) \{\mathbf{A}_n\} \quad (47) $$

and the nodal velocity vector by equation (17)

$$ \{\mathbf{V}\} = \{\mathbf{V}_n\} + \frac{1}{2}(1 - \gamma)\Delta t\{\mathbf{A}_n\} + \gamma\Delta t\{\mathbf{A}\} \quad (48) $$

The implicit dynamic resolution is summarized in Algorithm 2.

<caption>**Algorithm 2:** Implicit dynamic resolution</caption>

1 Initialization: $n = 0$, $\{U_0\}$, $\{V_0\}$, $\{A_0\}$, $d_0$;
2 **while** $n < n_{max}$ **do**
3 Compute nodal displacement prediction $\{U_p\}$ using (45);
4 $k = 0$;
5 $d^0 = d_n$;
6 **while** $k < k_{max}$ **and** $err_u > tol_u$ **and** $err_d > tol_d$ **do**
7 Assemble matrix and right-hand side of (46);
8 Modify matrix and right-hand side of (46) to apply boundary conditions;
9 Compute nodal displacement $\{U^{k+1}\}$ by solving (46);
10 Solve (23) to compute damage prediction $\bar{d}^{k+1}$;
11 Compute projections $\pi^u \bar{d}^{k+1}$ and $\pi^l \bar{d}^{k+1}$ using (24);
12 Solve (21b) where $\pi^u \bar{d}^{k+1} \neq \pi^l \bar{d}^{k+1}$ to find $d^{k+1}$;
13 Compute errors:
$$ err_u = \left| \frac{\{U^{k+1}\} - \{U^k\}}{\{U^{k+1}\} - \{U_n\}} \right|, \quad err_d = \left| \frac{d^{k+1} - d^k}{d^{k+1} - d_n} \right| \quad (49) $$
14 Compute nodal acceleration $\{A\}$ using (47) and nodal velocity $\{V\}$ using (48);
15 $n \leftarrow n + 1$;

# 6 Numerical examples

In this section, we present the results obtained with the Lip-field and the CZM. The considered example is a 1D bar of length $L = 2 \cdot 10^{-3}$ m and cross-sectional area $A = 2 \cdot 10^{-7}$ m² from [9, 7] made of dense alumina (Al₂O₃). The numerical values of the material parameters are given in Table 1. Following the reasoning of [69], the Lip-field regularization length is set to $\ell = 2.21 \cdot 10^{-6}$ m, which represents half the size of the expected fragment size from Drugan's analysis [7] for the highest strain rate of $\dot{\epsilon}_0 = 7.5 \cdot 10^6$ s⁻¹ tested in [69]. Both will be tested with explicit and implicit dynamic resolution, with a time step of $\Delta t = 0.99\Delta t_c$, where $\Delta t_c = h_e/c$ with $c = \sqrt{E/\rho}$ the elastic wave speed in the same material. Note that the critical time step increases when damage increases, so using a time step smaller than $\Delta t_c$ ensures the stability of explicit dynamics resolution. The values of the numerical parameters for the implicit dynamic resolution are $k_{max} = 10$, $tol_u = 10^{-6}$, $tol_d = 10^{-6}$.

A stochastic approach is employed in order to hasten mesh convergence of global quantities in the otherwise-uniform fragmentation simulations, following the observations of Molinari et al. [78] (random mesh spacing) and Stershic et al. [69] (random material properties). Here, the Young's modulus follows a Weibull distribution, which can be obtained as explained in [69] by:

$$ E(r) = E_0(-\ln(r))^{1/m} + E_{min} \quad (50) $$

<footnote n="4">⁴Reproducing the approach of [69], an effective Young's modulus is employed in the one-dimensional finite element analysis to capture the plane-strain conditions from [9]. Note that [69] used an uniaxial-strain modulus, whereas [9] presents a plane-strain geometry.</footnote>
<table><caption>Table 1: Material properties for numerical example</caption><thead><tr><th>Properties</th><th>Units</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Density</td><td>kg·m<sup>-3</sup></td><td>ρ</td><td>3.9 · 10<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>Pa</td><td>E</td><td>380 · 10<sup>9</sup></td></tr><tr><td>Fracture toughness</td><td>N/m</td><td>G<sub>c</sub></td><td>83.13</td></tr><tr><td>Critical energy release rate</td><td>Pa</td><td>Y<sub>c</sub></td><td>820 · 10<sup>3</sup></td></tr><tr><td>Critical tensile stress</td><td>Pa</td><td>σ<sub>c</sub></td><td>1.0 · 10<sup>9</sup></td></tr></tbody></table>

Where $r \sim U(0, 1)$ is a uniformy distributed random variable. The value of $m$ is set to 2, which corresponds to a Rayleigh distribution. The values of $E_0$ and $E_{min}$ are computed to get an average value equal to $E$ and a coefficient of variation equal to $CV = 0.01$, which gives

$$E_{min} = E(1.0 - 1.9130584 \cdot CV), \quad E_0 = 2.1586552 \cdot E \cdot CV \quad (51)$$

Finally, the minimization problem (21b) is solved with the `scipyminimize Python routine, using Sequential Least Squares Programming.

In Figure 5, snapshots of the damage fields for the Lip-field and the CZM are given at different time steps, for $\dot{e}_0 = 1 \cdot 10^5$ s$^{-1}$ and $h_e = \ell/10$. For the Lip-field, after a first stage where the damage evolves in a purely local way, following the variations of the Young modulus, the Lip-field regularization activates itself in the zones where damage localizes, materialized by the black, thick horizontal lines. The number and size of these zones increase, then decrease until their number corresponds to the final number of cracks. One can observe that the CZM can create cracks with only one finite element. On the other hand, the Lip-field tends to spread damage on zones with a width of about 20 finite elements (see for instance Figure 5 (g)), which is equal to $2\ell$.

Note that in the case of the CZM, even if damage tends to localize in bands of one-element width, the dissipated energy does not tend to zero with mesh refinement, as illustrated on Figure 6. This is of course expected from the crack-band type model used for the CZM. On Figure 6 are also plotted the results obtained with the different approaches for different mesh element sizes. On the left is plotted the dissipated energy $D$, computed as:

$$D = A \int_{0}^{L} Y_{c}h(d) \, dd \quad (52)$$

On the right is plotted the mean fragment size, obtained by taking the average value of the distance between two cracks. A crack location is defined as the position of the centroid of any element where $d > 0.98$. Each data point of Figure 6 has been obtained by taking the average of 20 runs with different random distributions of the Young modulus. A first observation is that none of the different computed dissipated energy tends to zero with mesh refinement, which is due to the dependency on the mesh element size of the softening function for the CZM, and due to the Lip-field regularization for the Lip-field. The convergence speed of the CZM is rather fast compared to the Lip-field, with results which are almost independent from the mesh refinement, while the dissipation of the Lip-field stabilizes for a mesh element size of about $h_e \approx \ell/35$.

Similar behavior was observed in previous studies for a Thick Level Set (TLS) model equivalent to a linear CZM, where fine meshes were necessary to reach convergence [69, 79]. Comparable studies of dynamic fragmentation using phase-field fracture modeling [62, 63] have also been performed; however, mesh-convergence rate was not reported in detail.

A solution [79] was proposed to improve the TLS convergence rate, which consisted in introducing cohesive zones prior to having the bulk damage variable going to one, thus avoiding strains going to infinite values in the bulk while preserving the amount of dissipated energy. This solution applied to the Lip-field approach would probably help increasing convergence speed, but is out of the scope of this paper.

The results obtained with a mesh element size of $h_e = \ell/10$ for different strain rates are compared to different references on Figure 7. The four tested approaches are able to reproduce the two regimes captured by the different references: a dissipated energy and mean fragment size that is strongly rate-dependent, while transitioning toward rate-independence at lower strain rates. The results obtained with the Lip-field match quite well the results obtained with the CZM. For high strain rates, the results obtained with the explicit and implicit approaches are close to each other, while some differences can be observed for low strain rates. It can be explained by the possible non-uniqueness of the solution of problem (19) with respect to the variable $(u, d)$, due to its non-convexity. Indeed, in explicit dynamics, a unique $u$ is computed for a given $u_n$, and therefore a unique $d$. For a given couple $(u_n, d_n)$ it gives a unique couple $(u, d)$ which is one of the solutions of (19). On the contrary, in implicit dynamic, $(u, d)$ is obtained by alternate
<table><caption>Table 1: Material properties for numerical example</caption><thead><tr><th>Properties</th><th>Units</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Density</td><td>kg·m<sup>-3</sup></td><td>ρ</td><td>3.9 · 10<sup>3</sup></td></tr><tr><td>Young's modulus</td><td>Pa</td><td>E</td><td>380 · 10<sup>9</sup></td></tr><tr><td>Fracture toughness</td><td>N/m</td><td>G<sub>c</sub></td><td>83.13</td></tr><tr><td>Critical energy release rate</td><td>Pa</td><td>Y<sub>c</sub></td><td>820 · 10<sup>3</sup></td></tr><tr><td>Critical tensile stress</td><td>Pa</td><td>σ<sub>c</sub></td><td>1.0 · 10<sup>9</sup></td></tr></tbody></table>

Where $r \sim U(0, 1)$ is a uniformy distributed random variable. The value of $m$ is set to 2, which corresponds to a Rayleigh distribution. The values of $E_0$ and $E_{min}$ are computed to get an average value equal to $E$ and a coefficient of variation equal to $CV = 0.01$, which gives

$$E_{min} = E(1.0 - 1.9130584 \cdot CV), \quad E_0 = 2.1586552 \cdot E \cdot CV \qquad (51)$$

Finally, the minimization problem (21b) is solved with the `scipyoptimize.minimize` Python routine, using Sequential Least Squares Programming.

In Figure 5, snapshots of the damage fields for the Lip-field and the CZM are given at different time steps, for $\dot{e}_0 = 1 \cdot 10^5$ s$^{-1}$ and $h_e = \ell/10$. For the Lip-field, after a first stage where the damage evolves in a purely local way, following the variations of the Young modulus, the Lip-field regularization activates itself in the zones were damage localizes, materialized by the black, thick horizontal lines. The number and size of these zones increase, then decrease until their number corresponds to the final number of cracks. One can observe that the CZM can create cracks with only one finite element. On the other hand, the Lip-field tends to spread damage on zones with a width of about 20 finite elements (see for instance Figure 5 (g)), which is equal to $2\ell$.

Note that in the case of the CZM, even if damage tends to localize in bands of one-element width, the dissipated energy does not tend to zero with mesh refinement, as illustrated on Figure 6. This is of course expected from the crack-band type model used for the CZM. On Figure 6 are also plotted the results obtained with the different approaches for different mesh element sizes. On the left is plotted the dissipated energy $D$, computed as:

$$D = A \int_{0}^{L} Y_{c}h(d) \, dd \qquad (52)$$

On the right is plotted the mean fragment size, obtained by taking the average value of the distance between two cracks. A crack location is defined as the position of the centroid of any element where $d > 0.98$. Each data point of Figure 6 has been obtained by taking the average of 20 runs with different random distributions of the Young modulus. A first observation is that none of the different computed dissipated energy tends to zero with mesh refinement, which is due to the dependency on the mesh element size of the softening function for the CZM, and due to the Lip-field regularization for the Lip-field. The convergence speed of the CZM is rather fast compared to the Lip-field, with results which are almost independent from the mesh refinement, while the dissipation of the Lip-field stabilizes for a mesh element size of about $h_e \approx \ell/35$.

Similar behavior was observed in previous studies for a Thick Level Set (TLS) model equivalent to a linear CZM, where fine meshes were necessary to reach convergence [69, 79]. Comparable studies of dynamic fragmentation using phase-field fracture modeling [62, 63] have also been performed; however, mesh-convergence rate was not reported in detail.

A solution [79] was proposed to improve the TLS convergence rate, which consisted in introducing cohesive zones prior to having the bulk damage variable going to one, thus avoiding strains going to infinite values in the bulk while preserving the amount of dissipated energy. This solution applied to the Lip-field approach would probably help increasing convergence speed, but is out of the scope of this paper.

The results obtained with a mesh element size of $h_e = \ell/10$ for different strain rates are compared to different references on Figure 7. The four tested approaches are able to reproduce the two regimes captured by the different references: a dissipated energy and mean fragment size that is strongly rate-dependent, while transitioning toward rate-independence at lower strain rates. The results obtained with the Lip-field match quite well the results obtained with the CZM. For high strain rates, the results obtained with the explicit and implicit approaches are close to each other, while some differences can be observed for low strain rates. It can be explained by the possible non-uniqueness of the solution of problem (19) with respect to the variable $(u, d)$, due to its non-convexity. Indeed, in explicit dynamics, a unique $u$ is computed for a given $u_n$, and therefore a unique $d$. For a given couple $(u_n, d_n)$ it gives a unique couple $(u, d)$ which is one of the solutions of (19). On the contrary, in implicit dynamic, $(u, d)$ is obtained by alternate
<figure><img src="image_5.png" /><figcaption>Figure 5: Snapshots of damage field for the Lip-field and CZM in explicit dynamics, for $\dot{\epsilon}_0 = 1 \cdot 10^5$ s$^{-1}$ and $h_e = \ell/10$. Black thick horizontal lines indicate where the Lipschitz constraint is active.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Snapshots of damage field for the Lip-field and CZM in explicit dynamics, for $\dot{\epsilon}_0 = 1 \cdot 10^5$ s$^{-1}$ and $h_e = \ell/10$. Black thick horizontal lines indicate where the Lipschitz constraint is active.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 6: Lip-field and CZM convergence study for $\dot{e}_0 = 7.5 \cdot 10^6$ s$^{-1}$. Dissipated energy is plotted on the left and mean fragment size on the right.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 7: Strain rate dependency for $h_e = \ell/10$. Comparison with Drugan [7], Zhou et al. [11], Grady [5], Glenn and Chudnovsky [6] and Stershic [69]. Dissipated energy is plotted on the left and mean fragment size on the right. Error bars indicate the standard deviation of the plotted quantity among the stochastic runs. Note that the magnitude of the error bars on the dissipated energy plot is too small to be visible.</figcaption></figure>

minimization (problem 21). This solution is one solution of (21) among others, which is not necessarily the same as the one obtained in explicit dynamics. It may possible to obtain another one by changing the resolution method used, changing the initial value of the initial iterate $d^0$ at each time step, etc. Different methods were tried when writing this paper, but none of them matched the explicit dynamics results, so they are not presented in this paper. However, the differences observed are within the range of the differences observed between the different references, so telling which solution is better between explicit and implicit dynamic is difficult. Finally, the damage field at computation end, that is, after a time long enough for the dissipated energy to stop increasing, is plotted on Figure 8 for different strain rates, for one particular realization of the Young Modulus.

The dependency of the results with $\ell$ is studied on Figure 9. The dissipated and mean fragment sizes are plotted for the Lip-field in explicit dynamics, for two different strain rates and different values of $\ell$. Noting $\ell_0 = 2.21 \cdot 10^{-6}$ m the value of $\ell$ which was used until now, the values $\ell = \ell_0$, $\ell = \ell_0/2$, $\ell = \ell_0/3$, $\ell = \ell_0/4$ are tested on two different cases

• Case A: The mesh element size evolves with $\ell$, so that the ratio $\frac{h_e}{\ell}$ is fixed and equal to $\frac{1}{5}$. It allows to keep a fixed resolution of the non-local zones.

• Case B: The mesh element size is equal to $\frac{\ell_0}{20}$, whatever the value of $\ell$.
<figure><img src="image_7.png" /><figcaption>Figure 6: Lip-field and CZM convergence study for $\dot{e}_0 = 7.5 \cdot 10^6$ s$^{-1}$. Dissipated energy is plotted on the left and mean fragment size on the right.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Strain rate dependency for $h_e = \ell/10$. Comparison with Drugan [7], Zhou et al. [11], Grady [5], Glenn and Chudnovsky [6] and Stershic [69]. Dissipated energy is plotted on the left and mean fragment size on the right. Error bars indicate the standard deviation of the plotted quantity among the stochastic runs. Note that the magnitude of the error bars on the dissipated energy plot is too small to be visible.</figcaption></figure>

minimization (problem 21). This solution is one solution of (21) among others, which is not necessarily the same as the one obtained in explicit dynamics. It may possible to obtain another one by changing the resolution method used, changing the initial value of the initial iterate $d^0$ at each time step, etc. Different methods were tried when writing this paper, but none of them matched the explicit dynamics results, so they are not presented in this paper. However, the differences observed are within the range of the differences observed between the different references, so telling which solution is better between explicit and implicit dynamic is difficult. Finally, the damage field at computation end, that is, after a time long enough for the dissipated energy to stop increasing, is plotted on Figure 8 for different strain rates, for one particular realization of the Young Modulus.

The dependency of the results with $\ell$ is studied on Figure 9. The dissipated and mean fragment sizes are plotted for the Lip-field in explicit dynamics, for two different strain rates and different values of $\ell$. Noting $\ell_0 = 2.21 \cdot 10^{-6}$ m the value of $\ell$ which was used until now, the values $\ell = \ell_0$, $\ell = \ell_0/2$, $\ell = \ell_0/3$, $\ell = \ell_0/4$ are tested on two different cases

* Case A: The mesh element size evolves with $\ell$, so that the ratio $h_e/\ell$ is fixed and equal to $1/5$. It allows to keep a fixed resolution of the non-local zones.

* Case B: The mesh element size is equal to $\ell_0/20$, whatever the value of $\ell$.
<figure><img src="image_8.png" /><figcaption>Figure 8: Snapshots of the damage field for the Lip-field in explicit dynamics, for $h_e = \ell/10$, for different values of the initial strain rate $\dot{\epsilon}_0$, at computation end</figcaption></figure>

Note that cases A and B are exactly the same for $\ell = \ell_0/4$. Case B was added to the study because one could expect in case A some parasite effects due to mesh refinement which may bias the interpretation of the results. However, one can observe on Figure 9 that for the dissipated energy and for each strain rate, the case A and B give similar results. Higher differences are observed for the mean fragment size, but the differences lies within the stochastic spread of the results. It seems that the dissipated energy is less sensitive to $\ell$ for the lowest strain rates.

Finally, regarding computing times, for $h_e = \ell/5$ (for a total of 4525 finite elements), the computation of $d^{k+1}$ (for one time step and one iteration) takes less than one second when using the projectors $\pi^l \bar{d}(x)$ and $\pi^u \bar{d}(x)$ of equations (24) to decrease the size of the domain where problem (21b) is solved. On the contrary, the resolution of (21b) over the entire domain, using a trust-region algorithm implemented in the `scipyedenminimize` Python routine (which is more adapted to large scale problems), takes around 350 seconds, which shows the interest of using $\pi^l \bar{d}(x)$ and $\pi^u \bar{d}(x)$.

7 Discussion

In this section, the equivalence between CZM and Lip-field is discussed. As explained in Section 4, the parameters of the Lip-field model have been derived in [1] to mimic the behavior of a CZM. To explain some of the results obtained
<figure><img src="image_9.png" /><figcaption>Figure 8: Snapshots of the damage field for the Lip-field in explicit dynamics, for $h_e = \ell/10$, for different values of the initial strain rate $\dot{\epsilon}_0$, at computation end</figcaption></figure>

Note that cases A and B are exactly the same for $\ell = \ell_0/4$. Case B was added to the study because one could expect in case A some parasite effects due to mesh refinement which may bias the interpretation of the results. However, one can observe on Figure 9 that for the dissipated energy and for each strain rate, the case A and B give similar results. Higher differences are observed for the mean fragment size, but the differences lies within the stochastic spread of the results. It seems that the dissipated energy is less sensitive to $\ell$ for the lowest strain rates.

Finally, regarding computing times, for $h_e = \ell/5$ (for a total of 4525 finite elements), the computation of $d^{k+1}$ (for one time step and one iteration) takes less than one second when using the projectors $\pi^\ell \bar{d}(x)$ and $\pi^\mathrm{u} \bar{d}(x)$ of equations (24) to decrease the size of the domain where problem (21b) is solved. On the contrary, the resolution of (21b) over the entire domain, using a trust-region algorithm implemented in the `scipy{optimize.minimize} Python routine (which is more adapted to large scale problems), takes around 350 seconds, which shows the interest of using $\pi^\ell \bar{d}(x)$ and $\pi^\mathrm{u} \bar{d}(x)$.

# 7 Discussion

In this section, the equivalence between CZM and Lip-field is discussed. As explained in Section 4, the parameters of the Lip-field model have been derived in [1] to mimic the behavior of a CZM. To explain some of the results obtained
<figure><img src="image_9.png" /><figcaption>Figure 9: Explicit Lip-field $\ell$-dependency for $\dot{\epsilon}_0 = 1 \cdot 10^4 \text{s}^{-1}$ and $\dot{\epsilon}_0 = 7.5 \cdot 10^6 \text{s}^{-1}$. Dissipated energy is plotted on the left and mean fragment size on the right. Note that the order or magnitude of the error bars on the dissipated energy plot are too small to be visible.</figcaption></figure>

in Section 6, it is necessary to have a closer look at how this equivalence was established. The idea, illustrated on Figure 10 (a) is to consider a 1D bar with a crack, modelled with a CZM on one hand, and with Lip-field on the other hand. The parameters of the Lip-field model are computed so that the displacement at the ends of the bar and the dissipated energy is the same as with the CZM.

However, this equivalence was established in the case where the Lipschitz constraint is active, that is the inequality in equation (20) is an equality (meaning that the slope of the damage field is exactly $1/\ell$). This is a necessary condition for an analytical solution to exist. Yet, zones where the slope of the damage field is strictly smaller than $1/\ell$ may appear in the general case as illustrated on Figure 10 (b). Dissipated energies between CZM and Lip-field are not the same anymore. This is especially true in fragmentation as we can see on Figure 5; the zones where the Lipschitz constraint is active, indicated by black thick lines, are very small compared to the total size of the structure.

Another source of differences between CZM and Lip-field is the case of multiple cracks, as illustrated on Figures 10 (c) and 10 (d). When two cracks are separated by a distance greater than $2\ell$ (and assuming that the Lipschitz constraint is always active), the equivalence between Lip-field and CZM is still true. However, when this distance is smaller than $2\ell$, the two Lip-field damaged zones overlap and the resulting dissipated energy is smaller than in the case where the cracks are far from each other, as illustrated on Figure 10 (d). It may explain why the computed dissipated energies seem less sensitive to $\ell$ for the lowest strain rate on Figure 9; the smaller the strain rate, the smaller the number of cracks, and therefore the smaller is the probability to get cracks emerging with a separation smaller than $\ell$. An obvious way to counter this effect would be to use a smaller value of $\ell$ for the highest strain rates, but it would require using smaller finite elements, and therefore further increasing computation times.

This section emphasizes that the equivalence between CZM and Lip-field may not be fully satisfied. In some cases it may lead to energies dissipated by Lip-field greater than the ones dissipated by CZM (situation from Figure 10 (b)) or smaller (situation from Figure 10 (d)). Both situations may frequently happen in fragmentation, which makes difficult for example to anticipate if Lip-field fragmentation computations will dissipate more or less energy than CZM fragmentation computations. However, it gives an interpretation to the differences observed between CZM and Lip-field computations.

# 8 Conclusion and future works

This paper presents the first dynamic application of the Lip-field regularization introduced in [1], to the case of 1D dynamic fragmentation. This problem is formulated as an optimization problem, from which we can obtain the classical explicit and implicit dynamic time discretizations. The optimization framework allows to look for the damage variable in a space such that it does not suffer from spurious localization. Parameters of the model are determined to mimic a cohesive zone model (CZM) with a linear traction-separation law, which makes the numerical properties easy to calibrate, knowing the tensile stress and fracture toughness of the material.
<figure><img src="image_10.png" /><figcaption>Figure 9: Explicit Lip-field $\ell$-dependency for $\dot{\epsilon}_0 = 1 \cdot 10^4 \text{ s}^{-1}$ and $\dot{\epsilon}_0 = 7.5 \cdot 10^6 \text{ s}^{-1}$. Dissipated energy is plotted on the left and mean fragment size on the right. Note that the order or magnitude of the error bars on the dissipated energy plot are too small to be visible.</figcaption></figure>

in Section 6, it is necessary to have a closer look at how this equivalence was established. The idea, illustrated on Figure 10 (a) is to consider a 1D bar with a crack, modelled with a CZM on one hand, and with Lip-field on the other hand. The parameters of the Lip-field model are computed so that the displacement at the ends of the bar and the dissipated energy is the same as with the CZM.

However, this equivalence was established in the case where the Lipschitz constraint is active, that is the inequality in equation (20) is an equality (meaning that the slope of the damage field is exactly $1/\ell$). This is a necessary condition for an analytical solution to exist. Yet, zones where the slope of the damage field is strictly smaller than $1/\ell$ may appear in the general case as illustrated on Figure 10 (b). Dissipated energies between CZM and Lip-field are not the same anymore. This is especially true in fragmentation as we can see on Figure 5; the zones where the Lipschitz constraint is active, indicated by black thick lines, are very small compared to the total size of the structure.

Another source of differences between CZM and Lip-field is the case of multiple cracks, as illustrated on Figures 10 (c) and 10 (d). When two cracks are separated by a distance greater than $2\ell$ (and assuming that the Lipschitz constraint is always active), the equivalence between Lip-field and CZM is still true. However, when this distance is smaller than $2\ell$, the two Lip-field damaged zones overlap and the resulting dissipated energy is smaller than in the case where the cracks are far from each other, as illustrated on Figure 10 (d). It may explain why the computed dissipated energies seem less sensitive to $\ell$ for the lowest strain rate on Figure 9; the smaller the strain rate, the smaller the number of cracks, and therefore the smaller is the probability to get cracks emerging with a separation smaller than $\ell$. An obvious way to counter this effect would be to use a smaller value of $\ell$ for the highest strain rates, but it would require using smaller finite elements, and therefore further increasing computation times.

This section emphasizes that the equivalence between CZM and Lip-field may not be fully satisfied. In some cases it may lead to energies dissipated by Lip-field greater than the ones dissipated by CZM (situation from Figure 10 (b)) or smaller (situation from Figure 10 (d)). Both situations may frequently happen in fragmentation, which makes difficult for example to anticipate if Lip-field fragmentation computations will dissipate more or less energy than CZM fragmentation computations. However, it gives an interpretation to the differences observed between CZM and Lip-field computations.

# 8 Conclusion and future works

This paper presents the first dynamic application of the Lip-field regularization introduced in [1], to the case of 1D dynamic fragmentation. This problem is formulated as an optimization problem, from which we can obtain the classical explicit and implicit dynamic time discretizations. The optimization framework allows to look for the damage variable in a space such that it does not suffer from spurious localization. Parameters of the model are determined to mimic a cohesive zone model (CZM) with a linear traction-separation law, which makes the numerical properties easy to calibrate, knowing the tensile stress and fracture toughness of the material.
<figure><img src="image_10.png" /><figcaption>Figure 10: Different cases where CZM and Lip-field may or may not be equivalent. (a) One crack and damaged zone where the Lipschitz constraint is active everywhere. (b) One crack and damage zone where the Lipschitz constraint may not be active everywhere. (c) Two cracks and two damaged zones more than $2\ell$ apart. (d) Two cracks and two damaged zones less than $2\ell$ apart.</figcaption></figure>

We would like to emphasize the simplicity of the proposed approach, that is a simple elastic potential without asymmetry between tension and compression plus a softening function computed to get the proper dissipated energy. The damage variable can be computed using any optimization package by providing the proper Lipschitz constraint. No particular treatment is required to ensure the irreversibility of the damage variable, which is directly taken into account when solving the optimization problem. This is an advantage over the phase-field approach for instance, where imposing damage irreversibility is not straightforward. Also, the Lip-field approach does not introduce any particular boundary conditions, which is not the case of phase-field where the gradient of the damage field must be orthogonal to the domain boundaries. The use of bounds, which allow to reduce the size of the domain where the damage variable must be computed, decreases the computing times significantly. This feature is also interesting in a parallel context, since the regularized damage field can be computed independently on each zone where the Lipschitz bound are different. This is even more true in dynamic fragmentation where lot of cracks are expected. Regarding computing costs, this paper is mainly a first step to demonstrate the capacity of Lip-field to work in dynamic, rather than to propose an optimal implementation. That is why a proper comparison with other approaches has not been done yet, but will be the object of a forthcoming paper.

For comparison, this paper also introduces a simplified version of a linear CZM consisting in a local damage model, akin to a crack-band model. This allows to run CZM computations without the tedious necessity to implement discrete displacement jumps. Results exhibited in this paper demonstrate the regularizing effect of the Lip-field approach, with dissipated energies converging to a finite value with finite element mesh refinement. We emphasize that the simplified CZM was used in this paper as an easy-to-implement, 1D reference solution, which behaves analytically the same as a classical CZM. However, it would likely suffer in 2D and 3D from the same mesh dependency problems as CZM and crack-band problem.

The numerical results obtained for the Lip-field and the CZM over a range of strain rates are compared to several references. Along with previous studies [11, 78, 69], this work may therefore serve as a benchmark for 1D dynamic fragmentation. Both explicit and implicit dynamic simulation results lie well within the variability range of these references and are able to reproduce the same trends. The differences observed between both approaches may be explained by the fact that at each time step, for a given displacement field and damage state, explicit dynamics will give a unique solution at the next time step, while in implicit dynamics, it may be possible to find different solutions depending on how the alternate minimization is solved. Because the results of both explicit and implicit simulations converge to the reference range, it is difficult to tell which one is better in terms of accuracy. Further, both approaches seem comparable in terms of mesh convergence rate and give results close to their CZM counterparts. Understanding the differences between explicit and implicit resolution remains an open question, that will be the subject of a forthcoming paper.

As the spatial convergence is slow, we expect that there would be a benefit to introduce cohesive zones with displacement jumps prior to the damage variable in the bulk going to one, as demonstrated by [79]. For 1D problems, the advantage would be to increase the mesh convergence rate. Regarding 2D and 3D problems, placing diffuse Lip-field damage around sharp displacement jumps would (i) provide a natural propagation criterion, (ii) help in having crack paths independent from mesh orientation, and (iii) allow modeling of complex crack patterns, representing a clear advantage over classical cohesive zone models. Future work will serve to verify this expectation.
<figure><img src="image_11.png" /><figcaption>Figure 10: Different cases where CZM and Lip-field may or may not be equivalent. (a) One crack and damaged zone where the Lipschitz constraint is active everywhere. (b) One crack and damage zone where the Lipschitz constraint may not be active everywhere. (c) Two cracks and two damaged zones more than $2\ell$ apart. (d) Two cracks and two damaged zones less than $2\ell$ apart.</figcaption></figure>

We would like to emphasize the simplicity of the proposed approach, that is a simple elastic potential without asymmetry between tension and compression plus a softening function computed to get the proper dissipated energy. The damage variable can be computed using any optimization package by providing the proper Lipschitz constraint. No particular treatment is required to ensure the irreversibility of the damage variable, which is directly taken into account when solving the optimization problem. This is an advantage over the phase-field approach for instance, where imposing damage irreversibility is not straightforward. Also, the Lip-field approach does not introduce any particular boundary conditions, which is not the case of phase-field where the gradient of the damage field must be orthogonal to the domain boundaries. The use of bounds, which allow to reduce the size of the domain where the damage variable must be computed, decreases the computing times significantly. This feature is also interesting in a parallel context, since the regularized damage field can be computed independently on each zone where the Lipschitz bound are different. This is even more true in dynamic fragmentation where lot of cracks are expected. Regarding computing costs, this paper is mainly a first step to demonstrate the capacity of Lip-field to work in dynamic, rather than to propose an optimal implementation. That is why a proper comparison with other approaches has not been done yet, but will be the object of a forthcoming paper.

For comparison, this paper also introduces a simplified version of a linear CZM consisting in a local damage model, akin to a crack-band model. This allows to run CZM computations without the tedious necessity to implement discrete displacement jumps. Results exhibited in this paper demonstrate the regularizing effect of the Lip-field approach, with dissipated energies converging to a finite value with finite element mesh refinement. We emphasize that the simplified CZM was used in this paper as an easy-to-implement, 1D reference solution, which behaves analytically the same as a classical CZM. However, it would likely suffer in 2D and 3D from the same mesh dependency problems as CZM and crack-band problem.

The numerical results obtained for the Lip-field and the CZM over a range of strain rates are compared to several references. Along with previous studies [11, 78, 69], this work may therefore serve as a benchmark for 1D dynamic fragmentation. Both explicit and implicit dynamic simulation results lie well within the variability range of these references and are able to reproduce the same trends. The differences observed between both approaches may be explained by the fact that at each time step, for a given displacement field and damage state, explicit dynamics will give a unique solution at the next time step, while in implicit dynamics, it may be possible to find different solutions depending on how the alternate minimization is solved. Because the results of both explicit and implicit simulations converge to the reference range, it is difficult to tell which one is better in terms of accuracy. Further, both approaches seem comparable in terms of mesh convergence rate and give results close to their CZM counterparts. Understanding the differences between explicit and implicit resolution remains an open question, that will be the subject of a forthcoming paper.

As the spatial convergence is slow, we expect that there would be a benefit to introduce cohesive zones with displacement jumps prior to the damage variable in the bulk going to one, as demonstrated by [79]. For 1D problems, the advantage would be to increase the mesh convergence rate. Regarding 2D and 3D problems, placing diffuse Lip-field damage around sharp displacement jumps would (i) provide a natural propagation criterion, (ii) help in having crack paths independent from mesh orientation, and (iii) allow modeling of complex crack patterns, representing a clear advantage over classical cohesive zone models. Future work will serve to verify this expectation.
Finally, further studies could explore the effect of using spatially-correlated random fields for material properties in fragmentation studies. In particular, it would be meaningful to examine the relationship between the material correlation length and the failure length scale as it pertains to fragment size distributions.

## Acknowledgments

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525.

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or United States Government.

## Appendix: A local damage model equivalent to a linear cohesive zone model

In this section, we derive the expression of $h_{CZM}$ so that the behavior of a 1D element modeled with a local damage model is the same as a 1D element modeled with a linear CZM. The opening at the tip of a 1D element modeled with a CZM with a displacement jump $w$ is:

$$u = \frac{h_e \sigma}{E} + w \qquad (53)$$

and for a 1D element with a damage variable $d$:

$$u = \frac{h_e \sigma}{E g(d)} \qquad (54)$$

so the expression of $w$ is:

$$w = \frac{h_e \sigma}{E} \left( \frac{1 - g(d)}{g(d)} \right) \qquad (55)$$

By injecting the expression of $w$ above to the linear CZM relation $\sigma = f(w)$, we get:

$$\sigma = \sigma_c - \lambda_c I(d) \qquad (56)$$

where $\lambda_c = \frac{\sigma_c h_e}{E w_c}$ and $I(d) = \frac{1-g(d)}{g(d)}$, which gives the expression of the element stress:

$$\sigma = \frac{\sigma_c}{1 + \lambda_c I(d)} \qquad (57)$$

The damage initiation condition gives:

$$Y_c = - \frac{\sigma_c^2}{2E} g'(0) \qquad (58)$$

and the damage growth condition is:

$$-\frac{E g'(d)}{2} \epsilon^2 = -\frac{g'(d)}{2E} \sigma^2 = Y_c H_{CZM}(d) \qquad (59)$$

By combining the above equations, the expression of $H_{CZM}$ is:

$$H_{CZM}(d) = - \frac{g'(d)}{((1 - \lambda_c)g(d) + \lambda_c)^2} \qquad (60)$$

and the expression of $h_{CZM}$ is:

$$h_{CZM}(d) = \frac{1}{(1 - \lambda_c)} \left( \frac{1}{(1 - \lambda_c)g(d) + \lambda_c} - 1 \right) \qquad (61)$$
Finally, further studies could explore the effect of using spatially-correlated random fields for material properties in fragmentation studies. In particular, it would be meaningful to examine the relationship between the material correlation length and the failure length scale as it pertains to fragment size distributions.

## Acknowledgments

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525.

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or United States Government.

## Appendix: A local damage model equivalent to a linear cohesive zone model

In this section, we derive the expression of $h_{CZM}$ so that the behavior of a 1D element modeled with a local damage model is the same as a 1D element modeled with a linear CZM. The opening at the tip of a 1D element modeled with a CZM with a displacement jump $w$ is:

$$ u = \frac{h_e \sigma}{E} + w \qquad (53) $$

and for a 1D element with a damage variable $d$:

$$ u = \frac{h_e \sigma}{E g(d)} \qquad (54) $$

so the expression of $w$ is:

$$ w = \frac{h_e \sigma}{E} \left( \frac{1 - g(d)}{g(d)} \right) \qquad (55) $$

By injecting the expression of $w$ above to the linear CZM relation $\sigma = f(w)$, we get:

$$ \sigma = \sigma_c - \lambda_c I(d) \qquad (56) $$

where $\lambda_c = \frac{\sigma_c h_e}{E w_c}$ and $I(d) = \frac{1-g(d)}{g(d)}$, which gives the expression of the element stress:

$$ \sigma = \frac{\sigma_c}{1 + \lambda_c I(d)} \qquad (57) $$

The damage initiation condition gives:

$$ Y_c = -\frac{\sigma_c^2}{2E} g'(0) \qquad (58) $$

and the damage growth condition is:

$$ -\frac{E g'(d)}{2} \epsilon^2 = -\frac{g'(d)}{2E} \sigma^2 = Y_c H_{CZM}(d) \qquad (59) $$

By combining the above equations, the expression of $H_{CZM}$ is:

$$ H_{CZM}(d) = -\frac{g'(d)}{((1 - \lambda_c)g(d) + \lambda_c)^2} \qquad (60) $$

and the expression of $h_{CZM}$ is:

$$ h_{CZM}(d) = \frac{1}{(1 - \lambda_c)} \left( \frac{1}{(1 - \lambda_c)g(d) + \lambda_c} - 1 \right) \qquad (61) $$
References

[1] N. Moës and N. Chevaugeon. Lipschitz regularization for softening material models: the Lip-field approach. *Comptes Rendus. Mécanique*, 349(2):415-434, 2021.

[2] N. Chevaugeon and N. Moës. Lipschitz regularization for fracture: the Lip-field approach. *arXiv preprint arXiv:2111.04771*, 2021.

[3] C. Denoual and F. Hild. A damage model for the dynamic fragmentation of brittle solids. *Computer Methods in Applied Mechanics and Engineering*, 183(3-4):247-258, 2000.

[4] N.F. Mott. Fragmentation of shell cases. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 189(1018):300-308, 1947.

[5] D.E. Grady. Local inertial effects in dynamic fragmentation. *Journal of Applied Physics*, 53(1):322-325, 1982.

[6] L.A. Glenn and A. Chudnovsky. Strain-energy effects on dynamic fragmentation. *Journal of Applied Physics*, 59(4):1379-1380, 1986.

[7] W. J. Drugan. Dynamic fragmentation of brittle materials: Analytical mechanics-based models. *Journal of the Mechanics and Physics of Solids*, 49(6):1181-1208, 2001.

[8] Godofredo T. Camacho and Michael Ortiz. Computational modelling of impact damage in brittle materials. *International Journal of Solids and Structures*, 33(20-22):2899-2938, 1996.

[9] O. Miller, L.B. Freund, and A. Needleman. Modeling and simulation of dynamic fragmentation in brittle materials. *International Journal of Fracture*, 96(2):101-125, 1999.

[10] A. Pandolfi, P. Krysl, and M. Ortiz. Finite element simulation of ring expansion and fragmentation: the capturing of length and time scales through cohesive models of fracture. *International Journal of Fracture*, 95(1):279-297, 1999.

[11] F. Zhou, J.-F. Molinari, and K.T. Ramesh. Effects of material properties on the fragmentation of brittle materials. *International Journal of Fracture*, 139(2):169-196, 2006.

[12] J. Pearson. A fragmentation model for cylindrical warheads. Technical report, Naval Weapons Center, China Lake, CA, 1990.

[13] A.L. Tonge and K.T. Ramesh. Multi-scale defect interactions in high-rate failure of brittle materials, Part II: Application to design of protection materials. *Journal of the Mechanics and Physics of Solids*, 86:237-258, 2016.

[14] A.P.S. Selvadurai. Fragmentation of ice sheets during impact. *Computer Modeling in Engineering and Science*, 52:259-277, 2009.

[15] J.A. Åström and D.I. Benn. Effective rheology across the fragmentation transition for sea ice and ice shelves. *Geophysical Research Letters*, 46(22):13099-13106, 2019.

[16] J.E. Bishop, M.J. Martinez, and P. Newell. Simulating fragmentation and fluid-induced fracture in disordered media using random finite-element meshes. *International Journal for Multiscale Computational Engineering*, 14(4), 2016.

[17] T. Hu, J. Guilleminot, and J.E. Dolbow. A phase-field model of fracture with frictionless contact and random fracture properties: Application to thin-film fracture and soil desiccation. *Computer Methods in Applied Mechanics and Engineering*, 368:113106, 2020.

[18] W. Jiang, B.W. Spencer, and J.E. Dolbow. Ceramic nuclear fuel fracture modeling with the extended finite element method. *Engineering Fracture Mechanics*, 223:106713, 2020.

[19] D.S. Dugdale. Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids*, 8:100-104, 1960.

[20] G.I. Baranblatt. The mathematical theory of equilibrium cracks formed by brittle fracture. *Zh. Prkl. Mekh. Tekh. Fiz.*, 25(4):3-56, 1961.

[21] A. Hillerborg, M. Modéer, and P-E. Petersson. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research*, 6:773-782, 1976.

[22] K. Park and G.H. Paulino. Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces. *Applied Mechanics Reviews*, 64(6), 2011.

[23] Raúl Radovitzky, Andrew Seagraves, Mike Tupek, and Ludovic Noels. A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method. *Computer Methods in Applied Mechanics and Engineering*, 200(1-4):326-344, 2011.
References

[1] N. Moës and N. Chevaugeon. Lipschitz regularization for softening material models: the Lip-field approach. *Comptes Rendus. Mécanique*, 349(2):415-434, 2021.

[2] N. Chevaugeon and N. Moës. Lipschitz regularization for fracture: the Lip-field approach. *arXiv preprint arXiv:2111.04771*, 2021.

[3] C. Denoual and F. Hild. A damage model for the dynamic fragmentation of brittle solids. *Computer Methods in Applied Mechanics and Engineering*, 183(3-4):247-258, 2000.

[4] N.F. Mott. Fragmentation of shell cases. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 189(1018):300-308, 1947.

[5] D.E. Grady. Local inertial effects in dynamic fragmentation. *Journal of Applied Physics*, 53(1):322-325, 1982.

[6] L.A. Glenn and A. Chudnovsky. Strain-energy effects on dynamic fragmentation. *Journal of Applied Physics*, 59(4):1379-1380, 1986.

[7] W. J. Drugan. Dynamic fragmentation of brittle materials: Analytical mechanics-based models. *Journal of the Mechanics and Physics of Solids*, 49(6):1181-1208, 2001.

[8] Godofredo T. Camacho and Michael Ortiz. Computational modelling of impact damage in brittle materials. *International Journal of Solids and Structures*, 33(20-22):2899-2938, 1996.

[9] O. Miller, L.B. Freund, and A. Needleman. Modeling and simulation of dynamic fragmentation in brittle materials. *International Journal of Fracture*, 96(2):101-125, 1999.

[10] A. Pandolfi, P. Krysl, and M. Ortiz. Finite element simulation of ring expansion and fragmentation: the capturing of length and time scales through cohesive models of fracture. *International Journal of Fracture*, 95(1):279-297, 1999.

[11] F. Zhou, J.-F. Molinari, and K.T. Ramesh. Effects of material properties on the fragmentation of brittle materials. *International Journal of Fracture*, 139(2):169-196, 2006.

[12] J. Pearson. A fragmentation model for cylindrical warheads. Technical report, Naval Weapons Center, China Lake, CA, 1990.

[13] A.L. Tonge and K.T. Ramesh. Multi-scale defect interactions in high-rate failure of brittle materials, Part II: Application to design of protection materials. *Journal of the Mechanics and Physics of Solids*, 86:237-258, 2016.

[14] A.P.S. Selvadurai. Fragmentation of ice sheets during impact. *Computer Modeling in Engineering and Science*, 52:259-277, 2009.

[15] J.A. Åström and D.I. Benn. Effective rheology across the fragmentation transition for sea ice and ice shelves. *Geophysical Research Letters*, 46(22):13099-13106, 2019.

[16] J.E. Bishop, M.J. Martinez, and P. Newell. Simulating fragmentation and fluid-induced fracture in disordered media using random finite-element meshes. *International Journal for Multiscale Computational Engineering*, 14(4), 2016.

[17] T. Hu, J. Guilleminot, and J.E. Dolbow. A phase-field model of fracture with frictionless contact and random fracture properties: Application to thin-film fracture and soil desiccation. *Computer Methods in Applied Mechanics and Engineering*, 368:113106, 2020.

[18] W. Jiang, B.W. Spencer, and J.E. Dolbow. Ceramic nuclear fuel fracture modeling with the extended finite element method. *Engineering Fracture Mechanics*, 223:106713, 2020.

[19] D.S. Dugdale. Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids*, 8:100-104, 1960.

[20] G.I. Baranblatt. The mathematical theory of equilibrium cracks formed by brittle fracture. *Zh. Prkl. Mekh. Tekh. Fiz.*, 25(4):3-56, 1961.

[21] A. Hillerborg, M. Modéer, and P-E. Petersson. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research*, 6:773-782, 1976.

[22] K. Park and G.H. Paulino. Cohesive zone models: A critical review of traction-separation relationships across fracture surfaces. *Applied Mechanics Reviews*, 64(6), 2011.

[23] Raúl Radovitzky, Andrew Seagraves, Mike Tupek, and Ludovic Noels. A scalable 3D fracture and fragmentation algorithm based on a hybrid, discontinuous Galerkin, cohesive element method. *Computer Methods in Applied Mechanics and Engineering*, 200(1-4):326-344, 2011.
[24] N. Moës and T. Belytschko. Extended finite element method for cohesive crack growth. *Engineering Fracture Mechanics*, 69(7):813-833, May 2002.

[25] Timon Rabczuk, Josef Eibl, and Lothar Stempniewski. Numerical analysis of high speed concrete fragmentation using a meshfree lagrangian method. *Engineering Fracture Mechanics*, 71(4-6):547-556, 2004.

[26] R. Pramanik and D. Deb. Implementation of smoothed particle hydrodynamics for detonation of explosive with application to rock fragmentation. *Rock Mechanics and Rock Engineering*, 48(4):1683-1698, 2015.

[27] B. Li, A. Pandolfi, and M. Ortiz. Material-point erosion simulation of dynamic fragmentation of metals. *Mechanics of Materials*, 80:288-297, 2015.

[28] F. Kun and H.J. Herrmann. A study of fragmentation processes using a discrete element method. *Computer Methods in Applied Mechanics and Engineering*, 138(1-4):3-18, 1996.

[29] Pai-Chen Guan, Sheng-Wei Chi, Jiun-Shyan Chen, Thomas R Slawson, and Michael J Roth. Semi-lagrangian reproducing kernel particle method for fragment-impact problems. *International Journal of Impact Engineering*, 38(12):1033-1047, 2011.

[30] X. Lai, B. Ren, H. Fan, S. Li, C.T. Wu, R.A. Regueiro, and L. Liu. Peridynamics simulations of geomaterial fragmentation by impulse loads. *International Journal for Numerical and Analytical Methods in Geomechanics*, 39(12):1304-1330, 2015.

[31] Huajian Gao and Patrick Klein. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bonds. *Journal of the Mechanics and Physics of Solids*, 46(2):187-218, 1998.

[32] PA Klein, JW Foulk, EP Chen, SA Wimmer, and HJ Gao. Physics-based modeling of brittle fracture: cohesive formulations and the application of meshfree methods. *Theoretical and Applied Fracture Mechanics*, 37(1-3):99-166, 2001.

[33] Xiaodan Ren, Jiun-Shyan Chen, Jie Li, TR Slawson, and MJ Roth. Micro-cracks informed damage models for brittle solids. *International Journal of Solids and Structures*, 48(10):1560-1571, 2011.

[34] Timon Rabczuk and Ted Belytschko. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. *International Journal for Numerical Methods in Engineering*, 61(13):2316-2343, 2004.

[35] Xiaohu Liu, Shaofan Li, and Ni Sheng. A cohesive finite element for quasi-continua. *Computational Mechanics*, 42(4):543-553, 2008.

[36] John T Foster, Stewart A Silling, and Weinong Chen. An energy based failure criterion for use with peridynamic states. *International Journal for Multiscale Computational Engineering*, 9(6), 2011.

[37] Minghua He and Shaofan Li. An embedded atom hyperelastic constitutive model and multiscale cohesive finite element method. *Computational Mechanics*, 49(3):337-355, 2012.

[38] L. Kachanov. Rupture Time Under Creep Conditions. *Otdelenie tekhnicheskich nauk*, pages 26-31, 1958.

[39] J.-L. Chaboche and J. Lemaitre. *Mécanique des matériaux solides*. Paris, Bordas, Paris, d'unod edition, 1988.

[40] Z.P. Bažant. Instability, Ductility and Size-effect in Strain-Softening Concrete. *Journal Of The Engineering Mechanics Division*, 102(2):331-343, 1976.

[41] E. Cosserat and F. Cosserat. *Théorie des corps déformables*. Herman, Paris, 1909.

[42] R. Chambon, D. Caillerie, and N. El Hassan. One-dimensional localisation studied with a second grade model. *European Journal of Mechanics - A/Solids*, 17(4):637-656, 1998.

[43] G. Pijaudier-Cabot and Z.P. Bažant. Non local damage theory. *Journal of Engineering Mechanics*, 113(10):1512-1533, 1987.

[44] D. Lasry and T. Belytschko. Localization limiters in transient problems. *International Journal of Solids and Structures*, 24(6):581-597, 1988.

[45] R.H.J. Peerlings, R. de Borst, W. A. M. Brekelmans, and J. H. P. De Vree. Gradient Enhanced Damage for Quasi-Brittle Materials. *International Journal for Numerical Methods in Engineering*, 39(19):3391-3403, oct 1996.

[46] M. Jirásek. Nonlocal models for damage and fracture: comparison of approaches. *International Journal of Solids and Structures*, 35(31-32):4133-4145, 1998.

[47] C. Giry, F. Dufour, and J. Mazars. Stress-based nonlocal damage model. *International Journal of Solids and Structures*, 48:3431-3443, 2011.

[48] A. Needleman. Material rate dependence and mesh sensitivity in localization problems. *Computer Methods in Applied Mechanics and Engineering*, 67(1):69-85, 1988.
[24] N. Moës and T. Belytschko. Extended finite element method for cohesive crack growth. *Engineering Fracture Mechanics*, 69(7):813–833, May 2002.

[25] Timon Rabczuk, Josef Eibl, and Lothar Stempniewski. Numerical analysis of high speed concrete fragmentation using a meshfree lagrangian method. *Engineering Fracture Mechanics*, 71(4-6):547–556, 2004.

[26] R. Pramanik and D. Deb. Implementation of smoothed particle hydrodynamics for detonation of explosive with application to rock fragmentation. *Rock Mechanics and Rock Engineering*, 48(4):1683–1698, 2015.

[27] B. Li, A. Pandolfi, and M. Ortiz. Material-point erosion simulation of dynamic fragmentation of metals. *Mechanics of Materials*, 80:288–297, 2015.

[28] F. Kun and H.J. Herrmann. A study of fragmentation processes using a discrete element method. *Computer Methods in Applied Mechanics and Engineering*, 138(1-4):3–18, 1996.

[29] Pai-Chen Guan, Sheng-Wei Chi, Jiun-Shyan Chen, Thomas R Slawson, and Michael J Roth. Semi-lagrangian reproducing kernel particle method for fragment-impact problems. *International Journal of Impact Engineering*, 38(12):1033–1047, 2011.

[30] X. Lai, B. Ren, H. Fan, S. Li, C.T. Wu, R.A. Regueiro, and L. Liu. Peridynamics simulations of geomaterial fragmentation by impulse loads. *International Journal for Numerical and Analytical Methods in Geomechanics*, 39(12):1304–1330, 2015.

[31] Huajian Gao and Patrick Klein. Numerical simulation of crack growth in an isotropic solid with randomized internal cohesive bonds. *Journal of the Mechanics and Physics of Solids*, 46(2):187–218, 1998.

[32] PA Klein, JW Foulk, EP Chen, SA Wimmer, and HJ Gao. Physics-based modeling of brittle fracture: cohesive formulations and the application of meshfree methods. *Theoretical and Applied Fracture Mechanics*, 37(1-3):99–166, 2001.

[33] Xiaodan Ren, Jiun-Shyan Chen, Jie Li, TR Slawson, and MJ Roth. Micro-cracks informed damage models for brittle solids. *International Journal of Solids and Structures*, 48(10):1560–1571, 2011.

[34] Timon Rabczuk and Ted Belytschko. Cracking particles: a simplified meshfree method for arbitrary evolving cracks. *International Journal for Numerical Methods in Engineering*, 61(13):2316–2343, 2004.

[35] Xiaohu Liu, Shaofan Li, and Ni Sheng. A cohesive finite element for quasi-continua. *Computational Mechanics*, 42(4):543–553, 2008.

[36] John T Foster, Stewart A Silling, and Weinong Chen. An energy based failure criterion for use with peridynamic states. *International Journal for Multiscale Computational Engineering*, 9(6), 2011.

[37] Minghua He and Shaofan Li. An embedded atom hyperelastic constitutive model and multiscale cohesive finite element method. *Computational Mechanics*, 49(3):337–355, 2012.

[38] L. Kachanov. Rupture Time Under Creep Conditions. *Otdelenie tekhnicheskich nauk*, pages 26–31, 1958.

[39] J.-L. Chaboche and J. Lemaitre. *Mécanique des matériaux solides*. Paris, Bordas, Paris, d'unod edition, 1988.

[40] Z.P. Bažant. Instability, Ductility and Size-effect in Strain-Softening Concrete. *Journal Of The Engineering Mechanics Division*, 102(2):331–343, 1976.

[41] E. Cosserat and F. Cosserat. *Théorie des corps déformables*. Herman, Paris, 1909.

[42] R. Chambon, D. Caillerie, and N. El Hassan. One-dimensional localisation studied with a second grade model. *European Journal of Mechanics - A/Solids*, 17(4):637–656, 1998.

[43] G. Pijaudier-Cabot and Z.P. Bažant. Non local damage theory. *Journal of Engineering Mechanics*, 113(10):1512–1533, 1987.

[44] D. Lasry and T. Belytschko. Localization limiters in transient problems. *International Journal of Solids and Structures*, 24(6):581–597, 1988.

[45] R.H.J. Peerlings, R. de Borst, W. A. M. Brekelmans, and J. H. P. De Vree. Gradient Enhanced Damage for Quasi-Brittle Materials. *International Journal for Numerical Methods in Engineering*, 39(19):3391–3403, oct 1996.

[46] M. Jirásek. Nonlocal models for damage and fracture: comparison of approaches. *International Journal of Solids and Structures*, 35(31-32):4133–4145, 1998.

[47] C. Giry, F. Dufour, and J. Mazars. Stress-based nonlocal damage model. *International Journal of Solids and Structures*, 48:3431–3443, 2011.

[48] A. Needleman. Material rate dependence and mesh sensitivity in localization problems. *Computer Methods in Applied Mechanics and Engineering*, 67(1):69–85, 1988.
[49] J.F. Dubé, G. Pijaudier-Cabot, and C. La Borderie. Rate dependent damage model for concrete in dynamics. *Journal of Engineering Mechanics*, 122(10):939-947, 1996.

[50] O. Allix and J.F Deü. Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading. *Engineering transactions*, 45(1):29-46, 1997.

[51] A. Suffis, T.A.A. Lubrecht, and A. Combescure. Damage model with delay effect analytical and numerical studies of the evolution of the characteristic damage length. *International Journal of Solids and Structures*, 40(13-14):3463-3476, 2003.

[52] R. Desmorat, M. Chambart, F. Gatuingt, and D. Guilbaud. Delay-active damage versus non-local enhance-ment for anisotropic damage dynamics computations with alternated loading. *Engineering Fracture Mechanics*, 77(12):2294-2315, 2010.

[53] J. Zghal and N. Moës. Analysis of the delayed damage model for three one-dimensional loading scenarii. *Comptes Rendus. Physique*, 21(6):527-537, 2020.

[54] A. Karma, D. Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Physical Review Letters*, 87(4):045501, 2001.

[55] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273-1311, 2010.

[56] G. A. Francfort and J-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319-1342, 1998.

[57] C.V. Verhoosel and R. de Borst. A phase-field model for cohesive fracture. *International Journal for Numerical Methods in Engineering*, 96(1):43-62, 2013.

[58] J-Y. Wu. A geometrically regularized gradient-damage model with energetic equivalence. *Computer Methods in Applied Mechanics and Engineering*, 328:612-637, 2018.

[59] B. Talamini, M.R. Tupek, A.J. Stershic, T. Hu, J.W. Foulk III, J.T. Ostien, and J.E. Dolbow. Attaining regulariza-tion length insensitivity in phase-field models of ductile failure. *Computer Methods in Applied Mechanics and Engineering*, 384:113936, 2021.

[60] M. J. Borden, C.V. Verhoosel, M. A. Scott, T.J.R. Hughes, and C.M. Landis. A phase-field description of dynamic brittle fracture. *Computer Methods in Applied Mechanics and Engineering*, 217:77-95, 2012.

[61] H. L. Ren, X. Y. Zhuang, C. Anitescu, and T. Rabczuk. An explicit phase field method for brittle dynamic fracture. *Computers and Structures*, 217:45-56, 2019.

[62] R.J.M. Geelen, Y. Liu, T. Hu, M.R. Tupek, and J.E. Dolbow. A phase-field formulation for dynamic cohesive fracture. *Computer Methods in Applied Mechanics and Engineering*, 348:680-711, 2019.

[63] A.G. Fischer and J.-J. Marigo. Gradient damage models applied to dynamic fragmentation of brittle materials. *International Journal of Fracture*, 220(2):143-165, 2019.

[64] N. Moës, C. Stolz, P.E. Bernard, and N. Chevaugeon. A level set based model for damage growth : the thick level set approach. *International Journal for Numerical Methods in Engineering*, 86(3):358-380, 2011.

[65] C. Stolz and N. Moës. A new model of damage: a moving thick layer approach. *International Journal of Fracture*, 174:49-60, 2012.

[66] P.E. Bernard, N. Moës, and N. Chevaugeon. Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings. *Computer Methods in Applied Mechanics and Engineering*, 233-236:11-27, August 2012.

[67] A. Salzman, N. Moës, and N. Chevaugeon. On use of the Thick Level Set method in 3D quasi-static crack simulation of quasi-brittle material. *International Journal of Fracture*, 202(1):21-49, 2016.

[68] K. Moreau, N. Moës, D. Picart, and L. Stainier. Explicit dynamics with a non-local damage model using the thick level set approach. *International Journal for Numerical Methods in Engineering*, 102(3-4):808-838, 2015.

[69] A. J. Stershic, J.E. Dolbow, and N. Moës. The Thick Level-Set model for dynamic fragmentation. *Engineering Fracture Mechanics*, 172:39-60, 2017.

[70] Michel Frémond and Claude Stolz. On alternative approaches for graded damage modelling. In *Models, Simulation, and Experimental Issues in Structural Mechanics*, pages 87-104. 2017.

[71] Nunziante Valoroso and Claude Stolz. Graded damage in quasi-brittle solids. *International Journal for Numerical Methods in Engineering*, (December 2021):1-32, 2022.
[49] J.F. Dubé, G. Pijaudier-Cabot, and C. La Borderie. Rate dependent damage model for concrete in dynamics. *Journal of Engineering Mechanics*, 122(10):939-947, 1996.

[50] O. Allix and J.F Deü. Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading. *Engineering transactions*, 45(1):29-46, 1997.

[51] A. Suffis, T.A.A. Lubrecht, and A. Combescure. Damage model with delay effect analytical and numerical studies of the evolution of the characteristic damage length. *International Journal of Solids and Structures*, 40(13-14):3463-3476, 2003.

[52] R. Desmorat, M. Chambart, F. Gatuingt, and D. Guilbaud. Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading. *Engineering Fracture Mechanics*, 77(12):2294-2315, 2010.

[53] J. Zghal and N. Moës. Analysis of the delayed damage model for three one-dimensional loading scenarii. *Comptes Rendus. Physique*, 21(6):527-537, 2020.

[54] A. Karma, D. Kessler, and H. Levine. Phase-Field Model of Mode III Dynamic Fracture. *Physical Review Letters*, 87(4):045501, 2001.

[55] C. Miehe, F. Welschinger, and M. Hofacker. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273-1311, 2010.

[56] G. A. Francfort and J-J. Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319-1342, 1998.

[57] C.V. Verhoosel and R. de Borst. A phase-field model for cohesive fracture. *International Journal for Numerical Methods in Engineering*, 96(1):43-62, 2013.

[58] J-Y. Wu. A geometrically regularized gradient-damage model with energetic equivalence. *Computer Methods in Applied Mechanics and Engineering*, 328:612-637, 2018.

[59] B. Talamini, M.R. Tupek, A.J. Stershic, T. Hu, J.W. Foulk III, J.T. Ostien, and J.E. Dolbow. Attaining regularization length insensitivity in phase-field models of ductile failure. *Computer Methods in Applied Mechanics and Engineering*, 384:113936, 2021.

[60] M. J. Borden, C.V. Verhoosel, M. A. Scott, T.J.R. Hughes, and C.M. Landis. A phase-field description of dynamic brittle fracture. *Computer Methods in Applied Mechanics and Engineering*, 217:77-95, 2012.

[61] H. L. Ren, X. Y. Zhuang, C. Anitescu, and T. Rabczuk. An explicit phase field method for brittle dynamic fracture. *Computers and Structures*, 217:45-56, 2019.

[62] R.J.M. Geelen, Y. Liu, T. Hu, M.R. Tupek, and J.E. Dolbow. A phase-field formulation for dynamic cohesive fracture. *Computer Methods in Applied Mechanics and Engineering*, 348:680-711, 2019.

[63] A.G. Fischer and J.-J. Marigo. Gradient damage models applied to dynamic fragmentation of brittle materials. *International Journal of Fracture*, 220(2):143-165, 2019.

[64] N. Moës, C. Stolz, P.E. Bernard, and N. Chevaugeon. A level set based model for damage growth : the thick level set approach. *International Journal for Numerical Methods in Engineering*, 86(3):358-380, 2011.

[65] C. Stolz and N. Moës. A new model of damage: a moving thick layer approach. *International Journal of Fracture*, 174:49-60, 2012.

[66] P.E. Bernard, N. Moës, and N. Chevaugeon. Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings. *Computer Methods in Applied Mechanics and Engineering*, 233-236:11-27, August 2012.

[67] A. Salzman, N. Moës, and N. Chevaugeon. On use of the Thick Level Set method in 3D quasi-static crack simulation of quasi-brittle material. *International Journal of Fracture*, 202(1):21-49, 2016.

[68] K. Moreau, N. Moës, D. Picart, and L. Stainier. Explicit dynamics with a non-local damage model using the thick level set approach. *International Journal for Numerical Methods in Engineering*, 102(3-4):808-838, 2015.

[69] A. J. Stershic, J.E. Dolbow, and N. Moës. The Thick Level-Set model for dynamic fragmentation. *Engineering Fracture Mechanics*, 172:39-60, 2017.

[70] Michel Frémond and Claude Stolz. On alternative approaches for graded damage modelling. In *Models, Simulation, and Experimental Issues in Structural Mechanics*, pages 87-104. 2017.

[71] Nunziante Valoroso and Claude Stolz. Graded damage in quasi-brittle solids. *International Journal for Numerical Methods in Engineering*, (December 2021):1-32, 2022.
[72] Joakim Petersson and Ole Sigmund. Slope constrained topology optimization. *International Journal for Numerical Methods in Engineering*, 41(8):1417–1434, 1998.

[73] D. E. Grady and M. L. Olsen. A statistics and energy based theory of dynamic fragmentation. *International Journal of Impact Engineering*, 29(1-10):293–306, 2003.

[74] H. Zhang and K. Ravi-Chandar. *On the dynamics of necking and fragmentation - I. Real-time and post-mortem observations in Al 6061-O*, volume 142. 2006.

[75] A. Parrilla Gómez, N. Moës, and C. Stolz. Comparison between Thick Level Set (TLS) and cohesive zone models. *Advanced Modeling and Simulation in Engineering Sciences*, 2(1):1–22, 2015.

[76] Z.P. Bažant. Crack Band Model for Fracture of Geomaterials. In *Proceedings of the Fourth International Conference on Numerical Methods in Geomechanics*, volume III, pages 1137–1152, 1982.

[77] Z.P. Bažant and B.H. Oh. Crack band theory for fracture of concrete. *Materials and Structures*, 16:155–177, 1983.

[78] J.-F. Molinari, G. Gazonas, R. Raghupathy, A. Rusinek, and F. Zhou. The cohesive element approach to dynamic fragmentation: the question of energy convergence. *International Journal for Numerical Methods in Engineering*, 69(3):484–503, 2007.

[79] B. Lé, N. Moës, and G. Legrain. Coupling damage and cohesive zone models with the Thick Level Set approach to fracture. *Engineering Fracture Mechanics*, 193:214–247, 2018.
[72] Joakim Petersson and Ole Sigmund. Slope constrained topology optimization. *International Journal for Numerical Methods in Engineering*, 41(8):1417–1434, 1998.

[73] D. E. Grady and M. L. Olsen. A statistics and energy based theory of dynamic fragmentation. *International Journal of Impact Engineering*, 29(1-10):293–306, 2003.

[74] H. Zhang and K. Ravi-Chandar. *On the dynamics of necking and fragmentation - I. Real-time and post-mortem observations in Al 6061-O*, volume 142. 2006.

[75] A. Parrilla Gómez, N. Moës, and C. Stolz. Comparison between Thick Level Set (TLS) and cohesive zone models. *Advanced Modeling and Simulation in Engineering Sciences*, 2(1):1–22, 2015.

[76] Z.P. Bažant. Crack Band Model for Fracture of Geomaterials. In *Proceedings of the Fourth International Conference on Numerical Methods in Geomechanics*, volume III, pages 1137–1152, 1982.

[77] Z.P. Bažant and B.H. Oh. Crack band theory for fracture of concrete. *Materials and Structures*, 16:155–177, 1983.

[78] J.-F. Molinari, G. Gazonas, R. Raghupathy, A. Rusinek, and F. Zhou. The cohesive element approach to dynamic fragmentation: the question of energy convergence. *International Journal for Numerical Methods in Engineering*, 69(3):484–503, 2007.

[79] B. Lé, N. Moës, and G. Legrain. Coupling damage and cohesive zone models with the Thick Level Set approach to fracture. *Engineering Fracture Mechanics*, 193:214–247, 2018.