# Dynamic damage modelling through phase-field approaches: Assessment, critical analysis and comparison

Lamia Mersel, Pascal Bouda, Jérémy Germain, Julien Réthoré

► To cite this version:

Lamia Mersel, Pascal Bouda, Jérémy Germain, Julien Réthoré. Dynamic damage modelling through phase-field approaches: Assessment, critical analysis and comparison. Comptes Rendus. Mécanique, 2025, 353 (G1), pp.687-724. 10.5802/crmeca.297. hal-05076962

HAL Id: hal-05076962

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

Submitted on 21 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Research article / *Article de recherche*

# Dynamic damage modelling through phase-field approaches: Assessment, critical analysis and comparison

*Modélisation dynamique de l'endommagement par l'approche champ de phase : Evaluation, analyse critique et comparaison*

Lamia Mersel<sup>®,*,a,b</sup>, Pascal Bouda<sup>®,*,c</sup>, Jérémy Germain<sup>®,*,a</sup> and Julien Réthoré<sup>®,*,b</sup>

<sup>a</sup> DMAS, ONERA, Lille, 59000, France

<sup>b</sup> Nantes Université, Ecole Centrale Nantes, CNRS, GeM, UMR 6183, 1 rue de la Noë, Nantes, 44321, France

<sup>c</sup> Université Paris-Saclay, CEA, Service d'Études Mécaniques et Thermiques, 91191 Gif-sur-Yvette, France

*E-mails*: lamia.mersel@onera.fr (L. Mersel), pascal.bouda@cea.fr (P. Bouda), jeremy.germain@onera.fr (J. Germain), julien.rethore@ec-nantes.fr (J. Réthoré)

**Abstract.** This work assesses the ability of various phase field approaches to model dynamic damage evolution up to fracture. The damage evolution can be formulated using several types of Partial Differential Equations (PDEs), such as elliptic, parabolic and hyperbolic damage PDEs. The latter enable the use of explicit time integration, which is a crucial feature if an implementation in solvers dedicated to transient analysis is targeted. The challenge of the damage irreversibility constraint through different strategies is also investigated. It is shown that the classical history variable does not strictly assure the damage irreversibility and can induce some numerical instabilities. However, combining this history variable with the damage viscosity term of the PDE might help to restore the irreversibility, even though explicit time stepping offers the opportunity to force this condition algorithmically. Finally, we propose a quantitative comparison of the solutions of each phase field dynamic approach through physically-based metrics. Eventually, the classical Kalthoff experiment is simulated and the predictions of different PDEs with different set of damage parameters are compared.

**Résumé.** Ces travaux évaluent l'aptitude de différents modèles à champs de phase pour modéliser l'évolution dynamique de l'endommagement jusqu'à rupture. L'évolution de l'endommagement peut être décrite par une équation à dérivée partielle (EDP) de type elliptique, parabolique ou hyperbolique. Ces dernières permettent l'utilisation de schémas d'intégration explicites, stratégie majoritairement employée dans les cas d'étude de dynamique transitoire. La condition d'irréversibilité de l'endommagement a également été imposée de différente manières. Il est montré que l'utilisation d'une variable d'histoire de la densité d'énergie élastique positive ne permet pas de garantir, strictement l'irréversibilité et peut induire des instabilités. Néanmoins, si cette stratégie est combinée à l'utilisation d'un terme de viscosité d'endommagement dans l'EPD gouvernant son évolution, il est possible d'assurer l'irréversibilité. Autrement, le cadre explicite de la stratégie de résolution offre la possibilité d'imposer l'irréversibilité de façon algorithmique. Ensuite,

<footnote>* Corresponding authors</footnote>
au travers des cas test d'une barre en traction et de l'expérience d'impact de Kalthoff et al., des études paramétriques ont été effectuées pour évaluer la sensibilité aux paramètres des EPDs. Pour finir, une étude comparative des solutions obtenues par les modèles de champs de phase en dynamique a été réalisée via l'utilisation de métrique basée sur des quantités physiques.

**Keywords.** Dynamic fracture, Damage phase field formulation, Explicit time integration, Sensitivity analysis.

**Mots-clés.** Rupture dynamique, formulation d'endommagement, intégration temporelle explicite, analyse de sensibilité.

**Electronic supplementary material.** Supplementary material for this article is supplied as a separate archive available from the journal's website under article's URL or from the author.

This article is a draft (not yet accepted)

# 1. Introduction

Phase field methods for fracture modelling are subjected to continuous improvements since the end of 20th century, especially in quasi-static framework. Among them is the well-known variational formulation for Griffith's theory for brittle materials. The latter is based on the regularization of the Mumford-Shah functional [1, 2], which is also known as the Ambrosio-Tortorelli approximation [3]. This approximated functional, leading to the resolution of elliptic Partial Differential Equations (PDEs), overcomes the loss of well-posedness of the boundary value problem to solve for the displacement. Like non-local damage models, this approach incorporates a damage gradient and a regularization length $l_c$ that governs the formation of a continuous damage growing band. When $l_c$ tends to 0 (i.e. the thickness of the damage area/volume vanishes), the global minimizer of the approximated functional tends, by $\Gamma$-convergence, towards the global minimizer of the Griffith functional and thus makes a bridge from the damage modelling to the fracture modelling.

Subsequently, these models have been extended to transient loading [4–8]. Dynamic damage phase-field-based modelling is described owing to a set of PDEs with an hyperbolic one for the displacement evolution and an elliptic one for the damage field if the classical formulation is adopted. In the latter, dynamic evolution of damage is solely driven by the elastodynamic equation for the displacement evolution. This strategy - which relies on an elliptic formulation for damage - is the most used in the literature and shows good agreement with the experimental crack patterns [5, 9, 10]. Generally speaking, an elliptic PDE can be considered as inappropriate to describe the dynamic nature of a field variable evolution since the notion of the information propagation speed does not exist and thus cannot be controlled [11]. Otherwise, considering the more general context of dynamic damage modelling, one can also jeopardize damage evolution main driver itself, and use direct time-dependent equations to govern the damage evolution. The parabolic type, akin to Ginzburg-Landau formulation [12], is one of the most used in the literature to model a dynamic damage evolution [6, 13–18]. These equations are often used to add delay effects through the use of a viscosity parameter, $\mu_d$, limiting the damage rate. However, this type of parabolic PDE typically predicts the information propagation at an infinite speed through the medium [19], as an elliptic one can be considered to do. Although this description is considered physically inappropriate in some scientific fields [19–21], this feature is usually neglected in the case of damage evolution that is considered to follow diffusive-like evolution and not wave-like transient evolution in time.

However, for heat transfer applications which involve fast-transient evolution of diffusive-like processes, a physical interpretation of the hyperbolic type equation was advocated [19]. Indeed,
the instantaneous propagation of the information intrinsic to the parabolic PDEs might fail in predicting rapid heat changes. This short process (an elapsed time close to the relaxation time of the solution) leads to a thermodynamic process far from equilibrium conditions, which is not taken into account by parabolic evolution. It was proposed a generalization of the Fourier's law to extend the model to a hyperbolic PDE to restore a limited speed propagation, similarly to a telegraph equation [22]. Since the 1990s, a less common hyperbolic PDE was also introduced for damage modelling [23–26]. In damage mechanics, the unbounded speed at which the information is travelling would then be expressed as a combination of fracture mechanics parameters and a damage inertia parameter $\rho_d$.

The selection of the aforementioned additional parameters ($\mu_d$, $\rho_d$) remains challenging. Some authors tried to develop identification methods of these additional terms. As an example, for damage parabolic PDE and in the context of a phase field damage model, a set of dimensionless parameters was identified based on the Linear Elastic Fracture Mechanics (LEFM) framework [16]. Finally, the method was able to numerically capture the oscillatory crack advance. In the hyperbolic type case, a parameter estimation strategy was proposed, based on a calibration with phenomenological observation [25]. In this work, a set of simulations was compared with phenomenological models leading to expressions of the damage inertia $\rho_d$ and the damage viscosity $\mu_d$. These formulations are controlled with numerous new coefficients and material parameters making the calibration tedious. Another approach to calibrate the hyperbolic damage PDE parameters might be to assume an over-damped system. For instance, a damage damping parameter expression depending on material properties and energetic quantities was developed in [27]. Besides, this expression can also be derived from a boundedness solution analysis of a generic hyperbolic equation [28].

From the computational point of view, the standard approach for solving the dynamic phase field involves the use an explicit time integration scheme like the central difference method for the elastodynamic displacement equation, and an iterative solver for the elliptic damage equation [7, 18]. Indeed, modelling the displacement due to elastic wave propagation motivates the use of explicit methods, since they are known to capture these phenomena properly. Although conditionally stable (this leading to the use of small time steps), the explicit time integration schemes can lead to a matrix-free formalism, avoiding an expensive linear system resolution required by implicit solvers. Thus, for transient dynamic analysis with an explicit solver for the displacement evolution, the iterative resolution of the elliptic damage problem must be solved at a very high frequency (typically the elastic wave velocity divided by the smallest element size in the mesh) which makes this strategy expensive in terms of computational resources. Therefore, a fully explicit strategy to solve the coupled system of equations is to be favoured, implying the use of parabolic or hyperbolic PDEs only.

In a nutshell, embedding a parabolic or an hyperbolic PDE to model damage evolution in a phase field approach considerably broaden its scope. However, the advantages and drawbacks of the different PDEs systems have been analyzed mostly individually so far, from both the numerical and the physical point of view. Thus, there is a dire need to put them in a common framework for a pragmatic comparative analysis and in depth physical interpretation of their predictions. However, in order to benefit for the potential of alternative formulation in terms of physical modelling, one has to ensure their implementation can be efficient and robust. This comes with the challenges of a fully explicit resolution strategy and its related numerical stability issue and damage irreversibility constraint. In this paper, these numerical issues are addressed and the model ability to replicate realistic crack evolution are discussed. Section 2 outlines the pre-
existing derivation of the coupled phase field equations to obtain parabolic and hyperbolic PDEs used in this study. Next, the implementation details for solving these different sets of PDEs are described in Section 3. A particular attention is drawn on the numerical stability and on the implementation of the damage irreversibility conditions. Comparison of damage formulations predictions are exposed in Sections 4 and 5. In the first one, a physically based criterion is developed to evaluate deviations from a standard solution, and a sensitivity study on additional parameters is conducted on a simple benchmark test case. Last, the classical Kalthoff and Winkler impact test is considered and the predictions of different damage formulations are compared.

## 2. Theory - Dynamic phase field for fracture modelling

In this section, the elliptic, parabolic and hyperbolic damage PDEs are to be derived depending on the considered energetic contributions. As a major ingredient of phase field models, the tension/compression energy split is also discussed. The discussion about damage irreversibility is left upon the next Section as different implementation can be adopted based upon different considerations.

### 2.1. General framework

Since the 1990's, several methods have emerged to derive dynamic phase field equations. One of these approaches is based on the principle of virtual power and thermodynamics principles [12, 24].

It relies on two main concepts: (a) the balance of the power of external, internal and inertial forces, (b) the first and the second law of thermodynamics, which is related to the Clausius-Duhem inequality.

Later on, the Hamiltonian Variational Principle was used [7, 15, 16, 26, 29]. This is done setting the Gateaux derivative of the action functional S:

$$S[T] = \int_{T} L^{dyn}(t, d, \underline{u}, \nabla d, \nabla \underline{u}, \dot{\underline{u}}) dt \quad (1)$$

to zero, where $L^{dyn}$ is the Lagrangian, $\underline{u}$ the displacement field, $\dot{\underline{u}}$ the velocity field, $t$ the time variable, $T$ the considered time interval, $d$ the damage variable, $\bar{d}$ its time derivative and $\nabla$ the gradient operator. In its general form, the Lagrangian $L^{dyn}$ is defined as follows:

$$L^{dyn}(t, \underline{u}, \dot{\underline{u}}, d, \dot{d}) = \mathcal{K}_{u} + \mathcal{K}_{d} - (\mathcal{E}_{el} + \mathcal{E}_{fr} - \mathcal{W}_{ext}(\underline{u})). \quad (2)$$

$\mathcal{K}_u$ and $\mathcal{K}_d$, are the kinetic energy related to the variation in $\underline{u}$ and $d$, respectively, $\mathcal{E}_{el}$ the strain energy, $\mathcal{E}_{fr}$ the fracture energy and $\mathcal{W}_{ext}$ the work of external forces to be defined.

Over $\Omega$, the considered bounded open subset of $R^{dim}$ ($dim$ being the dimension of the problem), the kinetic energies are defined as

$$\mathcal{K}_u = \int_{\Omega} \frac{\rho_u}{2} \dot{\underline{u}}^2 d\Omega, \quad (3)$$

$$\mathcal{K}_d = \int_{\Omega} \frac{\rho_d}{2} \dot{d}^2 d\Omega. \quad (4)$$

In these equations, $\rho_u$ is the mass density and $\rho_d$ is a so-called micro-damage inertia. The elastic strain energy is

$$\mathcal{E}_{el} = \int_{\Omega} \Psi \left( \underline{\varepsilon}(\underline{u}), d \right) d\Omega \quad (5)$$
with $\Psi$ the strain energy density which can be written as:

