<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRMECA,Unicode,biblatex,published]{cedram}
\addbibresource{crmeca20240702.bib}
\newenvironment{noXML}{}{}%
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\stbfrac#1#2{({#1})/({#2})}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{\noalign{\relax}\hline}
\def\botline{\\\hline}
\def\xsection{}
\allowdisplaybreaks
\def\ndash{\text{--}}
\def\inlinefig#1{\includegraphics{#1}}
\usepackage{etoolbox}
\usepackage{upgreek}
\usepackage{upgreek}
\usepackage[T1]{fontenc}
\def\jobid{crmeca20240702}
\usepackage{multirow} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\xmorerows#1#2{#2}
%\graphicspath{{/tmp/\jobid_figs/web/}} 
\graphicspath{{./figures/}} 
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\def\refinput#1{}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\csdef{Seqnsplit}{\\}
\def\Tg#1#2#3#4{}
\def\refinput#1{}
\def\back#1{}
\usepackage{algorithm}
\usepackage{algpseudocode}
\def\beqimg#1#2{}
\def\eeqimg#1#2{}
\let\noTagSeqn\relax
\DOI{10.5802/crmeca.297}
\datereceived{2024-08-28}
\daterevised{2025-01-03}
\datererevised{2025-04-07}
\dateaccepted{2025-04-22}
\ItHasTeXPublished
}
{\documentclass[crmeca]{article}
\usepackage{upgreek}
\usepackage[T1]{fontenc} 
\usepackage{stmaryrd}
\def\CDRdoi{10.5802/crmeca.297}
\def\0{\phantom{0}}
\def\newline{\break}
\def\centering{}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\def\newline{\unskip\break}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\def\stbfrac#1#2{({#1})/({#2})}
\let\beqimg\relax
\let\eeqimg\relax
}
\makeatother

\usepackage{moreverb}
\usepackage{mathrsfs}

\dateposted{2025-05-21}
\begin{document}

\begin{noXML}

\CDRsetmeta{articletype}{research-article}

\shortrunauthors

\title{Dynamic damage modelling through phase-field  approaches:
assessment, critical analysis and comparison}