$$ \Psi(\underline{\epsilon}(\underline{u}), d) = \frac{1}{2} g(d) \underline{\epsilon} : \mathcal{C} : \underline{\epsilon}, \qquad (6) $$

$\mathcal{C}$ being the elastic stiffness tensor and $g(d) = (1-d)^2$ a degradation function. The definition of the Cauchy stress tensor is derived from this definition:

$$ \underline{\sigma} = g(d)\mathcal{C} : \underline{\epsilon}. \qquad (7) $$

The fracture energy density $\mathcal{E}_{fr}(\ell_c, d, \nabla d)$ is as follows:

$$ \mathcal{E}_{el} + \mathcal{E}_{fr} = \int_{\Omega} \Psi(\underline{\epsilon}(\underline{u}), d) d\Omega + \int_{\Omega} G_c \gamma(\ell_c, d, \nabla d) d\Omega. \qquad (8) $$

It is a regularized form of the original integral of the fracture dissipated energy used within small strain framework [1]. This form introduces a crack density function $\gamma$ which depends on an internal length $\ell_c$, the damage field $d$ and its gradient $\nabla d$. Among the different form of this crack density function, we choose

$$ \gamma(\ell_c, d, \nabla d) = \frac{1}{4c_w} \left( \frac{w(d)}{\ell_c} + \ell_c |\nabla d|^2 \right) = \frac{1}{2} \left( \frac{d^2}{\ell_c} + \ell_c |\nabla d|^2 \right) \qquad (9) $$

proposed in [30]. This form is based on the AT2 model such as $w = d^2$ associated to a normalized coefficient $c_w = 4 \int_0^1 \sqrt{w(x)} dx = 1/2$. In the work of [10], a comparative study between the AT1 model $w = d$, $c_w = 3/8$ and AT2 model highlights some difficulties of the latter in reproducing physical observations. The lack of elastic domain leads to excessive damage over the bulk and results in length scale - $\ell_c$ - sensitivity. A strategy to mitigate this effect is to introduce an energy density threshold, see for example [31]. While, the AT1 model is usually more suitable for brittle failure, it may require additional techniques for implementation, as outlined in [32]. However, the AT2 model is favored in this work because this is the most widely used in the computational fracture mechanics community and thus will be considered in the following as the reference case. Last, the work of external forces: $\underline{b}$ the body forces and $\underline{t}$ the surface traction applied over the domain boundary $\partial\Omega_n$ submitted to Neumann conditions, is defined as:

$$ \mathcal{W}_{ext} = \int_{\Omega} (\underline{b} \cdot \underline{u}) d\Omega + \int_{\partial\Omega_n} (\underline{t} \cdot \underline{n}_s \underline{u}) dS. \qquad (10) $$

These definitions are complemented with the definition of the admissible space for the unknown displacement and damage fields. The displacement field $\underline{u}$ belongs to the admissible set

$$ \bar{\Omega}_u(\underline{v}) = \{ \underline{u} : T \times \Omega \to H^1(\Omega, R^{dim}), \underline{u} = \underline{v} \text{ on } \partial_v \Omega \}. \qquad (11) $$

with $H^1(\Omega, R^{dim})$ the Sobolev space of vector fields defined in $\Omega$, the solid domain.

The admissible damage space includes an irreversibility constraint on its evolution, changing the topological structure of the set from a linear vector space to a convex cone, defined as

$$ \bar{\Omega}_{\bar{d}}(t) = \{ d : T \times \Omega \to [0, 1], \bar{d} \le d \le 1, d = \bar{d} \text{ on } \partial_{\bar{d}} \Omega \} \qquad (12) $$

with $\bar{d}$ denoting the solution $d$ before $t$.

The variational problem is written through the coupled system:

$$ \delta S(\underline{u}, d; \delta \underline{u}) = 0 \text{ in } \forall \delta \underline{u} \in \bar{\Omega}_u, \qquad (13a) $$

$$ \delta S(\underline{u}, d; \delta d) \ge 0 \text{ in } \forall \delta d \in \bar{\Omega}_d. \qquad (13b) $$
The variational form of the energies are:

$$ \int_t \delta \mathcal{H}_u (\dot{u}) dt = - \int_t \int_{\Omega} (\rho_u \ddot{u} \delta u) d\Omega dt, \delta \mathcal{H}_d (\dot{d}) dt = - \int_t \int_{\Omega} (\rho_d \ddot{d} \delta d) d\Omega dt \quad (14) $$

for the kinetic energies,

$$ \int_t \delta \mathcal{E}_{el}(\underline{u}, d) dt = \int_t \int_{\Omega} \left( \underline{\sigma} : \delta \underline{\nabla} \underline{u} + \frac{\partial \Psi}{\partial d} \delta d \right) d\Omega dt, \quad (15) $$

for the strain energy,

$$ \int_t \delta \mathcal{E}_{fr}(d) dt = \int_t \int_{\Omega} \left( G_c \frac{\partial \gamma}{\partial d} \delta d + G_c \underline{\nabla} \cdot \left( \frac{\partial \gamma}{\partial \underline{\nabla} d} \delta d \right) - G_c \underline{\nabla} \cdot \left( \frac{\partial \gamma}{\partial \overline{\nabla} d} \right) \delta d \right) d\Omega dt \quad (16) $$

for the fracture energy and

$$ \int_t \delta \mathcal{W}_{ext} dt = \int_t \int_{\Omega} (b \delta \underline{u}) d\Omega dt + \int_t \int_{\partial \Omega_n} (t \cdot \underline{n}_s \delta \underline{u}) dS dt \quad (17) $$

for the work of external forces. This enables one to write the weak Cauchy problem described by the coupled system of equations in the most general case, respectively for the displacement and damage evolution:

$$ \left\{ \int_t \int_{\Omega} \left( -\rho_u \ddot{u} + \underline{\nabla} \cdot \underline{\sigma} + \underline{b} \right) \delta \underline{u} d\Omega dt = 0 \quad \text{in } \forall \delta \underline{u} \in \bar{\Omega}_u, \quad (18a) $$

$$ \int_{\partial \Omega} -\underline{\sigma} \cdot \underline{n}_s \delta \underline{u} dS + \int_{\partial \Omega_n} t \cdot \underline{n}_s \delta \underline{u}_d S = 0 \quad \text{on } \partial \bar{\Omega}_u, \quad (18b) $$

$$ \left\{ \int_t \int_{\Omega} \left( \rho_d \ddot{d} + \frac{\partial \Psi}{\partial d} + G_c \frac{\partial \gamma}{\partial d} - G_c \underline{\nabla} \cdot \left( \frac{\partial \gamma}{\partial \underline{\nabla} d} \right) \right) \delta d \, d\Omega \, dt \ge 0 \quad \text{in } \forall \delta d \in \bar{\Omega}_d, \quad (19a) $$

$$ \int_{\partial \Omega_d} \frac{\partial \gamma}{\partial \underline{\nabla} d} \delta d \cdot \underline{n}_d \, dS \ge 0 \quad \text{on } \partial \bar{\Omega}_d, \quad (19b) $$

From this general form of the problem, one can chose to discard some energy contributions to derive the strong form of the coupled problem. This leads to different PDEs for the displacement and damage field. In the reminder, we focus on the analysis of transient mechanical problem meaning that the kinetic energy relative to the displacement is considered. We then focus on the derivation of PDEs for the damage problem arising from different modelling choices.

## 2.2. Derivation of damage PDEs

### 2.2.1. Elliptic damage formulation: rate-independent model [33]

For this formulation, the damage kinetic energy is discarded. From the crack density function introduced in Eq. (9), the final strong form of the damage motion equation is given by:

$$ 0 \geq -\frac{\partial \Psi}{\partial d} - \frac{G_c}{\ell_c} (d - 2 \ell_c^2 \Delta d). \quad (20) $$

The damage PDE is reformulated with a damage yield function $f(Y, d) = Y - \frac{G_c}{\ell_c} (d - 2 \ell_c^2 \Delta d)$ where Y represents an energetic crack driving force defined by:

$$ Y = -\frac{\partial \Psi}{\partial d} = -g'(d) \frac{\partial \Psi}{\partial g}. \quad (21) $$

Then, a local form of the consistency condition, similar to the Karush-Kuhn-Tucker conditions, can be derived from the irreversibility condition ($\dot{d} \ge 0$) and the damage criterion ($f(Y, d) \le 0$):

$$ \dot{d} \ge 0, \quad f(Y, d) \le 0, \quad \dot{f}(Y, d) = 0. \quad (22) $$
It is also known as the conditions of consistency derived from the first order stability condition and the energy balance proposed in [34].

From the numerical point of view, the damage irreversibility (cf. Section 3.3) is handled within the damage formulation considering an equality form of the evolution equation $f(Y, d) = 0$.

### 2.2.2. Parabolic damage formulation: rate-dependent model [35, 36]

As for the elliptic formulation, the damage kinetic energy is discarded but the system is made non conservative by considering a quadratic dissipative potential

$$ \Psi_{nc} = \frac{\mu_d}{2} \dot{d}^2 \qquad (23) $$

that is incorporated in the principle of least action:

$$ \int_t \delta L^{dyn} dt \geq \int_t \int_{\Omega} -\delta \Psi_{nc}(\dot{d}) d\Omega dt = \int_t \int_{\Omega} -\frac{\partial \Psi_{nc}}{\partial \dot{d}} \delta d\Omega dt. \qquad (24) $$

This leads to the local damage rate-dependent evolution equation:

$$ \mu_d \dot{d} \geq -\frac{\partial \Psi}{\partial d} - \frac{G_c}{\ell_c} (d - 2 \ell_c^2 \Delta d) = f(Y, d). \qquad (25) $$

The parabolic damage equation is considered as a rate-dependent governing equation with a viscous contribution introducing a delay in the solution. Asymptotically, the rate-independent law is restored when $\mu_d$ vanishes.

### 2.2.3. Hyperbolic damage formulation

In this case, the damage kinetic energy is considered. This damage kinetic energy contributes to the degradation process and accounts for the effect of microscopic accelerating forces. The dissipative process can be taken into account as in the parabolic case, it contributes to damp the wave-like behavior of the system when its evolution is driven by such an hyperbolic PDE. The damage governing equation results in:

$$ \int_t \int_{\Omega} \left( \rho_d \ddot{d} + \frac{\partial \Psi}{\partial d} + G_c \frac{\partial \gamma}{\partial d} - G_c \underline{\nabla} \cdot \left( \frac{\partial \gamma}{\partial \underline{\nabla} d} \right) \right) \delta d \, d\Omega \, dt \geq \int_{\Omega} \int_t - (\mu_d \dot{d} \delta d) \, dt \, d\Omega, \forall \delta d \in \bar{\Omega}_d \qquad (26) $$

$$ \rho_d \ddot{d} + \mu_d \dot{d} \geq -\frac{\partial \Psi}{\partial d} - \frac{G_c}{\ell_c} (d - 2 \ell_c^2 \Delta d) = f(Y, d). \qquad (27) $$

## 2.3. Tension-compression asymmetry models

To take into account the tension/compression asymmetric phenomenon in the fracture process, many strain energy splits are proposed in the literature. The purpose of this split is to decompose the strain energy density into a contribution affected and thus contributing to the damage process and a second one independent on the damage variable:

$$ \Psi(\underline{\epsilon}, d) = \frac{1}{2} \left( g(d) \underline{\epsilon}^{+} : C : \underline{\epsilon}^{+} + \underline{\epsilon}^{-} : C : \underline{\epsilon}^{-} \right) = g(d) \Psi^{+}(\underline{\epsilon}) + \Psi^{-}(\underline{\epsilon}). \qquad (28) $$
The split is based on the definition of a positive and negative strain $\underline{\epsilon}^{\pm}$ and their corresponding stress $\underline{\sigma}^{\pm}$ and strain energy density $\Psi^{\pm}$. The main requirement for such a split is the orthogonality between the strain and stress tensors of opposite sign:

$$ \underline{\sigma}^{+} : \underline{\epsilon}^{-} = 0, \quad \text{and} \quad \underline{\sigma}^{-} : \underline{\epsilon}^{+} = 0. \qquad (29) $$

$\frac{\partial \Psi}{\partial d}$ can then be rewritten as

$$ \frac{\partial \Psi}{\partial d} = g'(d) \Psi^{+} = -2(1-d) \Psi^{+}. \qquad (30) $$

The most used in the numerical fracture analysis are the volumetric and deviatoric split [37–39], the spectral decomposition of the strain tensor [14, 40, 41] or the masonry-like material model [42]. These energy splits were initially developed for isotropic material. For anisotropic constitutive material models, the orthogonality condition is not guaranteed.

To handle anisotropic behaviors, an orthogonal decomposition model has been developed, based on the Kelvin decomposition of the stiffness tensor $C$ or the compliance tensor $S$ [43, 44]. In [43], it is defined as a quadratic decomposition of the free energy using the square root stiffness tensor $C$:

$$ \underline{\epsilon} : C : \underline{\epsilon} = \underline{\tilde{\epsilon}} : \underline{\tilde{\epsilon}} \qquad (31) $$

with $\underline{\tilde{\epsilon}} = C^{1/2} \underline{\epsilon}$.

The strain energy is then defined introducing a damage asymmetric behavior through a damage degradation function $g$ which is further applied on the active part of the energy density.

This approach has been applied in phase field modelling for fracture considering a single scalar field $d$:

$$ \Psi = \frac{1}{2} g(d) \underline{\tilde{\epsilon}}^{+} : \underline{\tilde{\epsilon}}^{+} + \underline{\tilde{\epsilon}}^{-} : \underline{\tilde{\epsilon}}^{-}. \qquad (32) $$

As a consequence of this split, the strain energy definition is modified following Eq. (28).

A detailed 2D analytical analysis of this energy split model was carried out in [45] leading to its first application in phase field fracture model in [46]. This approach opens up the possibility to extend some decomposition initially developed for isotropic material to anisotropic material [47, 48], such as the volumetric/deviatoric split model and the no-tension split model [39, 42]. Recently, another application of this decomposition in the context of heterogeneous material was run in a quasi-static regime in [49].

Although anisotropic materials are not considered herein, the orthogonal split model formulated with a scalar damage variable from Eq. (32) is used in this work as a generic framework for implementing damage asymmetric effect, in 2D case.

# 3. Resolution strategy and Implementation

## 3.1. Discretization

The weak form of the equations of the coupled problem is expressed as $\forall t \in T$:

$$ \int_{\Omega} \left( \rho_{u} \underline{u} \underline{u}^{*} + \underline{\sigma} : \underline{\nabla} \underline{u}^{*} \right) d\Omega = \int_{\Omega} \underline{b} \underline{u}^{*} d\Omega + \int_{\partial\Omega} (\underline{\sigma} \underline{n}_{s}) \underline{u}^{*} dS \quad \forall \underline{u}^{*} \in \tilde{\Omega}_{u}(\underline{0}), \qquad (33) $$
$$ \int_{\Omega} \left( \rho_d \ddot{d} d^* + \mu_d \dot{d} d^* + 2 d d^* \Psi^+ + G_c \left( \frac{d}{\ell_c} d^* + \ell_c \underline{\nabla} d \underline{\nabla} d^* \right) \right) d\Omega \\ = \int_{\Omega} 2\Psi^+ d^* d\Omega \quad \forall d^* \in \tilde{\Omega}_d, \qquad (34) $$

In these equations, $u^*$ and $d^*$ are test functions for the displacement and damage fields respectively that belong to the admissible space defined by

$$ \tilde{\Omega}_u(\underline{v}) = \{ \underline{u} : \Omega \to H^1(\Omega, R^{dim}), \underline{u} = \underline{v} \text{ on } \partial v \Omega \} \qquad (35) $$

and

$$ \tilde{\Omega}_{\bar{d}} = \{ d : \Omega \to R \}. \qquad (36) $$

Considering the standard finite element framework to solve the weak form of the equations, $\Omega$ is discretised with a set of finite elements supporting shape functions. The displacement and the damage field are approximated by

$$ \underline{u} = \underline{N} \underline{U} , \quad d = N \underline{D} \qquad (37) $$

with $\underline{N}$ and $\underline{N}$ gathering vectors and scalar shape functions. Vectors like $\underline{U}$ or $\underline{D}$ collect the values of the nodal degrees of freedom for each field variable. Using finite element shape functions as test functions, Equations (33)–(34) can therefore be collected in a matrix vector format such as

$$ \underline{M}_u \ddot{\underline{U}} + \underline{F}_u^i (\underline{U}, \underline{D}) = \underline{F}_u^e \qquad (38) $$

for the displacement problem and

$$ \underline{M}_d \ddot{\underline{D}} + \underline{C}_d \dot{\underline{D}} + \underline{K}_d \underline{D} + \underline{F}_d^i (\underline{U}, \underline{D}) = \underline{0} \qquad (39) $$

for the damage problem. Note that, the boundness and irreversibility conditions are now additional constraints to be applied on the numerical solution (see section 3.2 and section 3.3). In these equations, the subscripts $u$ or $d$ refer to type of problem to be solved and the superscripts $i$, $v$, or $e$ to the type of generalized finite element forces (internal, viscous, external). The expression of the finite element matrices and vectors are:

$$ \underline{M}_u = \int_{\Omega_e} \rho_u \underline{N}^T \underline{N} d\Omega \qquad (40) $$

$$ \underline{M}_d = \int_{\Omega_e} \rho_d \underline{N}^T \underline{N} d\Omega \qquad (41) $$

$$ \underline{C}_d \dot{\underline{D}} = \int_{\Omega_e} \mu_d \underline{N}^T \underline{N} \dot{\underline{D}} d\Omega = \underline{F}_d^v (\dot{\underline{D}}) \qquad (42) $$

$$ \underline{F}_u^i (\underline{U}, \underline{D}) = \int_{\Omega_e} \underline{B}^T : \underline{\sigma} d\Omega \qquad (43) $$

$$ \underline{K}_d = \int_{\Omega_e} \left( 2 \left( \Psi^+ + \frac{G_c}{2l_c} \right) \underline{N}^T \underline{N} + G_c l_c \underline{B}^T \underline{B} \right) d\Omega \qquad (44) $$

$$ \underline{F}_d^i (\underline{U}) = - \int_{\Omega_e} 2\underline{N} \Psi^+ d\Omega. \qquad (45) $$

where $\underline{B}$ and $\underline{B}$ contains the vector and scalar shape function derivatives.

At each new time increment $t_{n+1}$, a staggered resolution strategy, initially proposed by [50], is usually adopted in the literature for solving the discretized coupled problem. In the following, the displacement problem is solved before the damage problem. This has actually no influence on the final result. The only consequence is that one has to chose to show, in the post-processing step, the displacement at time $t_{n+1}$ either with the damage at time $t_n$ or $t_{n+1}$. Usually it is necessary to loop over the alternate resolution of both problems until convergence is reached
following a given criterion. We will see later on that this is not required in the considered cases. Concerning time discretization, we focus on the use of explicit time stepping strategies. For the displacement problem, a Central Difference method (CD) is used. Considering a time step $\Delta t$, this leads to the following algorithm where subscripts $n$ or $n+1$ define quantities at the previous time step and those to be updated:

<caption>**Algorithm 1** Central difference scheme for the displacement resolution</caption>

1: **Predictor**

2: $\dot{U}_{n+1}^{pred} \leftarrow \dot{U}_n + 0.5 \Delta t \ddot{U}_n$

3: $U_{n+1} \leftarrow U_n + \Delta t \dot{U}_n + 0.5 \Delta t^2 \ddot{U}_n$

4: **Solve motion equation at** $t_{n+1}$

5: $M_u \ddot{U}_{n+1} \leftarrow F_u^e(U_{n+1}) - F_u^i(D_n, U_{n+1})$

6: **Corrector**

7: $\dot{U}_{n+1} \leftarrow \dot{U}_{n+1}^{pred} + 0.5 \Delta t \ddot{U}_{n+1}$

It is worth noting that in this particular case of an explicit CD scheme, the displacement vector $U_{n+1}$ at a new time instant is defined from displacement, velocity and acceleration vectors at $t_n$ only. As a consequence, the damage at $t_{n+1}$ does not affect the computation of $U_{n+1}$. In practice, the mass matrix is condensed in a vector format (it is actually made diagonal) using a row-sum lumping method, so that the displacement update does not require a linear system solver. Concerning the damage problem, different options have to be considered depending on the chosen formulation. They are summarized in Table 1.

<table><caption>**Table 1.** Resolution strategy</caption><thead><tr><th rowspan="2">Strategy</th><th colspan="2">PDE type</th><th colspan="2">Time stepping</th></tr><tr><th>u</th><th>d</th><th>u</th><th>d</th></tr></thead><tbody><tr><td>1</td><td rowspan="3">Hyperbolic</td><td>Elliptic</td><td rowspan="3">CD (Algorithm 1)</td><td>—</td></tr><tr><td>2</td><td>Parabolic</td><td>Forward Euler (Algorithm 2)</td></tr><tr><td>3</td><td>Hyperbolic</td><td>Verlet (CD) (Algorithm 3)</td></tr></tbody></table>

For each combination that will be detailed in the next section, the use of Algorithm 1 makes the damage and the displacement problems independent. Therefore, there is no need to perform an alternate minimization resolution and the latter can be replaced by a one-pass staggered resolution.

## 3.2. Resolution of the damage problem

The various types of damage PDE for the phase field model presented in the previous Section require different resolution schemes (see Table 1) that will be presented in this subsection. Before all, it is worth noticing that damage field nodal values are shrunken between 0 and 1 to enforce boundedness.

### Elliptic damage formulation

For the elliptic damage formulation, an implicit method is used for the resolution of the damage problem which could be refactored as follows:

$$ K_d \mathbf{D} + \mathbf{F}_d^i(\mathbf{U}) = \mathbf{0}. \qquad (46) $$
This mixed explicit/implicit approach is disadvantageous for large scale problems. In this case, *implicit* refers to the fact that a linear system of equations has to be solved for the damage problem. This system of equations is solved at each time step, leading to large computational time because the explicit time stepping scheme for the displacement problem requires very small time steps. In that perspective, the other damage PDEs are solved using explicit time integration schemes, leading to a fully explicit resolution framework.

### Parabolic damage formulation

In this case, the damage problem reduces to

$$ C_d \dot{\mathbf{D}} + \mathbf{K}_d \mathbf{D} + \mathbf{F}_d^i(\mathbf{U}) = \mathbf{0}. \quad (47) $$

The problem is integrated with the forward Euler scheme, (Algorithm 2). A lumped damping matrix $C_d$ is used so that the linear system of equations is diagonal.

### Algorithm 2 Forward Euler explicit scheme

$$ \text{1: } \mathbf{C}_d \dot{\mathbf{D}}_{n+1} \leftarrow -\mathbf{K}_d \mathbf{D}_n - \mathbf{F}_d^i(\mathbf{U}_{n+1}) $$

$$ \text{2: } \mathbf{D}_{n+1} \leftarrow \mathbf{D}_n + \Delta t \dot{\mathbf{D}}_{n+1} $$

### Hyperbolic damage formulation

When considering the hyperbolic formulation, the problem to solve is driven by Eq. (39). Here, the latter is solved with the Verlet method (Algorithm 3), relying on the definition of the half-step velocity [51, §13.4.2]. The viscous force $C_d \dot{\mathbf{D}}$ is computed at mid-step ([52], § 7.5.1) without any simplification. Here, the mass matrix $M_d$ considered, is lumped by the row-sum method.

### Algorithm 3 Verlet method with predictor-corrector form [53]

1: **Predictor**

2: $\dot{\mathbf{D}}_{n+1}^{\text{pred}} \leftarrow \dot{\mathbf{D}}_n + 0.5 \Delta t \ddot{\mathbf{D}}_n$

3: $\mathbf{D}_{n+1} \leftarrow \mathbf{D}_n + \Delta t \dot{\mathbf{D}}_n + 0.5 \Delta t^2 \ddot{\mathbf{D}}_n$

4: **Solve motion equation at** $t_{n+1}$

5: $\mathbf{M}_d \ddot{\mathbf{D}}_{n+1} \leftarrow -\mathbf{C}_d \dot{\mathbf{D}}_{n+1}^{\text{pred}} - \mathbf{K}_d \mathbf{D}_{n+1} - \mathbf{F}_d^i(\mathbf{U}_{n+1})$

6: **Corrector**

7: $\dot{\mathbf{D}}_{n+1} \leftarrow \dot{\mathbf{D}}_{n+1}^{\text{pred}} + 0.5 \Delta t \ddot{\mathbf{D}}_{n+1}$

### 3.3. Damage irreversibility constraint approximation

The irreversible damage process is not intrinsically taken into account in the PDEs. In order to prevent damage healing while unloading, several methods were proposed. A comprehensive review of these methods can be found in the literature [54]. Since the targeted simulations are to rely on a fully explicit resolution, i.e. without an alternate minimization iterative loop, the aim is to favour methods that do not require linear systems resolution, such as ones involving augmented Lagrangian [55, 56].
## History variable

An alternative to the previous strategies is to modify the energetic crack driving force $Y$ (Eq. (21)). For instance, it was suggested replacing $\Psi^+$ with a local history variable $H$ such as $Y = g'(d)H(\underline{u}, t)$ to constrain the monotonic growth of the damage evolution [14]. The classical formulation is such as Eq. (48), where $\Psi_s$ is a threshold (set to 0 in the following):

$$H(\underline{u}, t > \tau) = \max_{\tau \in [0, t]} \Psi^{+}(\underline{\epsilon}(\underline{u}), \tau) \quad \text{or} \quad \max_{\tau \in [0, t]} (H(x, t), \Psi^{+}(\underline{\epsilon}(\underline{u}), \tau) - \Psi_s). \quad (48)$$

This strategy has been tested in the context of hyperbolic damage formulation [25,27]. It is worth noting that for the latter papers, no violation of the irreversibility constraint are reported. However, many other forms of the history variable, not studied here, are also proposed in the literature and are compared in [18, 57, 58], based on stress criterion (Rankine model, Tresca model and so on). This method is very efficient, but the classical form proposed in [14] does not allow for the irreversibility condition to be met exactly, in particular for hyperbolic damage formulation. This fact will be illustrated in Section 4. However, it will be used in the following work for the three damage PDEs as a reference strategy. While complex numerical methods have to be implemented to strictly satisfy the irreversibility condition in the case of the elliptic formulation, using parabolic and hyperbolic formulations offers the opportunity to meet this condition for the discrete problem at a very low cost.

## Parabolic damage formulation

In the case of the parabolic formulation, an alternative option to the history variable consists in taking the positive part of the yield function $f(Y, d)$ to include this condition in the strong form of the continuous Cauchy problem (Eq. (25)):

$$\mu_d \dot{\mathbf{d}} = \langle f(Y, d) \rangle_+ = \left\langle -g'(d)\Psi^+ - \frac{G_c}{\ell_c^2} (d - 2 \ell_c^2 \Delta d) \right\rangle_+ \quad (49)$$

Even if this modification cannot be included in the weak form of Eq. (25), the latter is implemented in the discretized system of equations (47) in a purely algorithmic way. Therefore, Step 1 of Algorithm 2 is updated as follows:

$$C_d \dot{\mathbf{D}}_{n+1} = \langle -\mathbf{K}_d \mathbf{D}_n - \mathbf{F}_d^i (\mathbf{U}_{n+1}) \rangle_+ \quad (50)$$

where the Macaulay brackets are applied component wise to vector of nodal forces.

## Hyperbolic damage formulation

Similar to the parabolic damage PDE, we propose for the hyperbolic damage one, a strategy that can be implemented easily in the Verlet algorithm. Indeed, when updating the damage velocity $\dot{\mathbf{D}}_{n+1}$, one can prevent negative values by considering only the positive part of the correction using Macaulay brackets:

$$\dot{\mathbf{D}}_{n+1} = \langle \dot{\mathbf{D}}_{n+1}^{\text{pred}} + 0.5\Delta t \ddot{\mathbf{D}}_{n+1} \rangle_+ \quad (51)$$

The acceleration is then modified accordingly such as $\ddot{\mathbf{D}}_{n+1} = \frac{2}{\Delta t}(\dot{\mathbf{D}}_{n+1} - \dot{\mathbf{D}}_{n+1}^{\text{pred}})$. In fact, this strategy is similar to what it is proposed for the parabolic case.

## 3.4. Numerical stability of the time integrators

When the resolution of a coupled set of equations involves (fully) explicit time integrators, the stability of each related schemes have to be ensured. In the current study, the critical time step is therefore such as:

$$\Delta t_{cr} = \min(\Delta t_{cu}, \Delta t_{cd}), \quad \Delta t = c_{sec}\Delta t_{cr}. \quad (52)$$
where $\Delta t_{cu}$ and $\Delta t_{cd}$ are the critical time steps for the displacement and the damage problem respectively, $c_{sec}$ is a safety coefficient and $\Delta t$ the time step which is eventually used. Assuming a 1D bar of length $L$ discretized with a P1 element and considering a lumped mass matrix, the critical time step related to the displacement PDE is the classical CFL expression

$$ \Delta t_{cu} = \frac{h_{min}}{c}, \quad \text{with} \quad c = \sqrt{\frac{K + 4\mu/3}{\rho}} \qquad (53) $$

where $h_{min}$ is the smallest element size, $\mu$, $K$ are the shear and bulk modulus respectively and $c$ is the compressional-wave speed (or P-wave). For the parabolic damage formulation, the critical time step associated to the forward Euler explicit scheme is obtained from a stability analysis of a single degree of freedom problem leading to the following expression:

$$ \Delta t_{cd} \approx \frac{\mu_d h_{min}^2}{G_c \ell_c}. \qquad (54) $$

Theoretically speaking, the estimation of the critical time step requires calculating the spectral radius of an amplification matrix which depends on the mass and stiffness matrices and, consequently, on the shape functions, the element types and so on. In practice, a solver can either employ a CFL tool to rigorously compute the spectral radius at each time step, which is computationally expensive, or rely on assumptions to approximate an expression. Therefore, a safety coefficient $c_{sec}$ can be used.

The critical time step associated to the hyperbolic damage formulation using the Verlet algorithm [59] is such as:

$$ \Delta t_{cd} = \frac{2(\sqrt{\xi^2 + 1} - \xi)}{\omega} = f(\rho_d, \mu_d) \qquad (55) $$

with

$$ \omega = \frac{2}{h_{min}} \sqrt{\frac{G_c \ell_c}{\rho_d}}, \quad \xi = \frac{\mu_d}{2\omega\rho_d} \qquad (56) $$

The latter are defined the natural frequency and the damping ratio of the mode related to the maximal eigenvalue of the amplification matrix. A more detailed stability analysis is provided in Appendix A and B.

### 3.5. Software framework

In this work, each resolution mostly relies on the open-source software FEniCS-dolfinx, a cross-platform software used for solving partial differential equations through finite element method [60]. Furthermore, the implicit strategy is implemented owing to the PETSc library [61]. Lastly, material behaviour knowledge is handled by the code generator MFront [62]. In practice, the latter is used to integrate the constitutive law, i.e. the orthogonal split model of the traction compression asymmetry of the crack for the linear elastic case.

## 4. Dynamic failure of a bar under tension

The purpose of this first example is to perform a critical comparative analysis between the hyperbolic and parabolic phase field damage equations with the standard elliptic damage formulation. The influence of additional parameters is evaluated as well as the efficiency of the irreversibility constraint implementations and the validity of the critical time step evaluation presented in the previous section.
## 4.1. Numerical model

The geometry and the boundary conditions of the tensile bar are described in Figure 1. Only one quarter of a bar submitted to a tensile stress at both ends is modeled for reasons of symmetry. Dimensions of the considered simulation domain are reported in Figure 1. Along the left edge, the horizontal displacement is set to zero, as well as the vertical displacement along the bottom edge. A positive stress $\sigma_a$ is applied in the horizontal direction on the right edge. Its value is kept constant from the first time step so that it is a one step function in time.

The simulations are performed with the mechanical properties of a random steel-like material given in Table 2, except for the Poisson's ratio that is set to 0 in order to obtain a 1D stress state in the bar. In the case where the material is elastic, the reaction force (to the symmetry condition at the left edge) remains 0 until it reaches a plateau at $t = t_L = L\sqrt{\rho_u/E} = 0.205$ ms. The value of this plateau that lasts $2t_L$ is the applied stress $\sigma_a$. Then the reaction force vanishes. The bar is meshed with 20×400 quadrangular elements using linear interpolation and a four Gauss points quadrature rule. The material strength is estimated by the theoretical quasi-static fracture strength formula $\sigma_R = \frac{3}{16}\sqrt{\frac{3EG_c}{\ell_c}}$. As expected in transient loading, inertia effects require a lower applied load magnitude to break the bar and the value we considered for $\sigma_a$ (150 MPa) is 0.62 times lower than this quasi-static strength.

<table><caption>**Table 2.** Modelling parameters.</caption><thead><tr><th>h<br>[m]</th><th>ℓ<sub>c</sub><br>[m]</th><th>E<br>[GPa]</th><th>G<sub>c</sub><br>[N/m]</th><th>ν</th><th>ρ<sub>u</sub><br>[kg/m<sup>3</sup>]</th><th>σ<sub>R</sub><br>[MPa]</th></tr></thead><tbody><tr><td>0.0025</td><td>0.0075</td><td>190</td><td>22130</td><td>0</td><td>8000</td><td>243.2</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Tensile bar test - geometry and boundary conditions</figcaption></figure>

The values of the time step for the simulations are presented in section 3.4, eq. (52).

## 4.2. Definition of some metrics for deviation assessment

Physically based metrics are introduced to perform a quantitative comparison between the parabolic/hyperbolic damage PDE and the elliptic damage PDE. It measures the deviation of the solution $x_{sol}$ with respect to a reference solution $x_{ref}$, based on the mean squares difference (MSD):

$$ \text{MSD}_{\text{rel}} = \frac{\text{MSD}(x_{\text{sol}})}{\text{MSD}(0)} \quad \text{with} \quad \text{MSD}(x_{\text{sol}}) = \frac{1}{N} \sum_{n=1}^{N} (x_{\text{sol}} - x_{\text{ref}})^2, \qquad (57) $$

where N represents the set of points over which the comparison is made.
The MSD function is computed for two physical quantities: (a) the normal component of the reaction $R_x(t)$, calculated on the clamped edge of the bar and (b) the damage profile along the bar at a fixed time step $t_f$, identical for all the simulations, $d(x, t_f)$. In the remainder of the section, the results obtained with the elliptic formulation will be considered as the reference.

### 4.3. $\mu_d$-viscosity sensitivity analysis of parabolic damage PDE

In this section, a sensitivity analysis upon the damage viscosity of the parabolic PDE is carried out. Two ways to constrain the damage irreversibility are considered, provided previously in Eq. (50) and Eq. (48).

#### 4.3.1. Parabolic formulation: damage irreversibility with Macaulay brackets

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Damage irreversibility with Macaulay brackets: sensitivity of the reaction $R_x(t)$ to $\mu_d$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Damage irreversibility with Macaulay brackets: sensitivity of the damage profile along the bar at a fixed time $t_f$ to $\mu_d$.</figcaption></figure>

The parabolic formulation (Eq. (49)) is first tested with the Macaulay brackets to ensure the damage irreversibility. Figure 2 illustrates the reaction at the clamped edge and Figure 3 the damage profiles along the bar, for various $\mu_d$ values. Besides, the time at which a zero stress is reached (it can be considered as a time to failure if it is lower than $3t_L$) increases with increasing $\mu_d$ values, denoting a delay in the crack onset. This delay is also pointed out in the damage profiles shown in Figure 3(a). At the fixed time $t_f$, damage profiles show distinct states: increasing $\mu_d$ produces lower damage value at the left edge and confirms the delay effect mentioned before.

For high $\mu_d$ values ($\ge 15.0$ Pa.s), some oscillations are obtained on the reaction evolution just before it reaches a plateau. These oscillations are typically obtained for a fully elastic simulation (without damage). This observation is consistent with the damage profiles in Figure 3 where high $\mu_d$ values leads to a delayed (in time) and more diffuse increase of $d$ in the vicinity of the left bar edge. It is also worth noting that for very low values of $\mu_d$ ($\le 2.0$ Pa.s), the time to failure (represented by the time at which the reaction vanishes) is shortened compared to the prediction of the elliptic formulation. Let us remind that for the elliptic formulation, the irreversibility constraint is prescribed using the local history variable whereas for the results of the parabolic formulation presented in this subsection, this condition is more rigorously imposed with the use of the Macaulay brackets. This is the reason why the results of the parabolic formulation do not collapse with those of the elliptic formulation even for a vanishing damping coefficient $\mu_d$. However, a good agreement is obtained between these two formulations for $\mu_d \approx 3.0$ Pa.s. This is confirmed by Figure 4 that presents the metrics defined in Section 4.2. Indeed, the discrepancy assessed on the reaction evolution, in Figure 4(a), clearly shows a minimal value with the mechanical response of the elliptic damage results, which occurs at approximately 3.6 Pa.s. Regarding the damage profile, Figure 4(b), the minimal value shifts slightly to a higher value of $\mu_d$ compared to the criteria computed with the reaction force. These plots allow for assessing the differences between the prediction of the elliptic and the parabolic formulation. They quantitatively confirms that, for low viscosity, there is a good agreement between those two formulations whereas increasing $\mu_d$ leads to higher discrepancies. Note that the typical upper value of $\mu_d$ below which the two formulations give consistent results is different when considering the reaction force evolution with time or the damage profile. This bound is estimated to 5 Pa.s for the reaction and 25
Pa.s for the damage profile.

For a practical perspective, one would not want to run simulations with a time step lower than the time step of the displacement problem. In the plots of Figure 4, a black vertical line is drawn at the value of $\mu_d$ for which the critical time step of the damage problem equals the time step of the displacement problem. Not using a time step below $\Delta t_{cd}$ means that only results from the parabolic formulation deviating from the prediction of the elliptic one by about 40% for the reaction and 25% for the damage profile can be obtained. On one side, this is not a problem as one does not want to replicate the results of the elliptic formulation with the parabolic one, the parabolic formulation being considered as an enrichment of the elliptic one. On the other side, one would have though to take the opportunity to use the parabolic formulation with low $\mu_d$ values (typically less than 5 Pa.s), to replicate the prediction of the elliptic formulation at a much lower computational cost. As explained in the previous section, running mixed explicit/implicit simulation comes at a very expensive cost as a linear system resolution is required at each time step whereas with the explicit time stepping used for the resolution of the parabolic formulation, the cost of the damage problem resolution is drastically decreased. Furthermore, the proposed strategy using Macaulay brackets in the discretized problem offers the opportunity to reach the irreversibility constraint more efficiently. However, Eq. (54) shows that the critical time step for the damage problem using the parabolic formulation scales linearly with $\mu_d$. Thus, decreasing $\mu_d$ until the discrepancy between the predictions of parabolic formulation and elliptic formulation is below 1% comes at a reasonable additional cost ($\sim 4\%$) compared to that of an explicit/implicit strategy.

<figure><img src="image_5.png" /><figcaption>**Figure 4.** MSD versus the damage viscosity $\mu_d$. (4(a)) used the reaction step responses. (4(b)) used the damage profiles taken at the time $t_f$ of the fully crack stage of the elliptic damage result.</figcaption></figure>

### 4.3.2. Parabolic formulation: damage irreversibility with local history variable

The simulations are now performed with the local history variable $H$ [14] for prescribing the irreversible evolution of $d$. The reaction force of the tensile test is reported for different values of the damage viscosity $\mu_d$ (Figure 5). The time evolution shows an abrupt transition of the behaviour of the numerical simulation results between 21.5 Pa.s and 21.52 Pa.s. Beyond the latter value, the step response is recovered but with an important delay effect as obtained previously with the Macaulay brackets. For $\mu_d$ below this typical value of 25 Pa.s, the response is completely
diverging. Note that even for larger viscosity coefficient, strong oscillations appear in the reaction force evolution after the time to failure is reached. To try and explain these numerical instabilities, the irreversibility of the damage evolution is assessed in Figure 6. In the plots presented in this Figure, space time maps of the damage velocity $\dot{d}$ are displayed for 4 different values of $\mu_d$. For $\mu_d$ = 21.52 Pa.s, $\dot{d}$ is positive at any location and at any time, meaning that irreversibility is ensured. For lower values of $\mu_d$, negative values are obtained meaning that damage irreversibility is violated. Starting with the lower value of $\mu_d$, negative value of damage velocity are obtained at the right end of the bar where the load is applied since the very beginning of the simulation. When $\mu_d$ is increased, the onset of this artifact is delayed until $t = 0.1$ ms for $\mu_d$ = 21.5 Pa.s. So, if damage damping (through the dissipative terms of the parabolic formulation parametrized by $\mu_d$) is strong enough, any violation of the irreversibility condition is smoothed out while for lower damping, this artifact becomes an instability. The occurrence of this instability is amplified by the coupling between the damage problem with the hyperbolic displacement problem that help transmitting the numerical instability in the bar. One could argue that this occurs because the stability condition of the forward Euler time integration is not satisfied. The safety factor in the critical time step evaluation (see Eq. (52)) is thus reduced down to 0.4, as shown in Figure 7. This improves the results but the reaction still exhibits significant high-frequency oscillations upon bar failure, despite the theoretical numerical stability of the damage time integration scheme for the chosen time step. In contrast, when using Macaulay brackets, no such oscillations are observed at failure. The reaction force drops to zero and remains null thereafter. The latter are thus due to the way the local history variable constrains irreversibility.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** History-variable-based simulations - Reaction forces at the clamped edge of the bar</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** $\langle \dot{d} \rangle_{-}$ maps in the $(x, t)$-plane for different $\mu_d$ values. The upper limit of color scale range is set 0 to easily assess damage irreversibility.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Reaction response for two different $\mu_d$. For each one, the simulation is performed with three safety coefficients to set the time step $\Delta t = c_{sec}\Delta t_{cd}$, $\Delta t_{cd}$ being lower than $\Delta t_{cu}$ in these cases.</figcaption></figure>
4.4. ($\rho_d$, $\mu_d$) sensitivity analysis of hyperbolic damage PDE

In order to assess hyperbolic damage formulation predictions, 2430 simulations are run, involving 54 and 45 values for micro-inertia $\rho_d$ and micro-damping $\mu_d$ respectively. The MSD (Eq. (57)), evaluated for the damage profile or reaction force, are presented using maps in the ($\rho_d$, $\mu_d$)-plane. Each rectangular colored element corresponds to one evaluation point ($\rho_d^i$, $\mu_d^j$). As for the parabolic formulation we also investigate the effect of the implementation of the irreversibility constraint: first using the Macaulay brackets on the discretized problem (Eq. (51)) and then using the local history variable (Eq. (48)).

4.4.1. Hyperbolic formulation: damage irreversibility with Macaulay brackets

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Damage MSD map and reaction MSD map for the hyperbolic formulation using Macaulay brackets.</figcaption></figure>

The analysis of the results starts by reporting on the MSD maps presented in Figure 8. For the two physical quantities used to assess the difference between the elliptic formulation and the hyperbolic one, region of the ($\rho_d$, $\mu_d$)-plane where the prediction of the two formulations are in good agreement are clearly distinguished: it is roughly defined by $\mu_d \le 5$ Pa.s and $\rho_d \le 1.0 \ e^{-4} \ kg/m$. Either when $\mu_d$ or $\rho_d$ is increased from these limits, the deviation increases. It appears that $\rho_d$ is the most influential parameter on the two MSD.

Before going into a detailed analysis of these deviation maps by examining the reaction evolution and damage profiles, one has to comment on the practical use of arbitrary parameters in fully explicit resolution of the hyperbolic displacement and damage formulation. Indeed, as mentioned above for the parabolic formulation, the critical time step of the damage problem that is illustrated in Figure 9 decreases below the critical time step of the displacement problem when decreasing $\rho_d$ and/or increasing $\mu_d$. The region of the ($\rho_d$, $\mu_d$)-plane where the critical time step of the damage formulation is lower than that of the displacement problem is delimited by a white
dashed line in Figure 8a. The same comments as for the parabolic formulation hold: decreasing $\rho_d$ and $\mu_d$ allows for replicating the prediction of the elliptic formulation but with a time step lower than that of the displacement problem. However, even if one has to decrease the time step and thus to increase the computational cost, the additional expense is much lower than that of an explicit/implicit resolution strategy.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Critical time step for the hyperbolic damage PDE as a function of $\rho_d$ for different values of $\mu_d$.</figcaption></figure>

In Figure 8b, six black vertical lines with Roman numerals are positioned at typical $\rho_d$ values. The reaction evolution and damage profiles for different $\mu_d$ values along these lines (grouped by three) are now analyzed. Starting from lines I ($\rho_d$ = 9.28 e⁻⁷ kg/m) and II ($\rho_d$ = 4.0 e⁻⁶ kg/m) in Figure 10(a), 10(b), 10(c) and 10(d), the plots confirm very good agreement with the elliptic formulation for low $\mu_d$ values. Increasing $\mu_d$ produces the same effect as for the parabolic formulation: a delay effect is introduced that increases the time to failure and produces oscillations in the reaction force similar to the elastic case for the highest $\mu_d$. These oscillations are caused by the delay in crack nucleation that is also observed on the damage profile at $t_f$ = 0.41 ms, the value of $d$ at the left bar edge being below 0.6 for $\mu_d$ ≥ 10 Pa.s. For line III ($\rho_d$ = 1.2 e⁻⁵ kg/m), due to the increase of $\rho_d$ the results start deviating significantly from the elliptic formulation predictions. However, line III in Figure 8b appears as a hill in the MSD maps. Our interpretation is that this line is really close to the white line along which the critical time step for the damage problem and the displacement problem are identical. This implicitly means that the elastic stress wave and
the wave transporting damage (governed by the same PDE type as the displacement does) travels at the same speed. Further, in these simulations damage initiate all along the bar from the beginning of the simulation due to the incident stress wave (this is confirmed by the reaction plateau that is below the plateau of the elastic simulation, see Figure 10(e)). The stress wave propagation is thus slightly lower than in the undamaged material for which the critical time step is estimated. In other words, we think that for this particular values of ($ρ_d, μ_d$) stress wave and damage wave travel at the same speed which means that the front of these waves coincide. The synchronous propagation of the damage and stress wave front then produces an increase in damage along the bar (see Figure 10f). This phenomenon is specific to the parameter values along line III. It significantly differs from what is observed for line I and II and also for line IV, V and VI as it will be commented later on.
<figure><img src="image_11.png" /><figcaption>(a) $\rho_d = 9.28 \times 10^{-7} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) $\rho_d = 9.28 \times 10^{-7} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(c) $\rho_d = 4 \times 10^{-6} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(d) $\rho_d = 4 \times 10^{-6} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(e) $\rho_d = 1.2 \times 10^{-5} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(f) $\rho_d = 1.2 \times 10^{-5} \text{ kg/m}$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Reaction force evolution and damage profile for ($\rho_d$, $\mu_d$) along line I, II and III of Figure 8b for the hyperbolic formulation with Macaulay brackets.</figcaption></figure>

Moving to lines IV ($\rho_d = 1.5 \times 10^{-4} \text{ kg/m}$), V ($\rho_d = 6.0 \times 10^{-4} \text{ kg/m}$) and VI ($\rho_d = 3.95 \times 10^{-3} \text{ kg/m}$), the MSD for the reaction force and the damage profile significantly increases (see Figure 8). First, the reaction plateau now starts at the same value as the plateau of the elastic simulation (see Figure 11(a), 11(c) and 11(e)). This is because the damage micro-inertia prevents damage along the bar under the stress amplitude of the incident wave, this can be interpreted as an increase of the
apparent strength of the material. Second, one observes oscillations in the reaction evolution after the first peak at $t_L$. These oscillations period increases when $\rho_d$ is increased from $1.5 \times 10^{-4}$ kg/m (line IV) to $6.0 \times 10^{-4}$ kg/m (line V) and then to $3.95 \times 10^{-3}$ kg/m (line VI) and their amplitude decreases when $\mu_d$ is increased. This confirms that for the range of $\rho_d$ from line IV to line VI, damage evolution is mainly governed by (micro-)inertia effects. Oscillations also appear in the damage profile observed in Figure 11(d) and 11(f). The period and amplitude reveals the same dependency to the parameters as those in the reaction evolution. Note that for the simulations with the parameters values corresponding to line V and VI, the speed of the damage wave is much lower than that of the stress wave. For line IV, while the reaction evolution (Figure 11(a)) follows the same trend as for line V and VI (Figure 11(c) and 11(e)), the damage profile are different: there is no oscillation at the left edge but a peak (of decreasing amplitude when $\mu_d$ is increased) appears near the right edge. This occurs when the damage wave travelling from left to right interacts with the stress wave travelling from right to left (its speed being slightly higher).
<figure><img src="image_18.png" /><figcaption>**Figure 11.** Reaction force evolution and damage profile for $(\rho_d, \mu_d)$ along line IV, V and VI of Figure 8b for the hyperbolic formulation with Macaulay brackets.</figcaption></figure>

### 4.4.2. Hyperbolic formulation: damage irreversibility with local history variable

The results of the simulations using the hyperbolic formulation and the local history variable (Eq. (48)) for prescribing damage irreversibility are now analysed. Starting with the MSD maps in Figure 12, one can observe different trends from what was obtained with the use of Macaulay brackets. First, deviations increase up to 300%. Second, there is no clear trend when varying $\rho_d$.
or $\mu_d$, but a region where higher values of MSD are observed. Having in mind the problems evidenced in the previous subsection for the parabolic formulation combined with the local history variable, the results are analyzed in this perspective.

<figure><img src="image_19.png" /><figcaption>**Figure 12.** Damage MSD map and reaction MSD map for the hyperbolic formulation using local history variable.</figcaption></figure>

In Figure 13, the negative part of the damage velocity $\langle \dot{d} \rangle_-$ is plotted as a function of $x$ and $t$ for different $\mu_d$ values and $\rho_d = 1.2 e^{-5}$ (line III). For the three lower values of $\mu_d$, it is clearly observed that irreversibility is not assessed. Going back to the MSD maps in Figure 12, this violation of the irreversibility constraint is associated with higher deviation from the elliptic solution. For the highest value of $\mu_d$ considered in Figure 13, irreversibility is recovered and the deviation is much lower and similar to what was obtained in Figure 8 when Macaulay brackets are used to impose damage irreversibility. Another assessment of damage irreversibility constraint is presented in Figure 14. This figure allows for a global assessment of irreversibility in the $(\rho_d, \mu_d)$ plane by presenting a map of the normalized space-time cumulative negative damage rate. It is clearly observed that except for $\rho_d \ge 1 e^{-3}$ kg/m or $\mu_d \ge 10$ Pa.s, irreversibility is not ensured using the local history variable. As for the parabolic formulation, these observations lead to the conclusion that the use of a local history variable for constraining damage to an irreversible evolution is not a robust strategy. Further, we observed that, as for the parabolic formulation, the main effect for the violation of irreversibility is a numerical instability in the explicit time stepping schemes. For the sake of conciseness, we do not report a detailed analysis of the reaction evolution and the damage profile as many of the results are corrupted.
<figure><img src="image_20.png" /><figcaption>**Figure 13.** $\langle \dot{d} \rangle_-$ maps in the $(x, t)$ plane of four $\mu_d$ values and $\rho_d = 1.2 \ e^{-5} \ \text{kg/m}$ (line III). The color scale is adapted to each maps $\dot{d}_{min} = -3 \ e^5, -2.5 \ e^5 - 5 \ e^3 - 4 \ e^{-3} \ \text{s}^{-1}$ for 13(a), 13(b), 13(c), 13(d) respectively.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 14.** Normalized space-time cumulative negative damage rate for the hyperbolic formulation using local history variable.</figcaption></figure>

## 4.5. Conclusion

The detailed analysis of the results obtained for the simulation of the dynamic failure of bar helps to draw some conclusions on several aspects:

(1) Irreversibility constraint:

* for both parabolic and hyperbolic formulation, it is mandatory that this condition is satisfied, otherwise numerical instabilities occur;

* using Macaulay brackets evals to be an efficient strategy to ensure damage irreversibility;

* the constraint is not completely guaranteed with the local history variable proposed in [14]. It has been shown for the parabolic and the hyperbolic damage formulation that with the local history variable, the robustness of the irreversibility constraint is parameter-dependent.

(2) Resolution strategy:

* using the Macaulay brackets to prevent damage healing, the results show that the critical time step defined in Section 3.4 is appropriate;

* in the context of explicit dynamic analysis of the displacement evolution, using an explicit damage resolution strategy is a good alternative to the classical explicit/implicit strategy. Indeed, the elliptic formulation predictions are recovered for low $\mu_d$ (and $\rho_d$) values. And, even if using damage formulation additional parameters leading to small deviation from the elliptic formulation requires a time step lower than the critical time step of the displacement problem, the corresponding additional cost is much lower than the linear system resolution required for the elliptic formulation and the explicit/implicit resolution;

* in practice, if keeping the time step to its value required for the displacement problem is required, the damage additional parameters are to be set carefully and
one will not benefit for the extended modelling offered by parabolic and hyperbolic
damage PDEs.

(3) Physical interpretation:

* the dissipative quadratic potential included in the parabolic and hyperbolic formulation induces a delay effect in the damage process but not in damage nucleation as the reaction plateau value does not depend on $\mu_d$ but the time to failure does. In other words, and as it could have been anticipated, this additional term limits damage rate. However, it has no impact on one could define as an apparent strength;

* the damage kinetic energy contribution of the hyperbolic formulation produces a behaviour typical of a wave propagation solution for damage evolution. $\rho_d$ controls the induced inertia effects. In the considered 1D example, oscillatory damage patterns are obtained when $\rho_d$ is set to increasing values. This effect can be observed as soon as the damage wave propagation speed is lower than that of the stress wave. In the opposite case, damage evolution is mainly governed by stress wave and the solutions tends to that of the elliptic formulation. There is a special case where the wave speed of damage and displacement are similar. In this case strong interaction between the stress and damage evolution are observed. Another consequence of an increasing damage inertia is the increase of stress level required for damage nucleation. Damage inertia increases the apparent strength of the material.

When parabolic or hyperbolic PDEs are used to govern damage evolution, the use of a fully explicit resolution strategy opens up new opportunities for fast failure numerical simulation using the phase field approach: efficient computational implementation, robust irreversibility constraint, enriched damage evolution phenomenology. In the next section, a more complex situation will be analyzed to assess the behaviour of parabolic and hyperbolic PDEs in a configuration that is more representative.

## 5. Kalthoff & Winkler test

The purpose of this section is to investigate the sensitivity of the crack path and the crack tip velocity upon the PDEs additional parameters in order to enrich the comparison of the predictions of the three PDEs. The Kalthoff et al. experiment offers the possibility to study these crack features [63]. An algorithm was developed to define the position of the crack tip from the damage field computed at each time step in order to provide for the crack path and the crack tip speed in a post-processing step. These are the most direct interpretable features that can be extracted from the simulations. The algorithm used to track the crack tip is presented in Appendix C.

The geometry, boundary conditions and the loading are given by Figure 15 and Eq. (58) respectively. The material properties are similar to the previous test and defined in Table 2 except for the Poisson's ratio that is set to 0.3. Half of the specimen is simulated due to symmetry and it is meshed with 994000 quadrangular elements with a linear interpolation. The parallel computing supported by FEniCS-dolfinx allows us to run our simulations on a cluster using 16 processors for around a 2-hour cost per simulation.
<figure><img src="image_22.png" /><figcaption>Figure 15. Kalthoff et al. experiment: model and geometry</figcaption></figure>

$$u_{\text{applied}} = \begin{cases} 0.5 \frac{V_0}{t_0} t^2 & t < t_0 = 10^{-6} \text{s}, V_0 = 16.5 \text{ m/s} \\ 0.5 V_0^2 + V_0(t - t_0) & \text{otherwise} \end{cases} \qquad (58)$$

5.1. Parabolic damage formulation results.

We first start investigating the prediction of the parabolic damage formulation. $\mu_d$ is varied between 110. Pa.s to 35.0 Pa.s. The Macaulay bracket are used to impose irreversibility to the damage evolution.

In Figure 16(a), the crack paths for different $\mu_d$ are represented at the same time instant and compared to the reference solution obtained with the elliptic damage PDE. A crack nucleates from the initial notch and then propagates with an angle of about 70°. This is consistent with numerous results in the literature. The zoom-in view clearly shows a delay on the crack advance with an increasing viscosity value and a slight deviation of the crack trajectories. The delay means that after the same simulated time, increasing the value of $\mu_d$ leads to shorter crack propagation distances.

This can also be demonstrated through the crack tip velocity evolution in Figure 16(b). In these plots, one can observe that the crack starts propagating at 0.02 ms. Its speed reaches a plateau at around 0.03 ms. Note that at the end of the elliptic simulation the crack tip speed vanishes. This is an artifact due to the crack tip detection algorithm that fails tracking after full failure of the plate. It is clear from Figure 16(b) that increasing the value of $\mu_d$ decreases the maximum crack velocity. However, the rising time of the velocity up to the plateau is kept constant varying $\mu_d$. This time is indeed a consequence of stress wave reflection in the specimen, it is thus independent from the damage modelling. One can argue that the lower value of maximum crack speed is a consequence of the lower velocity predicted with the parabolic formulation just after crack initiation and a constant duration (depending upon stress wave propagation only) of the crack
acceleration period.

These results confirm that decreasing $\mu_d$ leads to predictions that are consistent with those of the elliptic formulation. They also confirm that increasing $\mu_d$ produces a delay effect by decreasing the damage rate and consequently the crack speed.

<figure><img src="image_23.png" /><figcaption>**Figure 16.** Crack tip kinematics extracted at the same final time step. $c_R$, the Rayleigh wave speed</figcaption></figure>

## 5.2. Hyperbolic damage formulation results.

The sensitivity of the hyperbolic damage formulation prediction to the damage viscosity, $\mu_d$, and to the damage micro-inertia, $\rho_d$, is now analyzed. In Figures 17–18, the $\mu_d$-sensitivity of the crack path and velocity are plotted for six $\rho_d$ values: 1.0 $e^{-5}$, 6.458 $e^{-5}$, 1.077 $e^{-4}$, 5.0 $e^{-4}$, 9.6 $e^{-4}$, 5 $e^{-3}$ kg/m. The crack paths are plotted for the same time instant.

Observing all these results together, one can conclude that the influence of $\mu_d$ is the same as for the parabolic formulation. However, the higher the value of $\rho_d$ the lower the influence of $\mu_d$. For $\rho_d = 9.6 e^{-4}$ (Figure 18(d)), the curves representing the prediction for the 6 different values of $\mu_d$ collapse. Inertia effects are becoming dominant upon damping effects when $\rho_d$ is increased above $9.6 e^{-4}$.

Increasing the value of $\rho_d$ decreases the maximum crack speed and the crack length until preventing the crack to develop. It can be related to an inertial effect. If $\rho_d$ is higher than $1.077 e^{-4}$ kg/m (from Figures 17(e) and 17(f)), a crack arrest is obtained before the end of the simulation time. For ($\rho_d, \mu_d$) equals to ($1.077 e^{-4}$ kg/m, 30.0 Pa.s) (Figure 17(f)), and $\rho_d = 5.0 e^{-4}$ kg/m (Figure 17(b)), a jerky crack propagation is obtained. Then when $\rho_d$ is even increased, the crack definitively stops after a shorter and shorter time. From a global point of view, this means that for the same amount of energy provided by the external load, damage inertia increases the energy
dissipated in the damage process. Inertia can thus be seen as a toughening mechanism.

Fundamentally, damping penalises velocity while inertia penalises acceleration. This is observed
in our results through the three following points:

* when $\rho_d$ is increased, the crack nucleation time progressively increases: more energy is required to start propagation (to produce a non zero acceleration). This also refers to the strengthening mechanisms illustrated in the previous section;

* the transition between the acceleration phase and the constant velocity plateau becomes smoother and smoother when $\rho_d$ is increased;

* the crack initiation angle is decreased when $\rho_d$ is increased.

For a better assessment of this comments, Figure 19 presents the results for $\mu_d$ = 19 Pa.s and different values of $\rho_d$. The results for the highest value of $\rho_d$ is however a bit difficult to interpret, the propagation distance for this case being very small.
<figure><img src="image_24.png" /><figcaption>**Figure 17.** Crack path and crack tip velocity for $\rho_d = 1.0 \text{ e}^{-5}, 6.458 \text{ e}^{-5}, 1.077 \text{ e}^{-4} \text{ kg/m}$.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 18.** Crack path and crack tip velocity for $\rho_d = 5.0 \text{ e}^{-4}, 9.6 \text{ e}^{-4}, 5 \text{ e}^{-3} \text{ kg/m}$.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Figure 19. Crack path and crack tip velocity for $\mu_d$ = 19 Pa.s.</figcaption></figure>

**6. Discussion and Conclusion**

In the context of dynamic damage modelling through phase-field approaches, we propose a joint assessment of existing phase field approaches for the transient regime. The comparative study focuses on three aspects: (a) the resolution strategy in a fully explicit framework, (b) the irreversibility condition and (c) the physical interpretation.

The main focus is on fully explicit time integration scheme: both displacement and damage problem resolutions are performed with an explicit solver. While this approach is already deployed for models based on a parabolic damage formulation, its use with an hyperbolic damage PDEs is less common. Using an explicit time stepping scheme allows for considering condensed operators, operations are thus reduced to vector-vector operations, which improves efficiency and reduces time cost per time step compared to implicit methods. However, the stability analysis showed a dependence of the critical time step upon the PDEs parameters, in particular on the damage PDEs, that tends to reduce its value compared to the mechanical time step. Therefore, if one aims to provide an alternative to the classical explicit/implicit resolution with the elliptic damage formulation for large scale structures, the choice which is sometimes done to use the mechanical critical time step as a global one needs to be jeopardized. We provide for analytical expressions of the critical time step for the parabolic and hyperbolic formulation as well as the definition of additional parameters ($\rho_d, \mu_d$) ranges within which the mechanical time step can be used.

The challenges of the irreversibility constraint is also investigated in this paper. In particular, the lack of robustness in the local history variable [14], within an explicit framework is exposed for both the parabolic and the hyperbolic formulation. To circumvent this problem, we propose an algorithmic strategy for these two formulations, that allows to ensure irreversibility of the damage evolution at a very low computational cost. In the context of explicit dynamic analysis, these fully explicit resolution strategies open new perspectives for a computationally efficient and robust implementation of phase field damage models.
Last, deviations from a reference solution are evaluated using a metric-based approach. While it is not a reference, in the sense that it cannot be considered as the ground truth, the results from an elliptic damage PDE are used to provide insights for additional parameters selection. This metric-based comparison could be particularly valuable with an experimental reference, laying the groundwork for realistic model calibration. The phenomenology of the prediction of the three damage formulations are compared based on the results obtained for two classical examples. Our analysis reveals a delay effect controlled by the damage viscosity parameter for both parabolic and hyperbolic damage PDEs, with the hyperbolic formulation additionally capturing inertial effects. In addition, an inertial effect is captured with the hyperbolic damage PDEs. These observations highlight the ability of parabolic and hyperbolic damage models to effectively reproduce and control the dynamics of damage evolution, dynamics capabilities that the elliptic formulation tends to overestimates when the transient regime is significant.

## Conflicts of interest

The authors declare no competing financial interest.

## Dedication

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.

## Appendix A. Stability analysis of parabolic damage PDE

Recalling, the matrix form of the parabolic formulation,

$$ C_d \dot{\mathbf{D}} + \mathbf{K}_d \mathbf{D} = \mathbf{F}. \qquad (59) $$

for one degree of freedom (DoF), the matrix form is approximated by

$$ \left\{ \begin{array}{l} c = \mu_d h \\ k = \frac{G_c \ell_c}{h} \end{array} \right. \qquad (60) $$

The frequency is defined by the ratio between the coefficients Eq. (60) as

$$ \omega = \frac{k}{c} = \frac{G_c \ell_c}{\mu_d h^2} \qquad (61) $$

Let's consider two DoFs. Using the Euler explicit scheme

$$ \dot{\mathbf{D}}_{n+1} = \frac{1}{\Delta t} (\mathbf{D}_{n+1} - \mathbf{D}_n), \qquad (62) $$

Eq. (59) can be modified to introduce the amplification matrix

$$ \mathbf{D}_{n+1} = \mathbf{A} \mathbf{D}_n + \Delta t \mathbf{C}_d^{-1} \mathbf{F} \quad \text{with} \quad \mathbf{A} = (\mathbf{I} - \Delta t \mathbf{C}_d^{-1} \mathbf{K}_d) \qquad (63) $$

To ensure the stability of the time integration scheme, the norm of the spectral radius $\rho(\mathbf{A})$ must be strictly lower than 1: (Eq. (64)).

$$ |\rho(\mathbf{A})| < 1 \qquad (64) $$

In the case of a 1D bar element, a condensed damping matrix (row-sum method) and a stiffness matrix is expressed as:

$$ C_d = \frac{h \mu_d}{2} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \quad \text{and} \quad K_d = \frac{G_c \ell_c}{h} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} + \frac{h}{6} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}. \qquad (65) $$
The amplification matrix is now