\alttitle{Mod\'elisation dynamique de l'endommagement par
l'approche champ de phase : \'evaluation, analyse critique et comparaison}

\author{\firstname{Lamia} \lastname{Mersel}\CDRorcid{0009-0000-8894-1320}\IsCorresp}
\address{DMAS, ONERA, Lille, 59000, France}
\address{Nantes Universit\'e, Ecole Centrale Nantes, 
CNRS, GeM, UMR 6183, 1 rue de la No\"e, Nantes, 44321, France}
\email[L. Mersel]{lamia.mersel@ec-nantes.fr}

\author{\firstname{Pascal} \lastname{Bouda}\CDRorcid{0009-0006-2905-7504}}
\address{Universit\'e Paris-Saclay, CEA, Service d'\'Etudes M\'ecaniques et Thermiques, 
91191 Gif-sur-Yvette, France}
\email[P. Bouda]{pascal.bouda@cea.fr}

\author{\firstname{J\'er\'emy} \lastname{Germain}\CDRorcid{0009-0001-7033-2356}}
\addressSameAs{1}{DMAS, ONERA, Lille, 59000, France}
\email[J. Germain]{jeremy.germain@onera.fr}

\author{\firstname{Julien}\nobreakauthor\lastname{R\'ethor\'e}\CDRorcid{0000-0003-0690-8549}}
\addressSameAs{2}{Nantes Universit\'e, Ecole Centrale Nantes, 
CNRS, GeM, UMR 6183, 1 rue de la No\"e, Nantes, 44321, France}
\email[J. R\'ethor\'e]{julien.rethore@ec-nantes.fr}

\keywords{\kwd{Dynamic fracture} 
\kwd{Damage phase field formulation}
\kwd{Explicit time integration} 
\kwd{Sensitivity analysis}}

\altkeywords{\kwd{Rupture dynamique}
\kwd{Formulation d'endommagement}
\kwd{Int\'egration temporelle explicite}
\kwd{Analyse de sensibilit\'e}}

\begin{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. 
\end{abstract}

\begin{altabstract}
Ces travaux \'evaluent l'aptitude de diff\'erents mod\`eles \`a champs
de phase  pour mod\'eliser l'\'evolution dynamique de l'endommagement
jusqu'\`a rupture.  L'\'evolution de l'endommagement peut \^etre
d\'ecrite   par une \'equation \`a d\'eriv\'ee partielle (EDP) de type
elliptique, parabolique  ou hyperbolique. Ces derni\`eres permettent
l'utilisation de sch\'emas  d'int\'egration explicites, strat\'egie
majoritairement employ\'ee dans les  cas d'\'etude de dynamique
transitoire. La condition d'irr\'eversibilit\'e de l'endommagement   a
\'egalement \'et\'e impos\'ee de diff\'erente mani\`eres. Il est
montr\'e que l'utilisation  d'une variable d'histoire de la densit\'e
d'\'energie \'elastique positive ne permet  pas de garantir.
strictement l'irr\'eversibilit\'e et peut  induire des instabilit\'es. 
N\'eanmoins, si cette strat\'egie est combin\'ee \`a l'utilisation d'un
terme de viscosit\'e  d'endommagement dans l'EPD gouvernant son
\'evolution, il est possible  d'assurer l'irr\'eversibilit\'e.
Autrement, le cadre explicite de la strat\'egie  de r\'esolution offre
la possibilit\'e d'imposer l'irr\'eversibilit\'e de fa\c{c}on 
algorithmique. Ensuite, au travers des cas test d'une barre en traction
et de   l'exp\'erience d'impact de Kalthoff et al., des \'etudes
param\'etriques ont \'et\'e  effectu\'ees pour \'evaluer la
sensibilit\'e aux param\`etres des EPDs. Pour finir,  une \'etude
comparative des solutions obtenues par les mod\`eles de champs  de
phase en dynamique a \'et\'e r\'ealis\'ee via l'utilisation de
m\'etrique bas\'ee   sur des quantit\'es physiques.
\end{altabstract}

\maketitle

\end{noXML}

\newcommand*{\fig}[1]{Figure~\ref{#1}}
\newcommand*{\tab}[1]{Table~\ref{#1}}

\newcommand{\me}{\ensuremath{\mathrm{e}}}
\newcommand*{\eq}[1]{Equation~(\ref{#1})}

\newcommand*{\twobar}{=}
\newcommand*{\rttensor}[1]{\underset{\twobar}{#1}}

\newcommand*{\alg}[1]{Algorithm~\ref{#1}}
\def\d{{\mathrm d}}
\newcommand{\lc}{\ell_{c}}

\def\say#1{``#1''}

\section{Introduction}\label{sec0}
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~\cite{francfort_revisiting_1998, bourdin_variational_2008},
which is also  known as the Ambrosio--Tortorelli 
approximation~\cite{ambrosio1990approximation}. 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 $\lc$ that governs the formation of  a continuous
damage growing band.  When $\lc$ 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~\cite{larsen2010existence, borden_phase-field_2012, 
ziaei-rad_massive_2016, li_analyse_2016, ren_explicit_2019}. 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~\cite{borden_phase-field_2012,
bleyer2017dynamic, MANDAL2020}. 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~\cite{agrawal2017dependence}.  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~\cite{gurtin1996generalized},  is one of the most used in
the literature to model a dynamic  damage evolution~\cite{KUHN2010,
miehe_2010, schluter_phase_2014,  ziaei-rad_massive_2016,
chen2017instability, wang_modeling_2019,  rudshaug2024modeling}. 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~\cite{cattaneo1958forme}, as
an elliptic one can be considered to do.  Although this description is
considered physically inappropriate  in some scientific
fields~\cite{cattaneo1958forme, Jaume_Masoliver_1996,  GUILLEMET2000},
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~\cite{cattaneo1958forme}. 
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~\cite{donaghy2018heaviside}.   Since the
1990s, a less common hyperbolic PDE was also introduced for  damage
modelling~\cite{Frmond1996DamageGO, Frmond1993Damage, parrinello_2017,
FAVRIE2024}. 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~\cite{chen2017instability}.  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~\cite{parrinello_2017}.
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~\cite{kamensky2018}. Besides,  this expression can also be derived
from a boundedness solution analysis  of a generic hyperbolic
equation~\cite{ORTEGA1998625}.

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~\cite{li_analyse_2016,  rudshaug2024modeling}. 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~\ref{section:math_background}
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~\ref{section:numerical_strategy}. 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~\ref{section:bar} and~\ref{section:kalthoff}.  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\break compared.

\section{Theory---dynamic phase field for fracture
modelling}\label{section:math_background} %2
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. 

\subsection{General framework}\label{sec2.1}
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~\cite{Frmond1993Damage,
gurtin1996generalized}.

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~\cite{schluter_phase_2014, li_analyse_2016, chen2017instability, 
svolos2022fourth, FAVRIE2024}. This is done setting the Gateaux
derivative of the action functional  $S$:
{\begin{equation}\label{Action}
S[T] =  
\int_{T} L^{\mathrm{dyn}}(t, d, \underline{u},
\underline{\nabla} d, 
\rttensor{\nabla u},
\dot{d}, \dot{\underline{u}}) \, \d t
\end{equation}}\unskip
to zero, where $L^{\mathrm{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, $\dot{d}$ its time derivative and $\nabla$ the gradient
operator. In its general form, the Lagrangian $L^{\mathrm{dyn}}$ is
defined as follows: 
{\begin{equation}\label{Lagrangian_detaille}
L^{\mathrm{dyn}}(t, \underline{u}, \dot{\underline{u}}, d, \dot{d})  =
\mathcal{K}_u+ \mathcal{K}_d- 
( \mathcal{E}_{el} + \mathcal{E}_{\mathrm{fr}}  - 
\mathcal{W}_{\mathrm{ext}}(\underline{u})).
\end{equation}}\unskip
$\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}_{\mathrm{fr}}$ the fracture
energy and $\mathcal{W}_{\mathrm{ext}}$ the work of external forces to
be defined.

Over $\Omega$, the considered bounded open subset of $R^{\mathrm{dim}}$
($\mathrm{dim}$ being the dimension of the problem), the kinetic
energies are defined as
{\begin{eqnarray}\label{K}
\mathcal{K}_u 
&amp;=&amp;
\int_{\Omega} 
\frac{\rho_u}{2} \dot{\underline{u}}^2  
\d\Omega,
\Seqnsplit
\mathcal{K}_d 
&amp;=&amp; 
\int_{\Omega}  \frac{\rho_d}{2} \dot{d}^2  \d\Omega.
\end{eqnarray}}\unskip
In these equations, $\rho_u$ is the mass density and $\rho_d$ is a
so-called micro-damage inertia. The elastic strain energy is
{\begin{eqnarray}\label{E_pot}
\mathcal{E}_{el}
&amp;=&amp; 
\int_{\Omega} 
\Psi (\rttensor{\varepsilon}( \underline{u}), d  )  \d\Omega
\end{eqnarray}}\unskip
with $\Psi $ the strain energy density which can be written as:
{\begin{equation}
\Psi 
(\rttensor{\varepsilon}( \underline{u}), d  ) =
\tfrac{1}{2} g(d) \rttensor{\varepsilon} : \mathbb{C}:  \rttensor{\varepsilon},
\end{equation}}\unskip
$\mathbb{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:
{\begin{equation}
\rttensor{\sigma}  = g(d) \mathbb{C} : \rttensor{\varepsilon}.
\end{equation}}\unskip

The fracture energy density
$\mathcal{E}_{\mathrm{fr}}
(\lc, d, 
\underline{\nabla}d)$ 
is as follows:
{\begin{equation}\label{E_pot+E_frac}
\mathcal{E}_{el} + \mathcal{E}_{\mathrm{fr}}
= 
\int_{\Omega} 
\Psi (\rttensor{\varepsilon}( \underline{u}), d  )  \d\Omega+ 
\int_{\Omega}  G_c \gamma
( \lc, d,  \underline{\nabla}  d  )  \d\Omega.
\end{equation}}\unskip

It is a regularized form of the original integral of the fracture 
dissipated energy used within small strain
framework~\cite{francfort_revisiting_1998}. This form introduces a
crack  density function $\gamma$ which depends on an internal length 
$\lc$, the damage field $d$ and its gradient $\underline{\nabla} d$.
Among the different form of this crack density function, we choose
{\begin{equation}\label{fracture_density}
\gamma(\lc, {d}, \underline{\nabla} d) =   
{\frac{1}{4 c_w} \left(\frac{w(d)}{\lc} +
\lc |\underline{\nabla} d |^2 \right)}=
\frac{1}{2} 
\left(\frac{d^2}{\lc} +
\lc |\underline{\nabla} d |^2 \right)
\end{equation}}\unskip
proposed in~\cite{farrell2017linear}.  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)}~\d x = 1/2 $. In the work of~\cite{MANDAL2020}, 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---$\lc$---sensitivity. A strategy to mitigate this effect is to
introduce an energy density threshold, see for example~\cite{ttn2}.
While, the AT1 model is usually more suitable for brittle failure, it
may require additional techniques for implementation, as outlined
in~\cite{molnar2022thermodynamically}. 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:
{\begin{equation}\label{Wext}
\mathcal{W}_{\mathrm{ext}} =
\int_{\Omega } 
( \underline{b}  \cdot \underline{u} )  \, \d\Omega+ 
\int_{\partial \Omega_{n}} 
( \underline{t} \cdot \underline{n}_s\,  \underline{u}  ) 
\, \d S. 
\end{equation}}\unskip

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
{\begin{equation}\label{admissible_set_U}
\overline{\Omega}_{u}(\underline{v})=
\{\,\underline{u} : T \times \Omega \rightarrow H^1(\Omega, R^{\mathrm{dim}}), 
\underline{u} = \underline{v} 
\;\text{on}\;
\partial_{v}\Omega\}.
\end{equation}}\unskip
with $H^1(\Omega, R^{\mathrm{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 
{\begin{equation}\label{admissible_set_D}
\overline{\Omega}_{\bar{d}}(t)=
\{d : T \times \Omega \rightarrow [0, 1], ~
\bar{d} \leq d \leq 1, 
\; d = \bar{d} 
\; \text{on}\;
\partial_{\bar{d}}\Omega\}
\end{equation}}\unskip
with $\bar{d}$ denoting the solution $d$ before $t$.

\noindent 
The variational problem is written through the coupled system:
{\begin{subequations}\label{variational_both_equations} 
\begin{eqnarray}
\delta S(\underline{u}, d; \delta \underline{u}) 
&amp; = &amp;
0 \quad\mbox{in}\;
\forall \delta \underline{u} \in \overline{\Omega}_u, 
\label{eq:initOmega1}\Seqnsplit
\delta S(\underline{u}, d; \delta d) 
&amp; \geq &amp;
0\quad\mbox{in}\;
\forall \delta d \in \overline{\Omega}_d.
\label{eq:initOmega2}                  
\end{eqnarray}
\end{subequations}}\unskip
The variational form of the energies are:
{\begin{eqnarray}\label{dK}
\int_{t} 
\delta \mathcal{K}_u(\dot{\underline{u}}) \, \d t
= -  \int_{t} 
\int_{\Omega } 
( \rho_u \ddot{\underline{u}} \delta \underline{u}) \, \d\Omega  \,\d t , 
\delta \mathcal{K}_d(\dot{d}) \, \d t =
-  \int_{t} \int_{\Omega }
( \rho_d \ddot{d} \delta d ) \, \d\Omega  \, \d t 
\end{eqnarray}}\unskip
for the kinetic energies, 
{\begin{eqnarray}\label{dEel}
\int_{t} \delta \mathcal{E}_{el}(\underline{u}, d)  \, \d t
= \int_{t}   \int_{\Omega }  
\left(\rttensor{\sigma} : \delta 
\rttensor{\nabla u}
+ \frac{\partial \Psi}{\partial d} \delta d\right) 
\, \d\Omega  \, \d t , 
\end{eqnarray}}\unskip
for the strain energy,
{\begin{eqnarray}\label{dEfr}
\int_{t} \delta \mathcal{E}_{\mathrm{fr}}(d)  \, \d t 
= \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 \underline{\nabla} d} \right) \delta d
\right) \, \d\Omega \,  \d t 
\end{eqnarray}}\unskip
for the fracture energy and 
{\begin{equation}\label{weak_Wext}
\int_{t} \delta \mathcal{W}_{\mathrm{ext}} \, \d t
= \int_{t}  \int_{\Omega } 
( b \delta \underline{u} )  \, \d\Omega\, \d t+ 
\int_{t} \int_{\partial \Omega_{n}}  
( \underline{t} \cdot \underline{n}_s
 \, \delta \underline{u}  )  \, \d S  \,  \d t
\end{equation}}\unskip
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:
{\begin{subequations}\label{cauchy_pb_u} 
\begin{eqnarray}
&amp; \displaystyle
\int_{t} \int_{\Omega}   
( -\rho_u \ddot{\underline{u}} + 
\underline{\nabla} \cdot 
\rttensor{\sigma} + \underline{b})
\delta \underline{u} \, \d\Omega \, \d t 
=  0
\quad\mbox{in}\;
\forall \delta \underline{u} \in \overline{\Omega}_u,
\label{cauchy_pb_u_1}\Seqnsplit
&amp;\displaystyle
\int_{\partial \Omega} - 
\rttensor{\sigma} \cdot \underline{n}_s \delta \underline{u} 
\, \d S + 
\int_{\partial \Omega_{n}} 
\underline{t}\cdot \underline{n}_s \delta \underline{u} 
\, _d S =  0
\quad\mbox{on}\;
\partial \overline{\Omega}_u, 
\label{cauchy_pb_u_2}
\end{eqnarray}
\end{subequations}
\vspace*{2pt}
\begin{subequations}\label{hyperbolic_cauchy_pb_d}
\begin{eqnarray}
&amp;\displaystyle
\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 \, \d
t \geq 0
\quad\mbox{in}\;
\forall \delta d \in \overline{\Omega}_d, 
\label{hyperbolic_cauchy_pb_d_1}\Seqnsplit
&amp;\displaystyle
\int_{\partial \Omega_{d}} 
\frac{\partial \gamma}{\partial \underline{\nabla} d}
\delta d \cdot \underline{n}_d \, \d 
S\geq 0
\quad\mbox{on}\;
\partial \overline{\Omega}_d,
\label{hyperbolic_cauchy_pb_d_2}
\end{eqnarray}
\end{subequations}}\unskip
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.

\subsection{Derivation of damage PDEs}\label{sec2.2}

\subsubsection{Elliptic damage formulation: rate-independent
model~\cite{mielke2005evolution}}\label{sec2.2.1}
For this formulation, the damage kinetic energy is discarded. From the
crack density function introduced in \eq{fracture_density},  the final
strong form of the damage motion equation is given by:
{\begin{eqnarray}\label{elliptic_strong_form_d}
0 \geq  
- \frac{\partial \Psi}{\partial d}
- \frac{G_c }{\lc} 
(d - 2 \,\lc^2 \, \Delta d ).
\end{eqnarray}}\unskip
The damage PDE is reformulated with a damage yield  function $f(Y,
d)=Y- \sofrac{G_c }{\lc} 
(d - 2 \,\lc^2 
\, \Delta d )$ where
$Y$ represents an energetic  crack driving force defined by:
{\begin{equation}\label{crack_driving_force_0}
Y = - \frac{\partial \Psi}{ \partial d}
 = - g'(d) \frac{\partial \Psi}{\partial g}.
\end{equation}}\unskip

Then, a local form of the consistency  condition, similar to the
Karush--Kuhn--Tucker conditions,  can be  derived from the
irreversibility condition ($\dot{d}\geq0$) and the damage criterion
\mbox{($f(Y,d) \leq 0$):}
{\begin{eqnarray}\label{KKT_condition}
\Dot{d} \geq 0, 
\quad 
f(Y,d) \leq 0, 
\quad 
\Dot{d}f(Y,d) = 0.
\end{eqnarray}}\unskip
It is also known as the conditions of  consistency derived  from the
first order stability condition and the energy  balance proposed
in~\cite{pham2010approche}.

From the numerical point of view, the damage irreversibility (cf.\ 
Section~\ref{section:damage_irreversibility}) is handled within the
damage formulation considering an equality form of the evolution
equation  $f(Y,d) = 0$.

\subsubsection{Parabolic damage formulation: rate-dependent 
model~\cite{simo1987strain, cervera2003viscoelasticity}}
As for the elliptic formulation, the damage kinetic energy is discarded
but the system is made non conservative by considering a quadratic
dissipative potential
{\begin{equation}\label{dissipative_potential}
\Psi_{nc}= \frac{\mu_d}{2} \dot{d}^{2}
\end{equation}
that is incorporated in the principle of least action:
\begin{equation}\label{dActionplusdE_diss}
\int_{t} \delta L^{\mathrm{dyn}} \, \d t
\geq
\int_{t} \int_{\Omega } 
-\delta \Psi_{nc}(\dot{d})  \, \d\Omega  \, \d t
 = \int_{t} \int_{\Omega } 
 -\frac{\partial \Psi_{nc}}{\partial \dot{d}} 
 \delta d \, \d\Omega  \, \d t.
\end{equation}}\unskip
This leads to the local damage rate-dependent evolution equation:
{\begin{eqnarray}\label{parabolic_strong_form_d}
\mu_d \Dot{d}
\geq 
- \frac{\partial \Psi}{\partial d}
- \frac{G_c }{\lc} 
(d - 2 \,\lc^2 \, \Delta d )
 = f(Y,d).
\end{eqnarray}}\unskip

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.

\subsubsection{Hyperbolic damage formulation}\label{sec2.2.3}
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:
{\begin{eqnarray}
&amp; \displaystyle
\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 \, \d t
\geq
\int_{\Omega} \int_{t} - 
\left( \mu_d \dot{d} \,\delta d  \right)
\, \d t \, \d\Omega, 
\forall \delta d \in \overline{\Omega}_d
\label{Formulation_faible_d_Body_nonconservative}\Seqnsplit
&amp; \displaystyle
\rho_d \ddot{d} + \mu_d \Dot{d}
\geq 
- \frac{\partial \Psi}{\partial d}
- \frac{G_c }{\lc}
(d - 2 \,\lc^2 \, \Delta d )
 = f(Y,d).
\label{hyperbolic_strong_form_d}       
\end{eqnarray}}\unskip

\subsection{Tension-compression asymmetry models}\label{sec2.3}
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:
{\begin{equation}\label{energy_split}
\Psi(\rttensor{\varepsilon},d) = 
\tfrac{1}{2}
( g(d)\rttensor{\varepsilon}^+: \mathbb{C} : 
\rttensor{\varepsilon}^+ + 
\rttensor{\varepsilon}^-: \mathbb{C} : \rttensor{\varepsilon}^- )
= g(d) \Psi^+(\rttensor{\varepsilon}) +  
\Psi^-(\rttensor{\varepsilon}). 
\end{equation}}\unskip
The split is based on the definition of a positive and negative strain
$\rttensor{\varepsilon}^\pm$ and their corresponding stress
$\rttensor{\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: 
{\begin{equation}\label{stress_strain_pm_coax_isotropic}
\rttensor{\sigma}^{+}:\rttensor{\varepsilon}^{-} = 0, 
\quad \text{and} \quad
\rttensor{\sigma}^{-}:\rttensor{\varepsilon}^{+} = 0.
\end{equation}}\unskip
$\sofrac{\partial \Psi}{\partial d}$ can then be rewritten as 
{\begin{equation}
\frac{\partial \Psi}{\partial d} = g'(d) \Psi^+ = -2(1-d) \Psi^+.
\end{equation}}\unskip
The most used in the numerical fracture analysis are the  volumetric
and deviatoric split~\cite{LEMAITRE2000,  lancioni2009variational,
amor_regularized_2009}, the spectral  decomposition of the strain
tensor~\cite{PijaudierCabot2001ContinuousDM,  miehe_2010,
wu_cervera2018novel} or the masonry-like material
model~\cite{FREDDI20101154}. 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 $\mathbb{C}$ or the compliance tensor 
$\mathbb{S}$~\cite{desmorat2000dissymetrie, francois:hal-01004697}.
In~\cite{desmorat2000dissymetrie}, it is defined as a quadratic 
decomposition of the free energy using the square root stiffness tensor
$\mathbb{C}$: 
{\begin{equation}\label{square_root_decomposition}
\rttensor{\varepsilon} :\mathbb{C} : \rttensor{\varepsilon}
=  \tilde{\rttensor{\varepsilon}} : \tilde{\rttensor{\varepsilon}}
\end{equation}}\unskip
with $\tilde{\rttensor{\varepsilon}} = \mathbb{C}^{1/2}
\rttensor{\varepsilon}$. 

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$:
{\begin{equation}\label{transforemd_free_energy}
\Psi =  
\tfrac{1}{2} 
g(d) ~\rttensor{\tilde{\varepsilon}}^+ : 
\rttensor{\tilde{\varepsilon}}^+ + 
\rttensor{\tilde{\varepsilon}}^- : 
\rttensor{\tilde{\varepsilon}}^- .
\end{equation}}\unskip
As a consequence of this split, the strain energy definition is
modified following \eq{energy_split}.

A detailed 2D analytical analysis of this energy split model was 
carried out in~\cite{he_closed-form_2019} leading to its first 
application in phase field fracture model
in~\cite{nguyen_implementation_2020}. This approach opens up the
possibility to extend some decomposition  initially developed for
isotropic material to anisotropic  material~\cite{DIJK2020140,
ZIAEIRAD2023}, such as the  volumetric/deviatoric split model and the
no-tension split model~\cite{amor_regularized_2009, FREDDI20101154}.
Recently, another application of this decomposition in the context  of
heterogeneous material was run in a quasi-static regime 
in~\cite{thanh2023phase}.

Although anisotropic materials are not considered herein, the 
orthogonal split model formulated with a scalar damage variable from
\eq{transforemd_free_energy} is used  in this work as a generic
framework for implementing damage asymmetric effect, in 2D  case.

\section{Resolution strategy and
implementation}\label{section:numerical_strategy} %3

\subsection{Discretization}\label{sec3.1}
The weak form of the equations of the coupled problem is expressed as
$\forall t \in T$:
{\begin{eqnarray}
&amp;\displaystyle
\int_{\Omega} 
( \rho_u \ddot{\underline{u}}~ \underline{u}^* + 
\rttensor{\sigma}:
\underline{\nabla}~ \underline{u}^*) \, \d\Omega
= 
\int_{\Omega} \underline{b}~\underline{u}^*\, \d\Omega+
\int_{\partial \Omega}
(\rttensor{\sigma}~ \underline{n}_s)~\underline{u}^*
\, \d S
\quad 
\forall  \underline{u}^* 
\in \tilde{\Omega}_u(\underline{0}), 
\label{u_weak}\Seqnsplit
&amp;\displaystyle
\int_{\Omega} 
\left(\rho_d \ddot{d} ~ d^*+
\mu_d \dot{d} ~ d^*+ 
2~d~   d^*~\Psi^+ +
G_c \left(  \frac{d}{\lc} ~ d^* +
\lc \underline{\nabla}d \; \underline{\nabla} d^* \right)
\right)
  \, \d\Omega 
\nonumber\\
&amp;\displaystyle
=  \int_{\Omega} 2 \Psi^+ ~d^* ~\, \d\Omega 
\quad \forall d^* \in \tilde{\Omega}_d, 
\label{d_weak_hyp}              
\end{eqnarray}}\unskip
In these equations, $\underline{u}^*$ and $d^*$ are test functions for
the displacement and damage fields respectively that belong to the
admissible space defined by
{\begin{equation}\label{admissible_set_U_star}
\tilde{\Omega}_{u}(\underline{v})=
\{\,\underline{u} : \Omega \rightarrow H^1(\Omega, R^{\mathrm{dim}}), 
\underline{u} = \underline{v} 
\;\text{on}\;
\partial_{v}\Omega\}
\end{equation}}\unskip
and
{\begin{equation}\label{admissible_set_D_star}
\tilde{\Omega}_{\bar{d}}=
\{ d :  \Omega \rightarrow R\}.
\end{equation} }\unskip
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 
{\begin{equation}\label{discrete_disp_damage}
\underline{u} = \underline{\mathbf{N}} ~\mathbf{U}, 
\quad
d = \mathbf{N} ~\mathbf{D} 
\end{equation}}\unskip
with $\underline{\mathbf{N}}$ and $\mathbf{N}$ gathering vectors and 
scalar shape functions. Vectors like $\mathbf{U}$ or $\mathbf{D}$
collect the values of the nodal degrees of freedom for each field
variable. Using finite element shape functions as test functions, 
Equations~(\ref{u_weak})--(\ref{d_weak_hyp}) can therefore be collected
in a matrix vector format such as
{\begin{equation}
\mathbf{M}_u\ddot{\mathbf{U}}+
\mathbf{F}^i_u(\mathbf{U},\mathbf{D})=\mathbf{F}^e_u
\end{equation}}\unskip
for the displacement problem and 
{\begin{equation}\label{eq:damage_ef}
\mathbf{M}_d\ddot{\mathbf{D}}+\mathbf{C}_d\dot{\mathbf{D}}+
\mathbf{K}_d\mathbf{D}+\mathbf{F}^i_d(\mathbf{U})=\mathbf{0}
\end{equation}}\unskip
for the damage problem. Note that, the boundness and irreversibility
conditions are now additional constraints to be applied on the
numerical solution (see Sections~\ref{subsection:implementation} 
and~\ref{section:damage_irreversibility}). 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: 
{\begin{eqnarray}
\mathbf{M}_u 
&amp;=&amp;
\int_{\Omega_e} 
\rho_u \underline{\mathbf{N}}^T ~\underline{\mathbf{N}}
\, \d\Omega   
\label{mass_matrix_expression}\Seqnsplit
\mathbf{M}_d
&amp;=&amp; 
\int_{\Omega_e}   
\rho_d {\mathbf{N}}^T  ~\mathbf{N} 
\, \d\Omega 
\Seqnsplit
\mathbf{C}_d\dot{\mathbf{D}}
&amp;=&amp;
\int_{\Omega_e} 
\mu_d {\mathbf{N}}^T  \mathbf{N} \dot{\mathbf{D}} 
\, \d\Omega
= \mathbf{F}^{v}_d(\dot{\mathbf{D}})
\Seqnsplit
\mathbf{F}^i_u(\mathbf{U},\mathbf{D})
&amp;=&amp; 
\int_{\Omega_e} 
\rttensor{\mathbf{B}}^T :
\rttensor{\sigma} \, \d\Omega  
\label{f_int_expression}\Seqnsplit
\mathbf{K}_d
&amp;=&amp;
\int_{\Omega_e}  
\left(2\left( \Psi^{+} + \frac{G_c}{2l_c}\right) 
{\mathbf{N}}^T \mathbf{N} + 
G_c l_c ~ {\underline{\mathbf{B}}}^T \underline{\mathbf{B}} \right) 
\, \d\Omega  
\Seqnsplit
\mathbf{F}^i_d(\mathbf{U})
&amp;=&amp; 
- \int_{\Omega_e}  
2 \mathbf{N}~\Psi^{+}  
\, \d\Omega .
\end{eqnarray}}\unskip
where $\rttensor{\mathbf{B}}$  and $\underline{\mathbf{B}}$ contains
the vector and scalar shape function derivatives.

At each new time increment $t_{n+1}$, a staggered resolution strategy,
initially proposed  by~\cite{bourdin1998methode}, 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:

\beqimg\begin{eqimg*}
\begin{algorithm}[ht]
      \centering
      \caption{Central difference scheme for the displacement resolution}
      \begin{algorithmic}[1]
            \State \textbf{ \small Predictor}
            \State $\dot{U}_{n+1}^{\mathrm{pred}} \gets  \dot{U}_n + 0.5 \: \Delta t \: \Ddot{U}_n $
            \State $U_{n+1}  \gets U_n + \Delta t \dot{U}_n  + 0.5 \: \Delta t^2 \: \ddot{U}_n $
            \State  \textbf{ \small Solve motion equation at $t_{n+1}$}
            \State $ \mathbf{M}_u \ddot{U}_{n+1}  \gets 
            \mathbf{F}_u^e(U_{n+1})- \mathbf{F}_u^i(D_n, U_{n+1}) $
            \State  \textbf{ \small Corrector}
            \State $\dot{U}_{n+1} \gets \dot{U}_{n+1}^{\mathrm{pred}} +  0.5 \Delta t \: \ddot{U}_{n+1}$
      \end{algorithmic}
\label{cd_algorithm}      
\end{algorithm}
\eeqimg\end{eqimg*}

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 \tab{table:resolution_strategy}.

\begin{table} %t1
\caption{\label{table:resolution_strategy}Resolution strategy}
\begin{tabular}{ccccc}
\thead
\xmorerows{1}{Strategy} &amp; 
\multicolumn{2}{c}{PDE type} &amp; 
\multicolumn{2}{c}{Time stepping} \\
\cline{2-3} 
\cline{4-5} 
 &amp; $u$ &amp; $d$ &amp; $u$ &amp; $d$\\
\endthead
1 &amp; \morerows{2}{Hyperbolic} &amp; Elliptic &amp; \morerows{2}{CD (Algorithm~1)} &amp; --- \\
2 &amp;  &amp; Parabolic &amp;  &amp; Forward Euler (Algorithm~2) \\
3 &amp; &amp; Hyperbolic &amp; &amp; Verlet (CD) (Algorithm~3)
\botline
\end{tabular}
\end{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.    

\subsection{Resolution of the damage
problem}\label{subsection:implementation} %3.2
The various types of damage PDE for the phase field model presented in
the previous  section require  different resolution schemes (see
\tab{table:resolution_strategy}) 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. 

\subsubsection{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:
{\begin{equation}
\mathbf{K}_d\mathbf{D}+\mathbf{F}^i_d(\mathbf{U})=\mathbf{0}.
\end{equation}}\unskip
This mixed explicit/implicit approach is disadvantageous for large 
scale problems. In this case, \textit{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.

\subsubsection{Parabolic damage formulation}
In this case, the damage problem reduces to 
{\begin{equation}\label{eq:para_discrete}
\mathbf{C}_d\dot{\mathbf{D}}+\mathbf{K}_d\mathbf{D}+
\mathbf{F}^i_d(\mathbf{U})=\mathbf{0}.
\end{equation}}\unskip
The problem is integrated with the forward Euler scheme, 
(Algorithm~2). 
A lumped damping matrix $\mathbf{C}_d$ is used so
that the linear system of equations is diagonal. 

\beqimg\begin{eqimg*}
\begin{algorithm}[ht]
      \centering
      \caption{Forward Euler explicit scheme}
      \begin{algorithmic}[1]
            \State $\mathbf{C}_d \: \dot{\mathbf{D}}_{n+1}  \gets
            -\mathbf{K}_d\mathbf{D}_{n} - \mathbf{F}_{d}^i(\mathbf{U}_{n+1}) $
            \State $\mathbf{D}_{n+1}  \gets \mathbf{D}_n + \Delta t \: \dot{\mathbf{D}}_{n+1}$
      \end{algorithmic}
\label{forwardeuler}
\end{algorithm}
\eeqimg\end{eqimg*}

\subsubsection{Hyperbolic damage formulation}
When considering the hyperbolic formulation, the problem to solve is
driven by \eq{eq:damage_ef}. Here, the latter is solved with the Verlet
method (Algorithm~3), 
relying on the definition of  the
half-step velocity~\cite[Section~13.4.2]{schlick2010molecular}. The
viscous force $\mathbf{C}_d \dot{\mathbf{D}}$ is computed at  mid-step
\cite[Section~7.5.1]{geradin2014mechanical} without any
simplification. Here, the mass matrix $\mathbf{M}_d$ considered, is 
lumped by the row-sum method.

\beqimg\begin{eqimg*}
\begin{algorithm}[ht]
      \centering
      \caption{Verlet method with predictor--corrector form~\cite{kolman2020}}
      \begin{algorithmic}[1]
            \State \textbf{ \small Predictor}
            \State $\dot{\mathbf{D}}_{n+1}^{\mathrm{pred}} \gets  \dot{\mathbf{D}}_n + 0.5 \: \Delta t \: \Ddot{\mathbf{D}}_n $
            \State $\mathbf{D}_{n+1} \gets \mathbf{D}_n + \Delta t \dot{\mathbf{D}}_n  + 0.5 \: \Delta t^2 \: \ddot{\mathbf{D}}_n $
            \State  \textbf{ \small Solve motion equation at $t_{n+1}$}
            \State $ \mathbf{M}_d \ddot{\mathbf{D}}_{n+1}  \gets 
           -\mathbf{C}_d\dot{\mathbf{D}}_{n+1}^{\mathrm{pred}}-\mathbf{K}_d\mathbf{D}_{n+1}-\mathbf{F}^i_d(\mathbf{U}_{n+1})$
            \State  \textbf{ \small Corrector}
            \State $\dot{\mathbf{D}}_{n+1} \gets \dot{\mathbf{D}}_{n+1}^{\mathrm{pred}} +  0.5 \Delta t \: \ddot{\mathbf{D}}_{n+1} $
      \end{algorithmic}
      \label{Verlet_algorithm}
\end{algorithm}
\eeqimg\end{eqimg*}

\subsection{Damage irreversibility constraint
approximation}\label{section:damage_irreversibility} %3.3
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~\cite{delorenzis2020numerical}. 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~\cite{molnar2020open, nedjar256application}. 

\subsubsection{History variable}
An alternative to the previous strategies is to modify the energetic
crack  driving force $Y$ (\eq{crack_driving_force_0}). 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~\cite{miehe_2010}.  The classical
formulation is such as \eq{history_miehe},  where $\Psi_s$ is a
threshold (set to 0 in the following):
{\begin{equation}\label{history_miehe}
H(\underline{u}, t&gt;\tau)= 
\operatorname*{max}_{\tau \in[0,t]}
\Psi^+(\underline{\underline{\varepsilon}}(\underline{u}),\tau)
\quad \text{or}\;
\operatorname*{max}_{\tau \in[0,t]} 
(H(x, t), \Psi^+(\underline{\underline{\varepsilon}}
(\underline{u}),\tau)-\Psi_s).
\end{equation}}\unskip

This strategy has been tested in the context of hyperbolic damage
formulation~\cite{kamensky2018, parrinello_2017}. 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~\cite{BILGEN2019, HAI2020106821, rudshaug2024modeling}, based on
stress criterion (Rankine model, Tresca model and so on). This method
is very efficient, but the classical form proposed in~\cite{miehe_2010}
does not allow for the irreversibility condition to be met exactly, in
particular for hyperbolic damage formulation. This fact will be
illustrated in  Section~\ref{section:bar}.  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. 

\subsubsection{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{parabolic_strong_form_d}):
{\begin{equation}\label{parabolic_macc_irr}
\mu_d \dot{\mathbf{d}}= 
\langle  f(Y,d) \rangle _+
 = \left \langle - g'(d) \Psi^+
- \frac{G_c }{\lc} 
(d - 2 \,\lc^2 \, \Delta d) 
\right \rangle _+
\end{equation}}\unskip

Even if this modification cannot be included in the weak form of
\eq{parabolic_strong_form_d}, the latter is implemented in the
discretized system of  equations~(\ref{eq:para_discrete}) in a purely
algorithmic way. Therefore, Step 1 of 
Algorithm~2 is updated as
follows:
{\begin{equation}\label{parabolic_macc_disc}
\mathbf{C}_d \: \dot{\mathbf{D}}_{n+1} = 
\langle -\mathbf{K}_d\mathbf{D}_{n} -
\mathbf{F}_{d}^i(\mathbf{U}_{n+1}) \rangle_+
\end{equation}}\unskip
where the Macaulay brackets are applied component wise to vector of
nodal forces.

\subsubsection{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:
{\begin{equation}\label{hyperbolic_macc_disc}
\dot{\mathbf{D}}_{n+1}= 
\langle \dot{\mathbf{D}}_{n+1}^{\mathrm{pred}} +
0.5 \Delta t \: \ddot{\mathbf{D}}_{n+1} \rangle_+.
\end{equation}}\unskip
The acceleration is then modified accordingly such as
$\ddot{\mathbf{D}}_{n+1}=\sofrac{2}{\Delta t}
(\dot{\mathbf{D}}_{n+1}-\dot{\mathbf{D}}_{n+1}^{\mathrm{pred}})$. In
fact, this strategy is similar to what it is proposed for the parabolic
case.

\subsection{Numerical stability of the time
integrators}\label{sec:stability} %3.4
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:
{\begin{equation}\label{dt_u_damage}
\Delta t_{cr} =   \min(\Delta t_{cu}, \Delta t_{cd}), 
\quad 
\Delta t = c_{\mathrm{sec}} \Delta t_{cr} .
\end{equation}}\unskip
where $\Delta t_{cu}$ and $\Delta t_{cd}$ are the critical time steps
for the displacement and the damage problem respectively,
$c_{\mathrm{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 
{\begin{equation}\label{dtu_hyp}
\Delta t_{cu} = \frac{h_{\mathrm{min}}}{c}, 
\quad \text{with}\;
c = \sqrt{\frac{K+4\mu/3}{\rho}}
\end{equation}}\unskip
where $h_{\mathrm{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:
{\begin{equation}\label{dt_damage_parabolic}
\Delta t_{cd} \approx \frac{\mu_d h_{\mathrm{min}}^2}{G_c \lc}.
\end{equation}}\unskip

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_{\mathrm{sec}}$ can be used.

The critical time step associated to the hyperbolic damage formulation 
using the Verlet algorithm~\cite{serfHozHo2023method} is such as:
{\begin{equation}\label{dt_damage_hyperbolic}
\Delta t_{cd} = 
\frac{2(\sqrt{\xi^2 + 1} - \xi)}{\omega} = f(\rho_d, \mu_d)
\end{equation}}\unskip
with 
{\begin{equation}\label{xi_omegai_expression}
\omega = 
\frac{2}{h_{\mathrm{min}}} 
\sqrt{\frac{G_c \lc}{\rho_d}}, 
\quad
\xi = 
\frac{\mu_d}{2\omega \rho_d}.
\end{equation}}\unskip

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 
Appendices~\ref{appendix:stability_parabolic}
and~\ref{appendix:stability_hyperbolic}.

\subsection{Software framework}\label{sec3.5}
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~\cite{BarattaEtal2023}. Furthermore, the  implicit strategy is
implemented owing to  the PETSc library~\cite{balay2022petsc}. Lastly,
material behaviour knowledge  is handled by the code generator
MFront~\cite{helfer_MFront_2015}. 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. 

\section{Dynamic failure of a bar under tension}\label{section:bar} %4

{\vspace*{-2pt}}

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. 

{\vspace*{-4pt}}

\subsection{Numerical model}\label{sec4.1}

{\vspace*{-2pt}}

The geometry and the boundary conditions of the tensile bar are 
described in \fig{fig:model_bar}. 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  
\fig{fig:model_bar}. 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.

\begin{figure}
{\vspace*{-3pt}}
\includegraphics{fig01}
{\vspace*{-4pt}}
\caption{\label{fig:model_bar}Tensile bar  test---geometry and boundary
conditions.}
{\vspace*{-1pt}}
\end{figure}

The simulations are performed  with the mechanical properties of a
random steel-like material given in  \tab{table:material_parameter},
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_{\mathrm{L}}=L\sqrt{\rho_u/E}=0.205$~ms. The
value of this plateau that lasts $2t_{\mathrm{L}}$ is the applied stress
$\sigma_a$. Then the reaction force vanishes.   The bar is meshed with
$20\times400$ 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}
= \sofrac{3}{16} \sqrt{\sofrac{3EG_c}{\lc}}$.
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.

\begin{table} %t2
\caption{\label{table:material_parameter}Modelling parameters{\vspace*{-2pt}}}
\begin{tabular}{ccccccc} 
\thead
$h$ (m) &amp; $\lc$ (m)  &amp; $E$ (GPa) &amp; $G_c$ (N/m) &amp; $\nu$ &amp; $\rho_u$ (kg/m$^3$) &amp; $\sigma_R$ (MPa) \\
\endthead
0.0025 &amp; 0.0075 &amp; 190 &amp; 22130 &amp; 0 &amp; 8000 &amp; 243.2
\botline
\end{tabular}
{\vspace*{-1pt}}
\end{table}

The values of the time step for the simulations are presented  in
Section~\ref{sec:stability}, \eq{dt_u_damage}.

{\vspace*{-4pt}}

\subsection{Definition of some metrics for deviation
assessment}\label{subsection:error_metric_def} % 4.2

{\vspace*{-2pt}}

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_{\mathrm{sol}}$ with respect  to a reference solution
$x_{\mathrm{ref}}$, based on the mean squares difference  (MSD):
{\begin{equation}\label{msd_def}
\mathrm{MSD}_{\mathrm{rel}} = 
\frac{\mathrm{MSD}(x_{\mathrm{sol}})}{\mathrm{MSD}(0)}
\quad \mbox{with}\;
\mathrm{MSD}(x_{\mathrm{sol}}) = 
\frac{1}{N} 
\sum_{n=1}^{N}
(x_{\mathrm{sol}} - x_{\mathrm{ref}})^2,
\end{equation}}\unskip
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_{\mathrm{f}}$, 
identical for all the simulations, $d(x, 
t_{\mathrm{f}})$. In
the reminder of the section, the results obtained with the elliptic
formulation will be considered as the reference. 

\subsection{$\mu_d$-viscosity sensitivity analysis  of parabolic damage
PDE} %4.3
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
Equations~(\ref{parabolic_macc_disc}) and (\ref{history_miehe}).

\subsubsection{Parabolic formulation: damage irreversibility with
Macaulay brackets}\label{subsubsection:Macaulay_brackets} %4.3.1
The parabolic formulation (\eq{parabolic_macc_irr}) is first tested
with the  Macaulay brackets to ensure the damage irreversibility.
\fig{fig:Macaulay_reaction} illustrates the reaction  at the clamped
edge and \fig{fig:Macaulay_damage} 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_{\mathrm{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 \fig{fig:Macaulay_damage}a. At the fixed time
${t}_{\mathrm{f}}$, damage profiles show distinct states: increasing $\mu_d$
produces lower damage value at the left edge and confirms the delay
effect mentioned before.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:Macaulay_reaction}Damage irreversibility with
Macaulay brackets: sensitivity of the reaction $R_x(t)$ to
$\mu_d$.}
\end{figure}

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:Macaulay_damage}Damage irreversibility with
Macaulay brackets: sensitivity of the damage profile along the bar at a
fixed time $t_{\mathrm{f}}$ to $\mu_d$.}
\end{figure}

For high $\mu_d$ values $({\geq}15.0~\mathrm{Pa}{\cdot}\mathrm{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 \fig{fig:Macaulay_damage} 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$ (${\leq}2.0~\mathrm{Pa}{\cdot}\mathrm{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~\mathrm{Pa}{\cdot}\mathrm{s}$. This is
confirmed by \fig{fig:Macaulay_msd} that presents the metrics defined
in  Section~\ref{subsection:error_metric_def}. Indeed, the discrepancy
assessed on the  reaction evolution, in \fig{fig:Macaulay_msd}a, clearly shows
a  minimal value with the mechanical response of the elliptic  damage
results, which occurs at approximately 
$3.6~\mathrm{Pa}{\cdot}\mathrm{s}$. Regarding the damage
profile, \fig{fig:Macaulay_msd}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~\mathrm{Pa}{\cdot}\mathrm{s}$ for the reaction
and $25~\mathrm{Pa}{\cdot}\mathrm{s}$ for the damage profile.

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig:Macaulay_msd}MSD versus the damage viscosity
$\mu_d$. (a) Used the reaction step responses.
(b)~Used the damage profiles taken at the time  
$t_{\mathrm{f}}$ of
the fully crack stage of the elliptic damage result.}
\end{figure}

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 \fig{fig:Macaulay_msd}, 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~\mathrm{Pa}{\cdot}\mathrm{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{dt_damage_parabolic} 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}{\times}}4$) compared to
that of an explicit/implicit strategy. 

\subsubsection{Parabolic formulation: damage irreversibility with local
history variable}\label{subsubsection:irreversibility_localhistory}%4.3.2
The simulations are now performed with the local history  variable
$H$~\cite{miehe_2010} 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$ (\fig{fig:par_H_reaction}). The
time evolution shows an abrupt transition of the behaviour of the
numerical simulation results between 
$21.5~\mathrm{Pa}{\cdot}\mathrm{s}$ and 
$21.52~\mathrm{Pa}{\cdot}\mathrm{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~\mathrm{Pa}{\cdot}\mathrm{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
\fig{fig:irr0_miehe}. 
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~\mathrm{Pa}{\cdot}\mathrm{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~\mathrm{Pa}{\cdot}\mathrm{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{dt_u_damage}) is thus reduced down to 0.4, 
as shown in \fig{fig:Macaulay_reaction_stability}. 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.

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig:par_H_reaction}History-variable-based 
simulations---reaction forces at the clamped edge of the bar.}
\end{figure}

\begin{figure}
\includegraphics{fig06}
\caption{\label{fig:irr0_miehe}$\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.}
\vspace*{2pc}
\end{figure}

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig:Macaulay_reaction_stability}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_{\mathrm{sec}}\Delta t_{cd}$, $\Delta t_{cd}$ being lower than
$\Delta t_{cu}$ in these cases.}
\end{figure}

\subsection{$(\rho_d, \mu_d)$ sensitivity analysis of hyperbolic damage
PDE} %4.4
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{msd_def}), 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{hyperbolic_macc_disc}) and then using the local history
variable (\eq{history_miehe}).

\subsubsection{Hyperbolic formulation: damage irreversibility with
Macaulay brackets}\label{sec:hyperbolic_maccauley} %4.4.1
The analysis of the results starts by reporting on the MSD maps
presented in \fig{fig:msd_hyperbolic_mac}. 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\leq 5~\mathrm{Pa}{\cdot}\mathrm{s}$ and $\rho_d\leq
1.0 \times10^{-4}~\mathrm{kg}/\mathrm{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.

\begin{figure}
\includegraphics{fig08}
\caption{\label{fig:msd_hyperbolic_mac}Damage MSD map and reaction MSD
map for the hyperbolic formulation using Macaulay brackets.}
\end{figure}

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 \fig{fig:dt_d}
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 \fig{fig:msd_hyperbolic_mac}a. 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.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig:dt_d}Critical time step for the hyperbolic damage
PDE as a function of $\rho_d$ for different values of $\mu_d$.}
\end{figure}
 
In \fig{fig:msd_hyperbolic_mac}b, 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\times10^{-7}~\mathrm{kg}/\mathrm{m})$ and II 
$(\rho_d=4.0\times10^{-6}~\mathrm{kg}/\mathrm{m})$ in
\fig{fig:hyp_mac_reac_dam_1_2_3}a--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_{\mathrm{f}}=0.41$~ms, 
the value of $d$ at the left bar edge being below 0.6
for $\mu_d\geq10~\mathrm{Pa}{\cdot}\mathrm{s}$. For line III
$(\rho_d=1.2\times10^{-5}~\mathrm{kg}/\mathrm{m})$, due to the
increase of $\rho_d$ the results start deviating significantly from the
elliptic formulation predictions. However, line III in
\fig{fig:msd_hyperbolic_mac}b 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 \fig{fig:hyp_mac_reac_dam_1_2_3}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 ($\rho_d$, $\mu_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 \fig{fig:hyp_mac_reac_dam_1_2_3}f). 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.

\begin{figure}
\vspace*{2pt}
\includegraphics{fig10}
\vspace*{4pt}
\caption{\label{fig:hyp_mac_reac_dam_1_2_3}Reaction force evolution and
damage profile for ($\rho_d$, $\mu_d$) along line I, II and III of
\fig{fig:msd_hyperbolic_mac}b for the hyperbolic formulation with
Macaulay brackets.}
\vspace*{6pt}
\end{figure}

Moving to lines IV 
$(\rho_d=1.5\times10^{-4}~\mathrm{kg}/\mathrm{m})$, V $(\rho_d =
6.0\times10^{-4}~\mathrm{kg}/\mathrm{m})$ and VI $(\rho_d =
3.95\times10^{-3}~\mathrm{kg}/\mathrm{m})$, the MSD for the reaction
force and the damage profile significantly increases (see
\fig{fig:msd_hyperbolic_mac}). First, the reaction plateau now starts
at the same value as the plateau of the elastic simulation (see
\fig{fig:hyp_mac_reac_dam_4_5_6}a, c and 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_{\mathrm{L}}$. These
oscillations period increases when $\rho_d$ is increased from 
$1.5\times10^{-4}~\mathrm{kg}/\mathrm{m}$ (line IV) to 
$6.0\times10^{-4}~\mathrm{kg}/\mathrm{m}$ (line V) and then to
$3.95\times10^{-3}~\mathrm{kg}/\mathrm{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 \fig{fig:hyp_mac_reac_dam_4_5_6}d and~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
(\fig{fig:hyp_mac_reac_dam_4_5_6}a) follows the same trend as for
line V and VI (\fig{fig:hyp_mac_reac_dam_4_5_6}c and~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).
 
\begin{figure}
\includegraphics{fig11}
\caption{\label{fig:hyp_mac_reac_dam_4_5_6}Reaction force evolution and
damage profile for ($\rho_d$, $\mu_d$) along line IV, V and VI of
\fig{fig:msd_hyperbolic_mac}b for the hyperbolic formulation with
Macaulay brackets.}
\end{figure}

\subsubsection{Hyperbolic formulation: damage irreversibility with
local history variable}\label{sec:hyperbolic_history} % 4.4.2
The results of the simulations using the hyperbolic formulation and the
local history variable (\eq{history_miehe}) for prescribing damage
irreversibility are now analysed. Starting with the MSD maps in
\fig{fig:msd_hyperbolic_h}, 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. 

\begin{figure}
{\vspace*{2pt}}
\includegraphics{fig12}
{\vspace*{2pt}}
\caption{\label{fig:msd_hyperbolic_h}Damage MSD map and reaction MSD
map for the hyperbolic formulation using local history variable.}
{\vspace*{2pt}}
\end{figure}

In \fig{irreversibility_check}, 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\times10^{-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
\fig{fig:msd_hyperbolic_h}, this violation of the irreversibility
constraint is associated with higher deviation from the elliptic
solution. For the highest value of $\mu_d$ considered in
\fig{irreversibility_check}, irreversibility is recovered and the
deviation is much lower and similar to what was obtained in
\fig{fig:msd_hyperbolic_mac} when Macaulay brackets are used to impose
damage irreversibility. Another assessment of damage irreversibility
constraint is presented in~\fig{fig:msd_rate}. 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\geq
1\times10^{-3}~\mathrm{kg}/\mathrm{m}$ or 
$\mu_d\geq10~\mathrm{Pa}{\cdot}\mathrm{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.     

\begin{figure}
\includegraphics{fig13}
\caption{\label{irreversibility_check}$\langle \dot{d}\rangle_{-} $
maps in the $(x,t)$ plane of four $\mu_d$ values and
$\rho_d=1.2\times10^{-5}~\mathrm{kg}/\mathrm{m}$ (line  III). The
color scale is adapted to each maps $\dot{d}_{\mathrm{min}}
=-3\times10^5$,
$-2.5\times10^5\,-5\times10^3\, -4\times10^{-3}~\mathrm{s}^{-1}$ 
for~(a),  (b), (c), (d) respectively.}
\end{figure}

\begin{figure}
\includegraphics{fig14}
\caption{\label{fig:msd_rate}Normalized space--time cumulative negative
damage rate for the hyperbolic formulation using local history
variable.}
\end{figure}

{\vspace*{2pt}}

\subsection{Conclusion}\label{sec4.5}

{\vspace*{2pt}}

The detailed analysis of the results obtained for the simulation of the
dynamic failure of bar helps to draw some conclusions on several
aspects:
\begin{enumerate}[(1)]
\item 
Irreversibility constraint:
\begin{itemize}
\item 
for both parabolic and hyperbolic formulation, it is mandatory that
this condition is satisfied, otherwise numerical instabilities occur;
\item 
using Macaulay brackets  reveals to be an efficient strategy to ensure
damage irreversibility;
\item 
the constraint is not completely guaranteed with the local  history
variable proposed in~\cite{miehe_2010}. 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.
\end{itemize}
\item 
Resolution strategy:\Tg&lt;/p&gt;
\begin{itemize}
\item 
using the Macaulay brackets to prevent damage healing, the results show
that the critical time step defined in Section~\ref{sec:stability} is
appropriate; 
\item 
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;
\item 
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.
\end{itemize}
\item 
Physical interpretation:\Tg&lt;/p&gt;  
\begin{itemize}
\item 
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;
\item 
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.
\end{itemize}
\end{enumerate}\Tg&lt;/p&gt;  

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 \mbox{situation} will be analyzed to assess
the behaviour of parabolic and hyperbolic PDEs in a configuration that
is more  representative. 

\section{Kalthoff and Winkler test}\label{section:kalthoff} %5
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~\cite{kalthoff2000}. 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~\ref{section:crack_tracking}.

The geometry, boundary conditions and the loading are given by 
\fig{fig:model_kalthoff} and \eq{eq:kalthoff_loading} respectively. 
The material properties are similar to the previous test and defined 
in \tab{table:material_parameter} 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.
{\begin{equation}\label{eq:kalthoff_loading}
u_{\mathrm{applied}} =
\begin{cases}
0.5\dfrac{V_0}{t_0}t^2  &amp; 
t&lt;t_0=10^{-6}~\mathrm{s}, ~ V_0 = 16.5~\mathrm{m}/\mathrm{s}\\
0.5V_0^2 + V_0 (t-t_0) &amp; \text{otherwise}
\end{cases}
\end{equation}}\unskip

\begin{figure}
\includegraphics{fig15}
\caption{\label{fig:model_kalthoff}Kalthoff et~al.\ experiment: model
and geometry.}
\end{figure}

\subsection{Parabolic damage formulation results}\label{sec5.1}
We first start investigating the prediction of the parabolic damage
formulation.~$\mu_d$ is varied between 
$110.~\mathrm{Pa}{\cdot}\mathrm{s}$ to 
$35.0~\mathrm{Pa}{\cdot}\mathrm{s}$. The
Macaulay bracket are used to impose irreversibility to the damage
evolution.

In \fig{fig:crack_vel_pos_kalt_par}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\textdegree. 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.

\begin{figure}
\includegraphics{fig16}
\caption{\label{fig:crack_vel_pos_kalt_par}Crack tip kinematics
extracted at the same final time step. $c_{\mathrm{R}}$, 
the Rayleigh wave speed.}
\end{figure}

This can also be demonstrated through the crack tip velocity evolution 
in \fig{fig:crack_vel_pos_kalt_par}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  \fig{fig:crack_vel_pos_kalt_par}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.  

\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig17}
\vspace*{-4pt}
\caption{\label{fig:crack_vel_pos_kalt_hyp_s0}Crack path and crack tip
velocity for $\rho_d = 1.0\times10^{-5}$, $6.458\times10^{-5}$, 
$1.077\times10^{-4}~\mathrm{kg}/\mathrm{m}$.}
\end{figure}

\begin{figure}
\includegraphics{fig18}
\caption{\label{fig:crack_vel_pos_kalt_hyp_s1}Crack path and crack tip
velocity for $\rho_d = 5.0\times10^{-4}$,  $49.6\times10^{-4}$,
$45\times10^{-3}~\mathrm{kg}/\mathrm{m}$.}
\end{figure}

\subsection{Hyperbolic damage formulation results}\label{sec5.2}
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~\ref{fig:crack_vel_pos_kalt_hyp_s0}--\ref{fig:crack_vel_pos_kalt_hyp_s1}, 
the $\mu_d$-sensitivity of the crack path and velocity are plotted  for
six $\rho_d$ values: $1.0\times10^{-5}$, $6.458\times10^{-5}$,
$1.077\times10^{-4}$, $5.0\times10^{-4}$, $9.6\times10^{-4}$,
$5\times10^{-3}~\mathrm{kg}/\mathrm{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\times10^{-4}$ (\fig{fig:crack_vel_pos_kalt_hyp_s1}),
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\times10^{-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\times10^{-4}~\mathrm{kg}/\mathrm{m}$ (from
Figures~\ref{fig:crack_vel_pos_kalt_hyp_s0}e
and~f), a crack arrest is obtained
before the end of the simulation time. For $(\rho_d,\mu_d)$ equals to 
($1.077\times10^{-4}~\mathrm{kg}/\mathrm{m}$, 
$30.0~\mathrm{Pa}{\cdot}\mathrm{s}$)
(\fig{fig:crack_vel_pos_kalt_hyp_s0}f), and
$\rho_d=5.0\times10^{-4}~\mathrm{kg}/\mathrm{m}$
(\fig{fig:crack_vel_pos_kalt_hyp_s0}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:
\begin{itemize}
\item 
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; 
\item 
the transition between the acceleration phase and the constant velocity
plateau becomes smoother and smoother when $\rho_d$ is increased;
\item 
the crack initiation angle is decreased when $\rho_d$ is increased.
\end{itemize}
For a better assessment of this comments,
\fig{fig:crack_vel_pos_fixed_eta} presents the results for $\mu_d =
19~\mathrm{Pa}{\cdot}\mathrm{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.

\begin{figure}
{\vspace*{-4pt}}
\includegraphics{fig19}
{\vspace*{-4pt}}
\caption{\label{fig:crack_vel_pos_fixed_eta}Crack path and crack tip
velocity for $\mu_d = 19~\mathrm{Pa}{\cdot}\mathrm{s}$.}
{\vspace*{-2pt}}
\end{figure}

\section{Discussion and conclusion}\label{sec6}
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~\cite{miehe_2010}, 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.  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. 

\section*{Dedication}
The manuscript was written through contributions of all authors. All
authors have given approval to the final version of the manuscript.

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds
from any organization that could benefit from this article, and have
declared no affiliations other than their research organizations.

\back{}

\appendix{}

\def\appendixlabel{Appendix \Alph{section}}

\def\theequation{\thesection\arabic{equation}}

\section{Stability analysis of parabolic damage
PDE}\label{appendix:stability_parabolic} %Appendix A
Recalling, the matrix form of the parabolic formulation, 
{\begin{equation}\label{parabolique_matricielle}
  \mathbf{C}_d ~ \Dot{\mathbf{D}} + 
  \mathbf{K}_d  ~ \mathbf{D} =  \mathbf{F}.
\end{equation}}\unskip
for one degree of freedom (DoF), the matrix form is approximated by 
{\begin{equation}\label{c_k}
\left\{
\begin{array}{@{}ll}
 c = \mu_d h 
 \vspace*{2pt}\\
\displaystyle
 k =  \frac{G_c \lc }{h}
\end{array}
  \right.
\end{equation}}\unskip

The frequency is defined by the ratio between the coefficients \eq{c_k}
as 
{\begin{equation}\label{rapport_c_k}
\omega = \frac{k}{c} = 
\frac{G_c \lc}{\mu_d h^2}
\end{equation}}\unskip

Let's consider two DoFs. Using the Euler explicit scheme 
{\begin{equation}\label{euler_explicite}
\dot{\mathbf{D}}_{n+1} = 
\frac{1}{\Delta t} 
( \mathbf{D}_{n+1}-\mathbf{D}_{n}),
\end{equation}}\unskip

\eq{parabolique_matricielle} can be modified to introduce the
amplification matrix
{\begin{equation}\label{amplification_matrix_euler_explcite}
\mathbf{D}_{n+1}    
 =  \mathbf{A} ~ \mathbf{D}_n + 
 \Delta t ~ \mathbf{C}_d^{-1} \mathbf{F}
\quad \mbox{with}\;
\mathbf{A} = 
(\mathbf{I}
  - \Delta t ~ \mathbf{C}_d^{-1}  \mathbf{K}_d )
\end{equation}}\unskip

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{stability_euler_explicite}).
{\begin{equation}\label{stability_euler_explicite}
| \rho( \mathbf{A}) | &lt; 1
\end{equation}}\unskip

In the case of a 1D bar element, a condensed damping matrix {(row-sum
method)} and a stiffness matrix is  expressed as:
{\begin{equation}\label{stiff_damp_matrices_parabolic}
\mathbf{C}_d =  
\frac{h \mu_d}{2}
\begin{bmatrix}
1 &amp; 0 \\
0 &amp; 1
\end{bmatrix}
\quad \text{and} \quad
\mathbf{K}_d =  
\frac{G_c \lc}{h}
\begin{bmatrix}
1  &amp;  -1 \\
-1 &amp; 1
\end{bmatrix}
+
\frac{h}{6} 
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right)
\begin{bmatrix}
2 &amp; 1 \\
1 &amp; 2
\end{bmatrix}.
\end{equation}}\unskip
The amplification matrix is now
{\begin{equation}\label{amplification_matrix}
\mathbf{A}=
\begin{bmatrix}
\displaystyle
1-\frac{2}{\mu_d h}
\left( \frac{G_c\lc}{h} + \frac{h}{3}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) 
\Delta t
&amp;
\displaystyle
~~ -\frac{2}{\mu_d h}
\left( -\frac{G_c\lc}{h} + \frac{h}{6}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) \Delta t
\vspace*{4pt}\\
\displaystyle
-\frac{2}{\mu_d h}
\left( -\frac{G_c\lc}{h} + \frac{h}{6}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) \Delta t
&amp;
\displaystyle
~~ 1-\frac{2}{\mu_d h}
\left( \frac{G_c\lc}{h} + \frac{h}{3}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) \Delta t
\end{bmatrix}
\end{equation}}\unskip
and its eigenspectrum is given by the determinant
{\begin{eqnarray}
\det(\mathbf{A}- \lambda \mathbf{I}) 
&amp; =&amp; 
\left[1 - \lambda -\frac{2\Delta t }{\mu_d h}
\left( \frac{G_c\lc}{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\lc}{h} + \frac{h}{6}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) - \lambda \right]^2
\nonumber\\
&amp; = &amp;
\left[1 - \lambda -\frac{4\Delta t}{\mu_d h^2} 
\left(G_c\lc + \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].
\label{det_ampl_mat}
\end{eqnarray}}\unskip

The roots of the second order polynomial \eq{det_ampl_mat}, which 
results from $\det(\mathbf{A} - \lambda \mathbf{I})=0$, provides the 
maximum of the eigenvalues, 
{\begin{equation}
| \rho( \mathbf{A}) | =
\left| 1 -\frac{4\Delta t}{\mu_d h^2}
\left(G_c\lc +  \frac{h^2}{12} 
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) \right) \right|
\end{equation}}\unskip

The condition \eq{stability_euler_explicite} leads to the critical time
step  
{\begin{equation}\label{dt_cd_1_2}
\Delta t_{cd}^{1D} = 
\frac{\mu_d h^2}{ 2 G_c \lc + \frac{h^2}{6}
\left( 2 \Psi_+ + \frac{G_c}{l_c}\right) }  
\end{equation}}\unskip
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_+(\underbar{u})$ leads to a decrease
in the time  step assessment.

To avoid drastic decrease of $\Delta t_{cd}$, the following
approximation
{\begin{equation}\label{dt_cd_1_approx}
\Delta t_{cd}^{1D} = 
\frac{\mu_d h^2}{ 2 G_c \lc } (1 + \epsilon(h^2))
\end{equation}}\unskip
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 
{\begin{equation}\label{time_step_euler_explicit}
\Delta t_{cd}^{0D} = \frac{\mu_d h^2}{G_c \lc} 
\quad\text{and} \quad
\Delta t_{cd}^{1D} = \frac{\mu_d h_{\mathrm{min}}^2}{ 2 G_c \lc}
\end{equation}}\unskip

\section{Stability analysis of hyperbolic damage
PDE}\label{appendix:stability_hyperbolic} %Appendix B
Using the Verlet time integration scheme,  the discrete form of
\eq{d_weak_hyp} is such as:
{\begin{equation}\label{eq:hyp_damage_fem}
\mathbf{M}_d\ddot{\mathbf{D}}_{n}+
\mathbf{C}_d\dot{\mathbf{D}}_{n-1/2}+
\mathbf{K}_d\mathbf{D}_{n}+
\mathbf{F}^i_d(\mathbf{U})=\mathbf{0}  
\end{equation}}\unskip
with a central difference scheme for the acceleration and an upwind scheme
for the velocity: 
{\begin{equation}\label{eq:damage_dpp_and_dp_fem}
\ddot{\mathbf{D}}_{n} = 
\frac{\mathbf{D}_{n+1} - 2\mathbf{D}_{n} + \mathbf{D}_{n-1} }{2\Delta t^2}
\quad \text{and} \quad
\dot{\mathbf{D}}_{n-1/2} =
\frac{ \mathbf{D}_{n} - \mathbf{D}_{n-1} }{\Delta t}
\end{equation}}\unskip
Finally, the damage at $t_{n+1}$ results from the resolution of 
{\begin{equation}\label{d_n+1_verlet_0}
\mathbf{M}_d ~\mathbf{D}_{n+1} = 
\Delta t^2 \mathbf{F}_{d}(\mathbf{U}_{n})+
( 2 \mathbf{M}_d - \Delta t \mathbf{C}_d - 
\Delta t ^2 \mathbf{K}_d) \mathbf{D}_n+ 
( \Delta t \mathbf{C}_d - \mathbf{M}_d) 
\mathbf{D}_{n-1} 
\end{equation}}\unskip
what leads to the following expression: 
{\begin{equation}\label{d_n+1_verlet_1}
\mathbf{D}_{n+1} = 
\Delta t^2 \mathbf{M}_d^{-1} ~\mathbf{F}_{d}(\mathbf{U}_{n}) + 
( 2 \mathbf{I} - \Delta t \mathbf{M}_d^{-1} \mathbf{C}_d - 
\Delta t ^2 \mathbf{M}_d^{-1}\mathbf{K}_d) \mathbf{D}_{n}+ 
( \Delta t \mathbf{M}_d^{-1}\mathbf{C}_d - \mathbf{I} ) 
\mathbf{D}_{n-1}.
\end{equation}}\unskip

Considering, $\mathbf{X}_{n+1} = 
\begin{pmatrix} 
\mathbf{D}_{n+1} \\ 
\mathbf{D}_{n} 
\end{pmatrix}$, 
the numerical solution can be stated in the recurrent from 
{\begin{equation}\label{amplification_matrix_verlet_explicit}
\mathbf{X}_{n+1} = 
\mathbf{T} ~\mathbf{X}_{n} + \mathbf{P}_{n+1}
\end{equation}}\unskip
The matrix $\mathbf{T}$ and the vector $\mathbf{P}_{n+1}$ are the
amplification matrix and the external force
respectively
{\begin{equation}\label{X_T_P_verlet_explicit}
\mathbf{T} =  
\begin{bmatrix}
2 \mathbf{I} - \Delta t \mathbf{M}_d^{-1} \mathbf{C}_d - 
\Delta t ^2 \mathbf{M}_d^{-1}\mathbf{K}_d &amp; 
~~\Delta t \mathbf{M}_d^{-1}\mathbf{C}_d - \mathbf{I} \\
\mathbf{I} &amp; \mathbf{0}
\end{bmatrix},
\quad
\mathbf{P}_{n+1} = 
\begin{pmatrix} 
\Delta t^2 \mathbf{M}_d^{-1} ~\mathbf{F}_{d}(\mathbf{U}_{n}) \\ 
0
\end{pmatrix} 
\end{equation}}\unskip

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
\mathrm{i} \in \{1, \ldots, N\} 
\]}\unskip
where $\xi_i$, $\omega_i$ and $\underline{\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
eexpressions:
{\begin{equation}\label{modal_stiff_damp_matrices}
\mathbf{M}_d^{-1}\mathbf{C}_d =  
\begin{bmatrix}
2\xi_1\omega_1 &amp; 0 \\
0 &amp; 2\xi_2\omega_2
\end{bmatrix} 
\quad \text{and} \quad
\mathbf{M}_d^{-1}\mathbf{K}_d =  
\begin{bmatrix}
\omega_1^{2}  &amp;  0 \\
0 &amp; \omega_2^{2}
\end{bmatrix}.
\end{equation}}\unskip

Thus, the amplification matrix can be written as 
{\begin{equation}\label{amplification_modal_matrix_verlet_explicit}
\mathbf{T} =  
\begin{bmatrix}
2 - 2\Delta t \omega_1 \xi_1 - \Delta t ^2 \omega_1 ^2
&amp;  0 &amp;
\Delta t 2 \omega_1 \xi_1 - 1  &amp; 0
\vspace*{1pt}\\
0  &amp; 
2 - 2\Delta t \omega_2 \xi_2 - \Delta t ^2 \omega_2 ^2
&amp; 0 &amp; 
\Delta t 2 \omega_2 \xi_2 - 1
\vspace*{1pt}\\
1 &amp; 0 &amp; 0 &amp; 0\\
0 &amp; 1 &amp; 0 &amp; 0
\end{bmatrix}
\end{equation}}\unskip

As for the previous analysis on the parabolic PDE, the root
$(\lambda_1, \ldots\lambda_n)$ of the polynomial characteristics give
the eigenvalues of the amplification matrix $\mathbf{T}$:
{\begin{equation}\label{amplification_characteristic}
\det(\mathbf{T}-\lambda\mathbf{I})= 0 
\Longrightarrow 
\lambda^2 - 2 I_1 \lambda + I_2=0 
\quad \text{with}\; 
\left\{
\begin{array}{@{}ll}
\displaystyle
2~I_1 = \mathrm{tr}(\mathbf{T}) = 
\sum_{i=1}^{2} (2 - 2 \Delta t \xi_i \omega_i 
 - \Delta t^2 \omega_i ^2) 
 \vspace*{4pt}\\
\displaystyle
I_2 = \det(\mathbf{T}) = 
\prod_{i=1}^{2}(1-2\Delta t \xi_i \omega_i)
\end{array}
 \right.
\end{equation}}\unskip
given the stability condition 
{\begin{equation}\label{stability_verlet_explicit0}
|\max( \lambda_1,\ldots, \lambda_n) | &lt; 1 .
\end{equation}}\unskip

These eigenvalues are
{\begin{equation}\label{root_polynom_characteristic}
\left\{
\begin{array}{@{}ll}
\displaystyle
\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
\vspace*{4pt}\\
\displaystyle
\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{array}
\right.
\end{equation}}\unskip

Among the eigenvalues, $\lambda_1$ has the largest absolute value. The
stability conditions thus leads to
{\begin{equation}\label{stability_verlet_explicite}
| \lambda_1 | &lt; 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|&lt;1
\end{equation}}\unskip
with
{\begin{equation}\label{stability_final_condition_verlet_explicite0}
0&lt;\Delta t ^2 \omega_1^2 + 4 \Delta t \omega_1 \xi_1 - 4.
\end{equation}}\unskip
After some manipulations, one gets
{\begin{equation}\label{stability_final_condition_verlet_explicite1}
(\Delta t \omega_1)_{\pm} =  2( - \xi_1 \pm\sqrt{\xi_1^2 + 1}  ).
\end{equation}}\unskip
Eventually, the critical time step $\Delta t_{cd}$ to satisfy the
stability condition is an function of the damage damping and inertia:
{\begin{equation}\label{dt_verlet_explicite}
\Delta t_{cd} = 
\frac{2\left(\sqrt{\xi_1^2 + 1} - \xi_1 \right)}
{\omega_1} = f(\rho_d, \mu_d).
\end{equation}}\unskip

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, 
{\begin{eqnarray}
&amp; \displaystyle
\mathbf{M}_d^{-1}\mathbf{C}_d =  
\frac{\mu_d}{3\rho_d}
\begin{bmatrix}
2 &amp; 1 \\
1 &amp; 2
\end{bmatrix} 
\quad \rightarrow \quad
{\mathbf{M}_d^{-1}\mathbf{C}_d}^{DZ} =
\frac{\mu_d}{3\rho_d}
\begin{bmatrix}
3  &amp;  0 \\
0 &amp; 1
\end{bmatrix}
\label{damp_matrices_dz}\Seqnsplit
&amp; \displaystyle
\mathbf{M}_d^{-1}\mathbf{K}_d =  
\frac{2G_c \lc}{\rho_d h^2} 
\begin{bmatrix}
1  &amp;  -1 \\
-1 &amp; 1
\end{bmatrix} +
\frac{1}{3 \rho_d} 
\left(2 \Psi_+ + \frac{G_c}{l_c} \right)    
\begin{bmatrix}
2  &amp;  1 \\
1  &amp;  2
\end{bmatrix}
\nonumber\\
&amp; \displaystyle
\rightarrow \quad
{\mathbf{M}_d^{-1}\mathbf{K}_d}^{DZ} =  
\begin{bmatrix}
\displaystyle
\frac{4G_c \lc}{\rho_d h^2} +
\frac{1}{3 \rho_d} 
\left(2 \Psi_+ + \frac{G_c}{l_c} \right)  &amp;  0 \\
\displaystyle
0 &amp; 
\displaystyle
\frac{1}{\rho_d} 
\left(2 \Psi_+ + \frac{G_c}{l_c} \right) 
\end{bmatrix}
\label{stiff_matrices_dz}
\end{eqnarray}}\unskip

By identification with \eq{modal_stiff_damp_matrices}, the natural
frequencies and the damping ratio are defined by
{\begin{eqnarray}
&amp; \displaystyle
\omega_1^2 =  
\frac{4G_c \lc}{\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), 
\label{omegai}\Seqnsplit
&amp; \displaystyle
2\xi_1\omega_1 = 
3 \frac{\mu_d}{3\rho_d}, 
\quad
2\xi_2 \omega_2 =
\frac{\mu_d}{3\rho_d}
\label{xi}
\end{eqnarray}}\unskip

Considering the Taylor expansion of $\omega_1^2$ and $\omega_2^2$ with 
respect to $h$,
{\begin{equation}\label{TaylorExp_omegai}
\omega_1^2 = 
4 \frac{G_c \lc}{\rho_d h^2} + \epsilon(h^2),
\quad 
\omega_2^2 = 0, 
\end{equation}}\unskip
with finally
{\begin{equation}\label{xi_omegai_expression_appendix}
\omega_1 = 
\frac{2}{h}\sqrt{\frac{G_c \lc}{\rho_d}},
\quad
\xi_1 = \frac{\mu_d}{2\omega_1 \rho_d}
\end{equation}}\unskip

\section{Crack-tip tracking from a continuous damage
field}\label{section:crack_tracking} %Appendix C
A non-intrusive algorithm based on the clustering of the \say{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
\say{cracked}, pictured as  a colored empty circle in
\fig{fig:crack_advance_sketch}. At each time increment, the set of
\say{damaged} nodes is updated. Among these nodes, only ones downstream
the crack tip are tracked, while \say{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 \say{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 \say{damaged}
nodes clustering described above. A red box is plotted around the final
set of \say{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 \textit{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.

\begin{figure}
\vspace*{3pt}
\includegraphics{fig20}
\vspace*{5pt}
\caption{\label{fig:crack_advance_sketch}A tracking illustration of the
crack advance.}
\vspace*{8pt}
\end{figure}

\printbibliography
\refinput{crmeca20240702-reference.tex} 

\end{document}
               </pre></body></html>