$$A = \begin{bmatrix} 1 - \frac{2}{\mu_d h} \left( \frac{G_c \ell_c}{h} + \frac{h}{3} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \right) \Delta t & - \frac{2}{\mu_d h} \left( -\frac{G_c \ell_c}{h} + \frac{h}{6} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \right) \Delta t \\ - \frac{2}{\mu_d h} \left( -\frac{G_c \ell_c}{h} + \frac{h}{6} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \right) \Delta t & 1 - \frac{2}{\mu_d h} \left( \frac{G_c \ell_c}{h} + \frac{h}{3} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \right) \Delta t \end{bmatrix} \quad (66)$$

and its eigenspectrum is given by the determinant

$$\det(A - \lambda I) = \left[1 - \lambda - \frac{2\Delta t}{\mu_d h} \left(\frac{G_c \ell_c}{h} + \frac{h}{3} \left(2\Psi_+ + \frac{G_c}{l_c}\right)\right)\right]^2 - \left[\frac{2\Delta t}{\mu_d h} \left(-\frac{G_c \ell_c}{h} + \frac{h}{6} \left(2\Psi_+ + \frac{G_c}{l_c}\right)\right) - \lambda\right]^2 \\ = \left[1 - \lambda - \frac{4\Delta t}{\mu_d h^2} \left(G_c \ell_c + \frac{h^2}{12} \left(2\Psi_+ + \frac{G_c}{l_c}\right)\right)\right] \left[1 - \lambda - \frac{\Delta t}{\mu_d} \left(2\Psi_+ + \frac{G_c}{l_c}\right)\right]. \quad (67)$$

The roots of the second order polynomial Eq. (67), which results from $\det(A - \lambda I) = 0$, provides the maximum of the eigenvalues,

$$|\rho(A)| = \left| 1 - \frac{4\Delta t}{\mu_d h^2} \left( G_c \ell_c + \frac{h^2}{12} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \right) \right| \quad (68)$$

The condition Eq. (64) leads to the critical time step

$$\Delta t_{cd}^{1D} = \frac{\mu_d h^2}{2G_c \ell_c + \frac{h^2}{6} \left(2\Psi_+ + \frac{G_c}{l_c}\right)} \quad (69)$$

This expression makes the critical time step dependent to the positive part of the elastic energy density $\Psi_+$. Thus, when damage $d$ tends to 1, the increase of $\Psi_+(\underline{u})$ leads to a decrease in the time step assessment.

To avoid drastic decrease of $\Delta t_{cd}$, the following approximation

$$\Delta t_{cd}^{1D} = \frac{\mu_d h^2}{2G_c \ell_c} (1 + \epsilon(h^2)) \quad (70)$$

is considered.

Finally, the critical time step required to satisfy the stability of the Euler forward scheme for the parabolic damage PDE of the phase field model is

$$\Delta t_{cd}^{0D} = \frac{\mu_d h^2}{G_c \ell_c} \quad \text{and} \quad \Delta t_{cd}^{1D} = \frac{\mu_d h_{min}^2}{2G_c \ell_c} \quad (71)$$

## Appendix B. Stability analysis of hyperbolic damage PDE

Using the Verlet time integration scheme, the discrete form of Eq. (34) is such as:

$$M_d \ddot{D}_n + C_d \dot{D}_{n-1/2} + K_d D_n + F_d^i(U) = 0 \quad (72)$$

with a central difference scheme for the acceleration and an upwind scheme for the velocity:

$$\ddot{D}_n = \frac{D_{n+1} - 2D_n + D_{n-1}}{2\Delta t^2} \quad \text{and} \quad \dot{D}_{n-1/2} = \frac{D_n - D_{n-1}}{\Delta t} \quad (73)$$

Finally, the damage at $t_{n+1}$ results from the resolution of

$$M_d D_{n+1} = \Delta t^2 F_d(U_n) + (2M_d - \Delta t C_d - \Delta t^2 K_d) D_n + (\Delta t C_d - M_d) D_{n-1} \quad (74)$$

what leads to the following expression:

$$D_{n+1} = \Delta t^2 M_d^{-1} F_d(U_n) + (2I - \Delta t M_d^{-1} C_d - \Delta t^2 M_d^{-1} K_d) D_n + (\Delta t M_d^{-1} C_d - I) D_{n-1}. \quad (75)$$

Considering, $X_{n+1} = \begin{pmatrix} D_{n+1} \\ D_n \end{pmatrix}$, the numerical solution can be stated in the recurrent from

$$X_{n+1} = T X_n + P_{n+1} \quad (76)$$
The matrix **T** and the vector **P**<sub>n+1</sub> are the amplification matrix and the external force respectively

$$
\mathbf{T} = \left[ \begin{array}{c} 2\mathbf{I} - \Delta t \mathbf{M}_d^{-1} \mathbf{C}_d - \Delta t^2 \mathbf{M}_d^{-1} \mathbf{K}_d \quad \Delta t \mathbf{M}_d^{-1} \mathbf{C}_d - \mathbf{I} \\ \mathbf{I} \\ \mathbf{0} \end{array} \right], \mathbf{P}_{n+1} = \left( \begin{array}{c} \Delta t^2 \mathbf{M}_d^{-1} \mathbf{F}_d(\mathbf{U}_n) \\ \mathbf{0} \end{array} \right) \quad (77)
$$

Assuming linearly independent eigenpairs, a modal decomposition can be applied:

$$
\underline{\Phi}_i^T \mathbf{M}_d \underline{\Phi}_j^T = \delta_{ij}, \quad \underline{\Phi}_i^T \mathbf{C}_d \underline{\Phi}_j^T = 2\omega_i \xi_i \delta_{ij}, \quad \underline{\Phi}_i^T \mathbf{K}_d \underline{\Phi}_j^T = \omega_i^2 \delta_{ij}, \quad d = d_i \underline{\Phi}_i, \quad i \in \{1, \dots, N\}
$$

where $\xi_i$, $\omega_i$ and $\Phi$ represents the damping ratio, the natural frequencies and the eigenbasis vectors. $\underline{\Phi}_i^T \mathbf{M}_d \underline{\Phi}_j^T = \delta_{ij}$ represents a normalization of the system. Considering, the case of a 1D bar element, the modal decomposition gives the following expressions:

$$
\mathbf{M}_d^{-1} \mathbf{C}_d = \begin{bmatrix} 2\xi_1\omega_1 & 0 \\ 0 & 2\xi_2\omega_2 \end{bmatrix} \quad \text{and} \quad \mathbf{M}_d^{-1} \mathbf{K}_d = \begin{bmatrix} \omega_1^2 & 0 \\ 0 & \omega_2^2 \end{bmatrix}. \qquad (78)
$$

Thus, the amplification matrix can be written as

$$
T = \begin{bmatrix}
2 - 2\Delta t\omega_1\xi_1 - \Delta t^2\omega_1^2 & 0 & \Delta t 2\omega_1\xi_1 - 1 & 0 \\
0 & 2 - 2\Delta t\omega_2\xi_2 - \Delta t^2\omega_2^2 & 0 & \Delta t 2\omega_2\xi_2 - 1 \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0
\end{bmatrix}
\quad (79)
$$

As for the previous analysis on the parabolic PDE, the root ($\lambda_1, \dots, \lambda_n$) of the polynomial character-
istics give the eigenvalues of the amplification matrix **T**:

$$
\det(\mathbf{T} - \lambda \mathbf{I}) = 0 \implies \lambda^2 - 2I_1\lambda + I_2 = 0 \quad \text{with} \quad \begin{cases} 2 I_1 = \operatorname{tr}(\mathbf{T}) = \sum_{i=1}^2 (2 - 2\Delta t \xi_i \omega_i - \Delta t^2 \omega_i^2) \\ I_2 = \det(\mathbf{T}) = \prod_{i=1}^2 (1 - 2\Delta t \xi_i \omega_i) \end{cases} \quad (80)
$$

given the stability condition

$$
|\max(\lambda_1, \dots, \lambda_n)| < 1. \tag{81}
$$

These eigenvalues are

$$
\begin{cases}
\lambda_{1,2} = \frac{-\Delta t^2 \omega_1^2}{2} - \Delta t \omega_1 \xi_1 \pm \frac{\Delta t \omega_1}{2} \sqrt{(\Delta t \omega_1 + 2 \xi_1 - 2)(\Delta t \omega_1 + 2 \xi_1 + 2)} + 1 \\
\lambda_{3,4} = \frac{-\Delta t^2 \omega_2^2}{2} - \Delta t \omega_2 \xi_2 \pm \frac{\Delta t \omega_2}{2} \sqrt{(\Delta t \omega_2 + 2 \xi_2 - 2)(\Delta t \omega_2 + 2 \xi_2 + 2)} + 1
\end{cases}
\tag{82}
$$

Among the eigenvalues, $\lambda_1$ has the largest absolute value. The stability conditions thus leads to

$$
|\lambda_1| < 1 \rightarrow \left| \frac{-\Delta t^2 \omega_1^2}{2} - \Delta \omega_1 \xi_1 + \frac{\Delta t \omega_1}{2} \sqrt{(\Delta t \omega_1 + 2\xi_1)^2 - 4} + 1 \right| < 1 \quad (83)
$$

with

$$
0 < \Delta t^2 \omega_1^2 + 4 \Delta t \omega_1 \xi_1 - 4. \qquad (84)
$$

After some manipulations, one gets

$$
(\Delta t \omega_1)_\pm = 2(-\xi_1 \pm \sqrt{\xi_1^2 + 1}). \quad (85)
$$

Eventually, the critical time step $\Delta t_{cd}$ to satisfy the stability condition is an function of the damage damping and inertia:

$$
\Delta t_{cd} = \frac{2(\sqrt{\xi_1^2 + 1} - \xi_1)}{\omega_1} = f(\rho_d, \mu_d). \qquad (86)
$$

To specify, the expression of ($\xi_1, \omega_1$), a 1D bar analysis is considered with P1 elements. The matrices $\mathbf{M}_d^{-1}\mathbf{C}_d$ and $\mathbf{M}_d^{-1}\mathbf{K}_d$ are diagonalized and indexed by DZ indices,

$$
\mathbf{M}_d^{-1} \mathbf{C}_d = \frac{\mu_d}{3\rho_d} \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \rightarrow \mathbf{M}_d^{-1} \mathbf{C}_d^{\text{DZ}} = \frac{\mu_d}{3\rho_d} \begin{bmatrix} 3 & 0 \\ 0 & 1 \end{bmatrix} \quad (87)
$$
$$ 
\begin{aligned}
\mathbf{M}_d^{-1} \mathbf{K}_d &= \frac{2G_c \ell_c}{\rho_d h^2} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} + \frac{1}{3\rho_d} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix} \\
\rightarrow \mathbf{M}_d^{-1} \mathbf{K}_d^{DZ} &= \begin{bmatrix} \frac{4G_c \ell_c}{\rho_d h^2} + \frac{1}{3\rho_d} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) & 0 \\ 0 & \frac{1}{\rho_d} \left( 2\Psi_+ + \frac{G_c}{l_c} \right) \end{bmatrix}
\end{aligned}
\quad (88) $$

By identification with Eq. (78), the natural frequencies and the damping ratio are defined by

$$ \omega_1^2 = \frac{4G_c \ell_c}{\rho_d h^2} + \frac{1}{3\rho_d} \left( 2\Psi_+ + \frac{G_c}{l_c} \right), \quad \omega_2^2 = \frac{1}{\rho_d} \left( 2\Psi_+ + \frac{G_c}{l_c} \right), \quad (89) $$

$$ 2\xi_1\omega_1 = 3\frac{\mu_d}{3\rho_d}, \quad 2\xi_2\omega_2 = \frac{\mu_d}{3\rho_d} \quad (90) $$

Considering the Taylor expansion of $\omega_1^2$ and $\omega_2^2$ with respect to h,

$$ \omega_1^2 = 4 \frac{G_c \ell_c}{\rho_d h^2} + \epsilon(h^2), \quad \omega_2^2 = 0, \quad (91) $$

with finally

$$ \omega_1 = \frac{2}{h} \sqrt{\frac{G_c \ell_c}{\rho_d}}, \quad \xi_1 = \frac{\mu_d}{2\omega_1 \rho_d} \quad (92) $$

## Appendix C. Crack-tip tracking from a continuous damage field

A non-intrusive algorithm based on the clustering of the “damaged” nodes has been implemented to define successive positions of an assumed crack tip. A damage tolerance of 0.9 was obtained to consider a node as “cracked”, pictured as a colored empty circle in Figure 20. At each time increment, the set of “damaged” nodes is updated. Among these nodes, only ones downstream the crack tip are tracked, while “damaged” nodes developing upstream of the crack tip and laterally along the crack are ignored. A Principal Component Analysis (PCA) is then performed on the spatial distribution of the “cracked” nodes, this leading to main spatial directions (sketched as a grey orthogonal basis). The latter helps to approximate crack direction which is necessary for the “damaged” nodes clustering described above. A red box is plotted around the final set of “cracked” nodes and the barycentre of this set is identified as a new crack tip position. These barycentres along with their associated times, are used to construct a spline. The Python library *Splipy* provides the parametric equation of the curve and can compute its derivative analytically. Although the analytical expression cannot be explicitly provided, it is possible to evaluate it over a given time interval.
<figure><img src="image_27.png" /><figcaption>**Figure 20.** A tracking illustration of the crack advance</figcaption></figure>

**References**

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

[2] B. Bourdin, G. A. Francfort, J.-J. Marigo, "The Variational Approach to Fracture", *Journal of Elasticity* **91** (2008), no. 1, p. 5-148.

[3] L. Ambrosio, V. M. Tortorelli, "Approximation of functional depending on jumps by elliptic functional via t-convergence", *Communications on Pure and Applied Mathematics* **43** (1990), no. 8, p. 999-1036.

[4] C. J. Larsen, C. Ortner, E. Süli, "Existence of solutions to a regularized model of dynamic fracture", *Mathematical Models and Methods in Applied Sciences* **20** (2010), no. 07, p. 1021-1048.

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

[6] V. Ziaei-Rad, Y. Shen, "Massive parallelization of the phase field formulation for crack propagation with time adaptivity", *Computer Methods in Applied Mechanics and Engineering* **312** (2016), p. 224-253.

[7] T. Li, *Analyse de la rupture dynamique fragile via les modèles d'endommagement à gradient : principes variationnels et simulations numériques*, phdthesis, Université Paris-Saclay, 2016.

[8] H. Ren, X. Zhuang, C. Anitescu, T. Rabczuk, "An explicit phase field method for brittle dynamic fracture", *Computers & Structures* **217** (2019), p. 45-56.

[9] J. Bleyer, C. Roux-Langlois, J.-F. Molinari, "Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms", *International Journal of Fracture* **204** (2017), no. 1, p. 79-100.

[10] T. K. Mandal, V. P. Nguyen, J.-Y. Wu, "Evaluation of variational phase-field models for dynamic brittle fracture", *Engineering Fracture Mechanics* **235** (2020), p. 107169.

[11] V. Agrawal, K. Dayal, "Dependence of equilibrium Griffith surface energy on crack speed in phase-field models for fracture coupled to elastodynamics", *International Journal of Fracture* **207** (2017), p. 243-249.

[12] M. E. Gurtin, "Generalized Ginzburg-Landau and Cahn-Hilliard equations based on a microforce balance", *Physica D: Nonlinear Phenomena* **92** (1996), no. 3-4, p. 178-192.

[13] C. Kuhn, R. Müller, "A continuum phase field model for fracture", *Engineering Fracture Mechanics* **77** (2010), no. 18, p. 3625-3634.

[14] C. Miehe, M. Hofacker, F. Welschinger, "A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits", *Computer Methods in Applied Mechanics and Engineering* **199** (2010), no. 45-48, p. 2765-2778.

[15] A. Schlüter, A. Willenbücher, C. Kuhn, R. Müller, "Phase field approximation of dynamic brittle fracture", *Computational Mechanics* **54** (2014), no. 5, p. 1141-1161.
[16] C.-H. Chen, E. Bouchbinder, A. Karma, "Instability in dynamic fracture and the failure of the classical theory of cracks", *Nature Physics* **13** (2017), no. 12, p. 1186-1190.

[17] T. Wang, X. Ye, Z. Liu, D. Chu, Z. Zhuang, "Modeling the dynamic and quasi-static compression-shear failure of brittle materials by explicit phase field method", *Computational Mechanics* **64** (2019), no. 6, p. 1537-1556.

[18] J. Rudshaug, T. Børvik, O. S. Hopperstad, "Modeling brittle crack propagation for varying critical load levels: a dynamic phase-field approach", *International Journal of Fracture* **245** (2024), p. 1-17.

[19] C. Cattaneo, "Sur une forme de l'équation de la chaleur eliminant la paradoxe d'une propagation instantanée", *Comptes rendus de l'Académie des sciences* **247** (1958), p. 431-433.

[20] J. Masoliver, G. H. Weiss, "Finite-velocity diffusion", *European Journal of Physics* **17** (1996), no. 4, p. 190.

[21] P. Guillemet, J.-P. Bardon, "Conduction de la chaleur aux temps courts: les limites spatio-temporelles des modèles parabolique et hyperbolique", *International Journal of Thermal Sciences* **39** (2000), no. 9, p. 968-982.

[22] C. Donaghy-Spargo, "On Heaviside's contributions to transmission line theory: waves, diffusion and energy flux", *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **376** (2018), no. 2134, p. 20170457.

[23] M. Frémond, B. Nedjar, "Damage, gradient of damage and principle of virtual power", *International Journal of Solids and Structures* **33** (1996), p. 1083-1103.

[24] M. Frémond, B. Nedjar, "Damage and principle of virtual power", *Comptes rendus de l'Académie des sciences* **317** (1993), p. 857-864.

[25] A. F. Parrinello, *A Rate-Pressure-Dependent Thermodynamically-Consistent Phase Field Model for the Description of Failure Patterns in Dynamic Brittle Fracture*, Theses, University of Oxford, 2017.

[26] N. Favrie, A. Renaud, D. Kondo, "Hyperbolic modeling of gradient damage and one-dimensional finite volume simulations", *Computer Methods in Applied Mechanics and Engineering* **419** (2024), p. 116643.

[27] D. Kamensky, G. Moutsanidis, Y. Bazilevs, "Hyperbolic phase field modeling of brittle fracture: Part I—theory and simulations", *Journal of the Mechanics and Physics of Solids* **121** (2018), p. 81-98.

[28] R. Ortega, A. M. Robles-Pérez, "A Maximum Principle for Periodic Solutions of the Telegraph Equation", *Journal of Mathematical Analysis and Applications* **221** (1998), no. 2, p. 625-651.

[29] L. Svolos, H. M. Mourad, G. Manzini, K. Garikipati, "A fourth-order phase-field fracture model: Formulation and numerical solution using a continuous/discontinuous Galerkin method", *Journal of the Mechanics and Physics of Solids* **165** (2022), p. 104910.

[30] P. Farrell, C. Maurini, "Linear and nonlinear solvers for variational phase-field models of brittle fracture", *International Journal for Numerical Methods in Engineering* **109** (2017), no. 5, p. 648-667.

[31] T. Nguyen, J. Réthoré, M. Baietto, "Phase field modelling of anisotropic crack propagation", *European Journal of Mechanics - A/Solids* **65** (2017), p. 279-288.

[32] G. Molnár, A. Doitrand, A. Jaccon, B. Prabel, A. Gravouil, "Thermodynamically consistent linear-gradient damage model in Abaqus", *Engineering Fracture Mechanics* **266** (2022), p. 108390.

[33] A. Mielke, "Evolution of rate-independent systems", in *Handbook of differential equations: Evolutionary equations*, vol. 2, North-Holland : Elsevier, 2005, p. 461-559.

[34] K. Pham, J.-J. Marigo, "Approche variationnelle de l'endommagement: II. Les modèles à gradient", *Comptes Rendus Mécanique* **338** (2010), no. 4, p. 199-206.

[35] J. C. Simo, J. Ju, "Strain-and stress-based continuum damage models—I. Formulation", *International Journal of Solids and Structures* **23** (1987), no. 7, p. 821-840.

[36] M. Cervera, "Viscoelasticity and rate-dependent continuum damage models", *CIMNE, Monography N°-79, Technical University of Catalunya, Barcelona, Spain* (2003), p. 0-73.

[37] J. Lemaître, R. Desmorat, M. Sauzay, "Anisotropic damage law of evolution", *European Journal of Mechanics - A/Solids* **19** (2000), no. 2, p. 187-208.

[38] G. Lancioni, G. Royer-Carfagni, "The variational approach to fracture mechanics. A practical application to the French Panthéon in Paris", *Journal of elasticity* **95** (2009), p. 1-30.

[39] H. Amor, J.-J. Marigo, C. Maurini, "Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments", *Journal of the Mechanics and Physics of Solids* **57** (2009), no. 8, p. 1209-1229.

[40] G. Pijaudier-Cabot, R. De Borst, J. Mazars, "Continuous damage models for fracture of concrete", in *Proceedings of the International Conference on Fracture Mechanics of Concrete and Concrete Structures (FraMCoS-4)* (Cachan, France) (L. Swets & Zeitlinger, ed.), vol. 2, 2001.

[41] J.-Y. Wu, M. Cervera, "A novel positive/negative projection in energy norm for the damage modeling of quasi-brittle solids", *International Journal of Solids and Structures* **139** (2018), p. 250-269.

[42] F. Freddi, G. Royer-Carfagni, "Regularized variational theories of fracture: A unified approach", *Journal of the Mechanics and Physics of Solids* **58** (2010), no. 8, p. 1154-1174.

[43] R. Desmorat, "Dissymétrie de comportement élastique anisotrope couplé ou non à l'endommagement", *Comptes Rendus de l'Académie des Sciences-Series IIB-Mechanics* **328** (2000), no. 6, p. 445-450.
[44] M. François, "A damage model based on Kelvin eigentensors and Curie principle", *Mechanics of Materials* **44** (2012), p. 23-24.

[45] Q.-C. He, Q. Shao, "Closed-Form Coordinate-Free Decompositions of the Two-Dimensional Strain and Stress for Modeling Tension-Compression Dissymmetry", *Journal of Applied Mechanics* **86** (2019), no. 3, p. 031007.

[46] T.-T. Nguyen, J. Yvonnet, D. Waldmann, Q.-C. He, "Implementation of a new strain split to model unilateral contact within the phase field method", *International Journal for Numerical Methods in Engineering* **121** (2020), no. 21, p. 4717-4733.

[47] N. van Dijk, J. Espadas-Escalante, P. Isaksson, "Strain energy density decompositions in phase-field fracture theories for orthotropy and anisotropy", *International Journal of Solids and Structures* **196-197** (2020), p. 140-153.

[48] V. Ziaei-Rad, M. Mollaali, T. Nagel, O. Kolditz, K. Yoshioka, "Orthogonal decomposition of anisotropic constitutive models for the phase field approach to fracture", *Journal of the Mechanics and Physics of Solids* **171** (2023), p. 105143.

[49] V. B. Thanh, "Phase field method with strain orthogonal decomposition for modelling of damage in heterogeneous materials obtained by X-ray computed tomography images", *Transport and Communications Science Journal* **74** (2023), p. 20-34.

[50] B. Bourdin, *Une Méthode variationnelle en mécanique de la rupture*, PhD Thesis, Université Paris 13, 1998, 162 p. pages.

[51] T. Schlick, *Molecular modeling and simulation: an interdisciplinary guide*, vol. 2, New-York : Springer, 2010.

[52] M. Géradin, D. J. Rixen, *Mechanical vibrations: theory and application to structural dynamics*, Chichester ; New York ; Weinheim : John Wiley & Sons, 2014.

[53] R. Kolman, J. González, J. Kopačka, M. Mračko, R. Cimrman, S. S. Cho, A. Tkachuk, J. Plešek, M. Okrouhlík, D. Gabriel, I. Nčmece, I. Ševčíke K.C. Park, "Explicit time integration in finite element method for structural dynamic, wave propagation and contact-impact problems: A recent progress", Programs and algorithms of numerical mathematics (PANM), Institute of Mathematics of the Czech Academy of Sciences, 2020.

[54] L. De Lorenzis, T. Gerasimov, "Modeling in Engineering Using Innovative Numerical Methods for Solids and Fluids", in *Numerical implementation of phase-field models of brittle fracture*, CISM, International Centre for Mechanical Sciences, Udine, vol. 599, Springer, 2020, p. 75-101.

[55] G. Molnár, A. Gravouil, R. Seghir, J. Réthoré, "An open-source Abaqus implementation of the phase-field method to study the effect of plasticity on the instantaneous fracture toughness in dynamic crack propagation", *Computer Methods in Applied Mechanics and Engineering* **365** (2020), p. 113004.

[56] N. Boumediene, "Modélisation de l'endommagement avec gradient: Application à la dynamique rapide", *Bulletin des Laboratoires des Ponts et Chaussées* (2005), no. 256-57, p. 243-249.

[57] C. Bilgen, K. Weinberg, "On the crack-driving force of phase-field models in linearized and finite elasticity", *Computer Methods in Applied Mechanics and Engineering* **353** (2019), p. 348-372.

[58] L. Hai, J.-Y. Wu, J. Li, "A phase-field damage model with micro inertia effect for the dynamic fracture of quasi-brittle solids", *Engineering Fracture Mechanics* **225** (2020), p. 106821.

[59] D. Serfőzó, B. Pere, "A method to accurately define arbitrary algorithmic damping character as viscous damping", *Archive of Applied Mechanics* **93** (2023), no. 9, p. 3581-3595.

[60] J. Dokken, I. Baratta, J. Dean, M. Habera, J. Hale, C. Richardson, M. Rognes, N. Sime, G. Wells, "DOLFINx: The next generation FEniCS problem solving environment", *Zenodo* (2023).

[61] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener et al., "PETSc Users Manual", Tech. report, Argonne National Laboratory, Illinois, 2022.

[62] T. Helfer, B. Michel, J.-M. Proix, M. Salvo, J. Sercombe, M. Casella, "Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform", *Computers & Mathematics with Applications* **70** (2015), no. 5, p. 994-1023.

[63] J. F. Kalthoff, "Modes of dynamic shear failure in solids", *International Journal of fracture* **101** (2000), no. 1, p. 1-31.