Clément Jeremy Bleyer
email: jeremy.bleyer@epfl.ch@enpc.fr

Jean-François Molinari

Jérémy Bleyer

Clément Roux-Langlois

Jean-François Molinari

J Bleyer

J-F Molinari

Dynamic crack

Keywords: Dynamic fracture, Crack branching, Brittle materials, Phase-field model, Damage-gradient model

come

Introduction

Understanding the various mechanisms governing the dynamic propagation of a crack in a brittle medium is still a challenge. The difficulty lies in the strong interaction between stress concentrations at the crack tip, various non-linear phenomena occurring in the process zone, material heterogeneities at potentially different scales and dynamic stress redistribution due to waves emitted by the moving crack tip and reflections at the boundaries. For a single propagating crack, linear fracture mechanics relies on the balance between a crack driving force, the dynamic energy release rate, and a crack resisting force given by the fracture energy, which is assumed to be a material property. It also predicts that the Rayleigh wave speed c R is the limiting velocity of a mode I propagating crack [START_REF] Stroh | A theory of the fracture of metals[END_REF]. However, various experimental results have shown that this simple picture of dynamic brittle fracture is far from being complete. In particular, the existence of limiting velocities below c R , microscopic and macroscopic crack branching phenomena and a dependence of the fracture energy on the crack velocity have been observed experimentally.

Limiting speed and dynamic instabilities

In the absence of branching or before its occurrence, experiments have reported the following features of dynamic crack propagation [START_REF] Ravi-Chandar | An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching[END_REF][START_REF] Broberg | How fast can a crack go?[END_REF]:

cracks accelerate to a constant limiting velocity; this limiting velocity depends on the experimental setup for a given material; the stress intensity factor and surface roughness increase during the constant velocity phase.

Although the limiting velocity may depend on loading, boundary conditions or geometry, it has been measured to be in the range of 0.5-0.65c R for glass, 0.6-0.7c R for PMMA and 0.35-0.45c R for Homalite-100 [START_REF] Guerra Amaro | Dynamic fracture in brittle amorphous materials : Dissipation mechanisms and dynamically-induced microcracking in PMMA[END_REF].

Experiments have also reported that branching can be suppressed by constraining the crack path to a weaker interface or by drilling an array of holes ahead of the crack [START_REF] Washabaugh | A reconciliation of dynamic crack velocity and Rayleigh wave speed in isotropic brittle solids[END_REF]. In these cases, crack velocities can reach up to 90% of the Rayleigh wave speed.

In [START_REF] Gao | Surface roughening and branching instabilities in dynamic fracture[END_REF], Gao proposed a principle maximizing the fracture energy flux ( Γa = G(v a ) • v a with G the dynamic energy release rate and v a the apparent crack velocity) i.e. the energy being absorbed into the fracture process per unit time per unit length of the crack front. A wavy-crack model is introduced in which the crack tip follows a wavy path corresponding to a local (microscopic) crack tip speed which may be very high (close to c R ) whereas its apparent (macroscopic) speed v a is only a fraction of c R .

Experiments have evidenced different regimes for crack propagation controlled by dynamic instabilities [START_REF] Fineberg | Instability in dynamic fracture[END_REF]. At low velocities, the crack surface appears almost flat (mirror regime), velocity oscillations and an increase of surface roughening in the form of conic marks start to appear at higher velocities (mist regime), at even higher velocities microbranching and important surface roughening occur (hackle regime). A coalescence of microcracks in the crack process zone has also been proposed as a potential mechanism triggering this instability [START_REF] Ravi-Chandar | An experimental investigation into dynamic fracture: II. Microstructural aspects[END_REF]. Recent experiments over a wide range of crack velocities [START_REF] Guerra | Understanding fast macroscale fracture from microcrack post mortem patterns[END_REF][START_REF] Dalmas | Damage mechanisms in the dynamic fracture of nominally brittle polymers[END_REF] have confirmed that these instabilities are at the origin of the limiting speed: before attaining c R , such instabilities prevent the crack to further accelerate and the picture of a single propagating crack is not appropriate anymore.

Branching criterion

Yoffe showed that the circumferential stress at the crack tip reaches a maximum at a 60 • angle when the crack tip speed exceeds 0.6c s [START_REF] Yoffe | The moving Griffith crack. The London, Edinburgh[END_REF]. Unfortunately, this argument cannot explain the critical velocities and branching angles observed in the experiments.

Crack branching can be viewed as a consequence of an excess of available energy flowing to the crack tip which cannot be dissipated by a single crack propagation. Indeed, stress intensity factors (SIF) and energy release rates (ERR) increase prior to branching then drop since energy is now used to propagate two branches (explaining the smooth aspect of the crack surface just after branching). This suggests the existence of a critical value of SIF or ERR related to branching [START_REF] Seelig | On the interaction and branching of fast running cracksa numerical investigation[END_REF][START_REF] Ravi-Chandar | Dynamic fracture[END_REF]. Because of the dependence of SIF on crack velocity [START_REF] Freund | Dynamic fracture mechanics[END_REF], this critical value could also be related to a critical crack tip velocity. However, there is no universal relation between dynamic SIF and crack tip velocity [START_REF] Bouchbinder | The dynamics of rapid fracture: instabilities, nonlinearities and length scales[END_REF]. Besides, experiments indicate that branches propagate at roughly the same speed as the main crack before branching, therefore weakening the assumption that the crack tip velocity is the relevant criterion for crack branching [START_REF] Ravi-Chandar | An experimental investigation into dynamic fracture: II. Microstructural aspects[END_REF].

Velocity-toughening mechanism of the fracture energy

Experiments on PMMA for a strip geometry [START_REF] Sharon | Energy dissipation in dynamic fracture[END_REF] (similar to the one considered later in this paper) reported a strong increase of the apparent fracture energy with the crack velocity, ranging from approximately 1,000 J/m 2 for v = 0.2c R up to approximately 8,000 J/m 2 for 0.68c R . However, this strong increase is attributed to an important increase of relative surface area created by a crack advance through branching instabilities. The amount of fracture surface is shown to be linearly related to the energy flux into the crack tip so that the local fracture energy can still be considered as constant (equal to 1,100 J/m 2 ) for v ≥ v c = 0.35c R , where v c was identified as the critical velocity associated with the microbranching instability.

In contrast to these dynamic measurements of PMMA fracture energy, quasi-static experiments usually report values in the range of 300-400 J/m 2 . In [START_REF] Guerra | Understanding fast macroscale fracture from microcrack post mortem patterns[END_REF][START_REF] Dalmas | Damage mechanisms in the dynamic fracture of nominally brittle polymers[END_REF], a wedgesplitting experimental setup enabled to measure fracture energy at very low velocities. Between 0.11c R and 0.18c R , the fracture energy increases abruptly between 400 J/m 2 and 1,200 J/m 2 i.e. from a value consistent with quasi-static measurements to a value consistent with the results of [START_REF] Sharon | Energy dissipation in dynamic fracture[END_REF]. Interestingly, in the slow crack regime, a single propagating front is observed whereas, for v ≥ v a = 0.18c R nucleations of microdefects propagating at v a ahead of the main crack are observed. The collective dynamics of the main crack and microdefects yield an apparent crack velocity higher than v a .

These different results highlight the fact that the dynamic behaviour of crack propagation in PMMA is extremely complex, resulting from microdefects nucleation, microbranching instabilities, an increase of relative surface area and maybe a thickening of the damage zone. Experiments have also shown that microcracking and surface roughening are complex 3D phenomena [START_REF] Goldman | Origin of the microbranching instability in rapid cracks[END_REF] and that geometrical non-linearities at the crack tip may have to be taken into account to overcome LEFM deficiencies [START_REF] Livne | Breakdown of linear elastic fracture mechanics near the tip of a rapid crack[END_REF].

Finite element models of dynamic crack propagation

Unfortunately, an exhaustive modeling of all the previously mentioned phenomena would certainly require very fine 3D computations in a heterogeneous medium, which are still out of reach. Different classes of models have been proposed to simulate some of these dynamic crack features in brittle materials at a macro-scale.

Cohesive zone models [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF] represent the crack propagation process by considering a potential opening between two bulk elements. The cohesive law defines the constitutive relation between the surface traction and the relative opening displacement. By constraining the crack propagation along the element edges, this approach suffers from a mesh dependency at small scales but seems able to capture some features of crack branching patterns [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Falk | A critical evaluation of cohesive zone models of dynamic fracture[END_REF][START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] and the existence of a limiting velocity.

The extended finite element method (XFEM) enriches the finite element interpolation by adding either singularities or strong discontinuities for elements cut by the crack [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. The level-set is an efficient method to describe the crack path, and is now widely used [START_REF] Stolarska | Modelling crack growth by level sets in the extended finite element method[END_REF]. Although XFEM enables to have a discrete representation of the crack, one main drawback is that additional branching criteria and velocity toughening models are needed as input of the crack propagation algorithm to be able to obtain branched configurations [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Xu | Modeling of dynamic crack branching by enhanced extended finite element method[END_REF].

Non-local approaches regroup a wide range of models which represent the discrete crack by a continuous damage field and a regularization length to remove any mesh dependency. One can mention the non-local integral approach [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Jirasek | Nonlocal models for damage and fracture: comparison of approaches[END_REF][START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF], gradient-enhanced models [START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF], eigenerosion [START_REF] Pandolfi | An eigenerosion approach to brittle fracture[END_REF], peridynamics [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF], thick-level sets [START_REF] Moës | A level set based model for damage growth: the thick level set approach[END_REF] and phase-field approaches, which will be described later in more details. Generally, such non-local approaches do not require additional criteria to obtain branched patterns [START_REF] Karma | Unsteady crack motion and branching in a phase-field model of brittle fracture[END_REF][START_REF] Henry | Study of the branching instability using a phase field model of inplane crack propagation[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF].

In [START_REF] Henry | Fractographic aspects of crack branching instability using a phase-field model[END_REF], the Karma-Kessler-Levine phase-field model [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF] is used to study crack branching in 3D and its relation to fractographic patterns induced by the instability. Depending on the simulation parameters, 2D crack fronts or complex 3D patterns have been obtained. It is highlighted that crack velocities for 3D patterns are usually smaller than those observed for 2D patterns and are also associated with a higher total fracture surface.

Dynamic crack propagation using the phase-field approach has also been investigated in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Li | Variational approach to dynamic brittle fracture via gradient damage models[END_REF][START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]. In some cases, branching angles and crack velocities have been reported to be close to experimental observations whereas strange patterns have also been observed in other situations. In particular, the choice of the tensioncompression splitting has been shown to have an important influence on the observed pattern [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]. In [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF], crack branching is associated with a critical value of the crack surface velocity.

In [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF], mode I crack branching has been investigated using a peridynamic formulation of brittle fracture. Results have shown a dependence on loading conditions and sample geometry of branching angles and crack tip velocities. The authors proposed a stress-wave pile-up mechanism and damage spreading to the crack faces as the origin of crack branching.

It is to be noted that peridynamics rely on the same mechanical ingredients as phase-field modeling of brittle fracture, namely a linear elastic continuum model, a dissipation mechanism controlled by a material fracture energy and an intrinsic regularization length. It is therefore quite interesting that both approaches are able to reproduce characteristic aspects of crack branching without any additional criterion.

Objectives and organization of the manuscript

The purpose of the present work is to address the simulation of dynamic crack propagation using variational phase-field models. These models have emerged in the last decade as a promising tool for brittle fracture simulation as they do not require any a priori knowledge of the crack path or topology while being much less sensitive to the mesh discretization compared to cohesive zone models for instance. A large number of previous works on phase-field approaches for dynamic fracture mostly concentrated on numerical issues and qualitative aspects of the branched patterns. To our knowledge, a close inspection of the physical mechanisms at the origin of dynamic crack propagation has not been realized yet in the phase-field framework, especially regarding the conditions leading to crack branching. We aim here at demonstrating that such models can reproduce experimentally observed features of brittle dynamic fracture such as the existence of a limiting velocity and an increase of apparent fracture energy with the velocity.

In section 2, the general formulation of the phasefield model used in this work is briefly recalled. Constitutive modeling choices as well as numerical aspects are also discussed. In section 3, the numerical simulation of dynamic crack propagation in a pre-strained plate is investigated and the occurrence of crack branching and the existence of a velocity-toughening mechanism are more particularly examined. Finally, section 4 is devoted to the simulation of crack propagation in a heterogeneous medium. Situations in which the crack is constrained to propagate along a weak plane as well as propagation in presence of distant heterogeneities are investigated.

Variational phase-field models of brittle fracture

The phase-field models proposed in the literature as regularized models of brittle fracture can be split in two categories. The first category is based on a Ginzburg-Landau phase transformation evolution equation [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Henry | Study of the branching instability using a phase field model of inplane crack propagation[END_REF][START_REF] Henry | Fractographic aspects of crack branching instability using a phase-field model[END_REF]. Although this first approach enabled to produce interesting results regarding crack branching, mixedmode instabilities and surface roughening, some issues concerning the evolution of the phase field away from the crack tip, even if the sample has been brought to equilibrium, have been raised in [START_REF] Bourdin | A timediscrete model for dynamic fracture based on crack regularization[END_REF]. Besides, this formulation seems more distant to Griffith's theory of fracture although some physical parameters can be related to a fracture energy and a regularization length.

The second category of phase-field models has been developed as a regularization of free discontinuity problems [START_REF] Ambrosio | Approximation of functional depending on jumps by elliptic functional via tconvergence[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] arising in the variational approach to fracture proposed in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. In this case, a continuous function d varying between 0 (sound material) and 1 (fully cracked material) is introduced to obtain a smooth formulation of the Griffith fracture functional by a crack density function depending on d and ∇d. This functional is parametrized by the Griffith fracture energy G c and an internal length scale l 0 , which is related to the distance over which the phase field varies from 0 to 1 in a localized damaged band. The majority of the phase-field models adopted in the literature are based on this formulation, either in a rate-independent form [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Bourdin | A timediscrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF] or in a rate-dependent form including viscous dissipation of crack propagation as a numerical regularization of the time evolution problem [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF].

Continuum descriptions of cohesive fractures have also been proposed using the phase-field approach [START_REF] May | A numerical assessment of phase-field models for brittle and cohesive fracture: Γ -convergence and stress oscillations[END_REF][START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], continuum damage-gradient [START_REF] Lorentz | Modelling large crack propagation: from gradient damage to cohesive zone models[END_REF] or the thick-level set (TLS) method [START_REF] Moës | A level set based model for damage growth: the thick level set approach[END_REF], see also [START_REF] Cazes | Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture[END_REF] for a comparison between the TLS and the phase-field approach. Recent works also proposed to extend the phase-field approach to the case of ductile fracture in an elastoplastic-damageable formulation [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF].

Phase-field models based on a regularization of the variational approach to fracture share strong links with continuum damage-gradient models [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. In fact, the phase-field model based on [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] belongs to the general class of gradient damage models for a specific choice of damage constitutive laws. These damage-gradient models have been recently extended to the dynamic propagation in [START_REF] Li | Variational approach to dynamic brittle fracture via gradient damage models[END_REF][START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]. Due to the variational formulation of these models and their generality, we will adopt this point of view in the remainder of this work.

Finally, let us mention that, for crack-like localized solutions, quasi-static damage-gradient models converge to the Griffith theory of brittle fracture when the internal length scale l 0 is much smaller than the characteristic length L of the structure [START_REF] Sicsic | From gradient damage laws to Griffiths theory of crack propagation[END_REF]. Although no formal proof has been published yet concerning the dynamic case, it seems reasonable to expect that such a result will still hold in this case.

Formulation of the variational approach to dynamic fracture

We briefly recall here the general formulation of the variational phase-field or gradient damage model of dynamic fracture. We invite the reader to refer to [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF] and to references mentioned therein for more details on the construction of such models.

In the context of isotropic damage models, a continuous scalar field d, taking values in the [0; 1] interval, is introduced as a smooth representation of a transition between intact material (d = 0) and a fully cracked material (d = 1) (see Figure 1). Denoting by u a kinematically admissible displacement field and by u its associated velocity, the variational gradient damage model relies on the definition of the following energies:

the elastic strain energy: where Ω represents the continuum (Fig. 1), ε = 1 2 (∇u + ∇u T ) is the linearized strain tensor associated with u, ψ being the elastic potential depending on d; the kinetic energy:

E el (u, d) = Ω ψ(ε, d) dΩ (1)
E kin ( u) = Ω 1 2 ρ u • u dΩ (2)
where ρ is the material density; the non-local damage or fracture energy:

E f rac (d) = Ω w(d) + W 0 l 2 0 ∇d • ∇d dΩ (3)
where w(d) is the local damage energy density, W 0 is homogeneous to an energy density and l 0 is homogeneous to a distance.

Finally, an action-integral over a time interval I = [t 1 , t 2 ] is introduced as follows:

A = I (E el (u, d) + E f rac (d) -E kin ( u) -W ext (u)) dt (4)
where W ext (u) represents the work done by the external forces on the displacement field u.

The variational formulation of the gradient damage model, giving the space-time evolution of the mechanical fields (u, d), is then obtained by assuming that:

the damage variable evolution is irreversible: ḋ ≥ 0 for all t; the variation action-integral is positive with respect to arbitrary variations of admissible displacements and damage evolution; energy is dissipated only via the fracture energy

E f rac (d).
The derivation of these principles, especially the second one, shows that the elastodynamic equilibrium equation is obtained for a damage-dependent stress state σ = ∂ ε ψ(ε, d):

div σ + f = ρü (5)
with appropriate boundary conditions and where f is a given body force density.

Besides, the variation of (4) in terms of admissible damage evolution results in the total energy satisfying a minimum principle constrained by the irreversibility condition of damage evolution:

E el (u, d) + E f rac (d) ≤ E el (u, d) + E f rac ( d) ∀ d ∈ D(d) (6)
where D(d) is the admissible space for damage variations from a given damage state 0 ≤ d ≤ 1 defined by:

D(d) = d s.t. 0 ≤ d ≤ d ≤ 1 (7)
In practice, the wave equation ( 5) and the damage minimum principle ( 6) are solved numerically at every time step.

Constitutive modeling choices

The previous formulation is very general and different damage models can be obtained depending on the choice of the damageable elastic energy density ψ(ε, d) and the local damage density w(d).

For the elastic strain energy density, we chose the classical form [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]:

ψ(ε, d) = (1 -d) 2 ψ + (ε) + ψ -(ε) (8) ψ + (ε) = κ 2 tr ε 2 + + µε d : ε d (9) ψ -(ε) = κ 2 tr ε 2 - ( 10
)
where κ and µ are the compressibility and shear elastic moduli, ε d is the deviatoric strain tensor and ± = ( ± | |)/2 denotes the positive (resp. negative) part of . This model enables to distinguish between tension and compression, the damage degradation impacting the deviatoric strain as well as the positive part of the volumetric strain, while the negative part of the volumetric strain is not impacted. This is not the only choice for tension-compression splitting, see [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF] for an extensive discussion on this subject. For the problem considered in this paper, the impact of the splitting choice is limited as compression stress waves reflected at the boundaries are generated only by the crack tip advance and are of small amplitude.

As regards the choice of the damage energy density, the most widely used model [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Cazes | Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture[END_REF] is given by 1 :

w(d) + W 0 l 2 0 ∇d • ∇d = G c 2l 0 d 2 + l 2 0 ∇d • ∇d (11)
where G c is the Griffith fracture energy. Although appealing from a numerical point of view due to its purely quadratic nature, this model predicts no elastic domain in which the material behaves elastically with zero damage [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. Therefore, when comparing the macroscopic total dissipated energy with respect to a crack length evolution, a non-negligible part of this dissipated energy actually comes from small damaging of the bulk domain due to the applied stress state. This pre-damaging of the surrounding bulk domain will also tend to reduce the values of the different wave speeds compared to a non-damaged domain. Besides, the localized 1D solution for d of a bar in tension corresponds to an exponentially decreasing function with infinite support (Fig. 2).

For these different reasons, we chose to use the following damage energy density proposed in [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]:

w(d) + W 0 l 2 0 ∇d • ∇d = 3G c 8l 0 d + l 2 0 ∇d • ∇d (12)
1 with sometimes a non-essential rescaling of the internal length l 0

The absence of a square exponent in the damage variable leads to a pure elastic phase associated with a critical stress σ c = 3G c E/(8l 0 ) for a 1D homogeneous traction test. In this case, the 1D profile corresponds to a portion of parabola with an exactly zero damage state for |x| ≥ 2l 0 (Figure 2).

Let us also mention that the normalization of both models is chosen so that, when considering the localized damage solution for a 1D bar in traction, the surface fracture energy corresponds exactly to G c .

Numerical aspects

The numerical code used in this work is inspired from the open-source implementation of the dynamic damagegradient model [START_REF] Li | FEniCS (Dynamic) Gradient Damage[END_REF] based on the FEniCS project for automated resolution of PDE's [START_REF] Logg | Automated solution of differential equations by the finite element method: The FEniCS book[END_REF] and the PETSc library [START_REF] Balay | PETSc Web page[END_REF].

Classical linear finite element interpolations on triangles are used for the displacement, velocity, acceleration and damage fields. As regards temporal discretization, an explicit Newmark scheme has been chosen for the update of accelerations, velocities and displacements. The damage problem is formulated as a bound-constrained quadratic optimization problem to ensure damage irreversibility. More details on the numerical implementation can be found in [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF].

The mesh size has usually been taken approximately 4-5 times smaller than l 0 , thus reducing the numerical overestimation of the fracture energy [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]. Time steps have been chosen sufficiently small to satisfy the conditional stability of the explicit scheme. Mesh size and time steps have been varied to ensure that converged results have been obtained. Besides, it has also been observed that energy conservation is ensured at a satisfying accuracy for small time steps as mentioned in [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF].

Although the following results will be presented for a constant regularization length of l 0 = 0.1 mm, let us highlight that this value has been also varied. Diminishing this value resulted in a decrease of the damaged band widths as expected but crack patterns as well as crack velocities remained essentially the same. Indeed, as long as l 0 L and the solution is localized in a damaged band, reducing l 0 will only produce a self-similar thinner damage profile in the direction normal to crack propagation whereas the crack will still propagate according to a similar energy release rate. [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] for different imposed displacements (taken from [START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF]). For low loading levels, single crack propagation is observed. For higher values, microbranches appear during propagation. For even higher loadings, macroscopic branches are formed.

3 Dynamic crack branching in homogeneous medium

Position of the problem

We investigate the problem of dynamic crack propagation in a pre-strained PMMA rectangular plate geometry. This problem has been previously investigated experimentally in [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] (it is also quite close to the experiments of [START_REF] Sharon | Energy dissipation in dynamic fracture[END_REF]) and from a numerical point of view using cohesive elements [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] and a non-local integral damage model [START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF]. Contrary to more traditional numerical benchmark problems in dynamic branching such as those investigated in [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF] where a constant stress echelon loading is applied, in this particular geometry, no stress waves are produced by the loading but only by crack propagation.

Indeed, the experimental setup of [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] consists, first, in moving the upper and lower surfaces of the plate in a quasi-static manner until reaching a desired displacement. While maintaining the upper and lower boundaries at this fixed value, a small sharp crack is created at the middle of the sample using a razor edge. The crack then accelerates until reaching a steady-state propagation regime. Another advantage of this problem is that the initially stored energy is well defined and can be related to the fracture energy dissipated by crack advance. This pre-strained strip configuration has been used in many experiments to measure the dependence of fracture energy on crack velocity. For all these reasons, we advocate this problem as a benchmark for assessing numerical models of dynamic brittle fracture.

The geometry and boundary conditions are represented in Figure 3 and are identical to those used in [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF][START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF]: the plate is 32 mm wide and 16 mm high with a 4 mm pre-notch and modeled as a 2D plane-stress medium. It is pre-strained by applying uniform displacements ±∆U in the vertical direction on the bottom and top surfaces. First, a quasi-static computation is done to reach the pre-strained state, then an explicit dynamic computation is carried out to simulate the crack propagation.

PMMA is assumed to be isotropic linear elastic with E = 3.09 GPa, ν = 0.35, ρ = 1180 kg/m 3 , the Rayleigh wave speed is c R = 906 m/s. As regards the damage parameters, a value of G c = 300 J/m 2 is retained, corresponding to quasi-static measurements and the regularization length is l 0 = 0.1 mm.

A comment about the infinite strip configuration

Let us mention that the previous problem is close to the infinite strip configuration, see [START_REF] Freund | Dynamic fracture mechanics[END_REF] for instance. For this translation-invariant problem, a steady state propagation occurs with a fracture energy Γ equal to the initially stored energy per unit length, given by W = 2E(∆U ) 2 /h where h is the plate height. This configuration seems thus appropriate to measure the possibly velocity-dependent fracture energy.

An approximate equation of motion for the infinite strip problem has been derived in [START_REF] Marder | New dynamical equation for cracks[END_REF] using a perturbation approach:

v = 1 -Γ (v)/W h/2 c 2 d (1 -v 2 /c 2 R ) 2 (13)
As W can be arbitrarily chosen, when W > max v Γ (v) the crack is supposed to accelerate up to the Rayleigh wave speed but with an increasing effective mass so that acceleration becomes always more difficult when approaching c R [START_REF] Bouchbinder | Branching instabilities in rapid fracture: Dynamics and geometry[END_REF]. The identification of Γ (v) to W can then only be valid for a true steady state. For high velocities, an even small value of acceleration can correspond to an important difference between Γ (v) and W , leading to an overestimation of Γ (v).

Numerical computation of the macroscopic fracture energy

The crack tip position has been recorded by considering the top-rightmost node of the finite element mesh for which the damage field is greater than 0.9, its horizontal position is denoted by a(t) with the origin taken at the prenotch tip (a(0) = 0). Hence, in the case of macroscopic branching, only the upper-branch tip is followed. Unless stated otherwise, the crack tip velocity v corresponds to the horizontal component of the instantaneous velocity vector of the crack tip ȧ. This velocity is computed by fitting an affine function for a few simulation outputs around a(t). The damage dissipation rate Γ , which can be interpreted as a macroscopic fracture energy, has been computed as the derivative of dissipated energy per unit crack advance

Γ = - dE el da - dE kin da = dE f rac da (14)
since the boundaries of the domain are fixed during crack propagation. Γ has then been estimated by fitting an affine function for a few simulation outputs around E f rac (a). Let us also mention that the apparent dynamic energy release rates can also be computed using the G -θ method (virtual domain extension) generalized to such damage-gradient models [START_REF] Li | Variational approach to dynamic brittle fracture via gradient damage models[END_REF]. Despite slight differences, the values computed using such a method are consistent with those derived from (14).

Dynamic crack branching results

Simulations have been undertaken for different values of the imposed displacement ∆U . Evolutions of crack velocities are reported in Figure 4 in which it can be observed that cracks accelerate up to a limiting velocity which is below the Rayleigh wave speed c R , around 0.68c R . The acceleration is stronger for higher loadings and it decreases when the crack speed approaches the final velocity, which is consistent with experimental observations and the infinite-strip equation of motion [START_REF] Bourdin | A timediscrete model for dynamic fracture based on crack regularization[END_REF].

For small values of the loading (∆U < 0.04 mm), no macroscopic branching has been observed during the simulation. For higher loadings, branching has been observed and it occurs sooner the higher the loading. phase-field distributions corresponding to different loading levels have been represented in Figure 5.

Let us observe that these simulations yield crack velocities that are higher than those observed experimentally for the same loading. As a consequence, macroscopic branching also occurs at lower loading levels than in the experiments. For instance, experiments with ∆U = 0.06 mm reported a steady state velocity of 0.37c R with a single crack propagation.

This discrepancy between simulations and experiments has always been observed in previous works using cohesive zone models [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] or non-local integral damage models [START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF]. Both works mentioned the lack of ratedependency of the constitutive model as the origin of this difference. Inclusion of a viscous dissipation energy, as initially introduced in [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF] and proposed in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] as a purely artificial regularization to improve numerical stability, would enable to introduce a rate-dependent effect reducing the value of the computed velocities. However, this aspect would be the purpose of another work and we will not attempt here to reproduce quantitatively the measured velocities but we will concentrate on other aspects of dynamic crack propagation such as the transition from single crack propagation to branching. Besides, we expect that the rate-dependency will have less influence for fast cracks reaching a steadystate for which the material response will be quasiinstantaneous contrary to slow or accelerating cracks. As the branching mechanism occurs at a high limiting velocity, rate-dependent effects may be less important in this context.

A striking result is that the damaged band widens as the crack propagates, which has also been previously reported in other phase-field simulations [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF] but also with other models such as a process-region cell model [START_REF] Johnson | Process region changes for rapidly propagating cracks[END_REF] or peridynamics [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF]. This "widening" is mild for low loading levels but it is clearly visible for higher levels, especially before a branching event. For this reason, it is not easy to detect the exact origin of branching as it does not seem to be an abrupt phenomenon but rather a progressive one, exhibiting a continuous transition from an increasingly wider band to almost two crack tips propagating horizontally. At some point, when the two tips are sufficiently established, they start to screen and repel each other, leading to a true bifurcation and propagation of two isolated branches with a smaller damage band width. If the branching angle is measured as the angle made between the straight parts of the two isolated branches, a value close to 30 • is found, which is consistent with experiments and other simulations [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF]. However, we have also observed that different values of the branching angle can be obtained by varying the domain geometry (it is smaller for plates with a smaller height h).

The evolution of Γ (t) as a function of the crack tip position a(t) enables to better understand how the damaged band width increases with crack propagation. Indeed, denoting by l(t) the total damaged surface (per unit length in the transverse direction), let us recall that Fig. 6 The normalized damage dissipation rate Γ/G c increases during crack propagation. For the low loading level (∆U = 0.035 mm), the evolution is regular with G c ≤ Γ ≤ 1.5G c and no branching is observed. This increase leads to an effective thickening of the crack width. For a higher loading (∆U = 0.045 mm), branching is observed slightly after Γ ≥ 2G c , corresponding to two crack tips. After branching, the dissipation associated with a single crack tip (Γ/2) is close to its initial value, slightly above G c .

E f rac (t) = G c Ω 3 8l 0 d + l 2 0 ∇d • ∇d dΩ = G c • l(t) (15)
since G c is uniform in the whole sample. Therefore,

Γ (t) = G c dl da
, so that the normalized damage dissipation rate Γ/G c is equal to the increase in total fracture surface l(t) per apparent crack surface increase (if one interprets a(t) as an apparent crack surface). This quantity has been represented in Figure 6 for two different loading values: ∆U = 0.035 mm for which no macroscopic branching takes place and ∆U = 0.045 mm for which macroscopic branching is observed. After a first phase of initiation, the crack advances with values of Γ slightly above G c . The damage dissipation rate then increases during crack advance with an almost constant rate depending on the loading level. As mentioned before, this increase is directly related to the increase of total over apparent fracture surface due to the damaged band widening. It is to be noted that, for the duration of the simulation with ∆U = 0.035 mm, Γ increased from G c to 1.5G c and no branching has been observed. For the higher loading (∆U = 0.045 mm), branching has been observed slightly after Γ exceeded 2G c , a value corresponding to two crack tips, although detection of the exact instant of branching is difficult (this is symbolized by the gray zone in Fig. 6, the corresponding values of Γ have not been represented during the transition from single to branched crack). After branching, the total damage dissipation rate continued to increase and when looking only at the contribution of a single branch (given by Γ/2 due to the symmetry of the branched pattern, see Fig. 5(d)), values close to the initial ones (slightly above G c ) are obtained. This last observation is related to the fact that the damage band width of the branched parts are smaller than the main crack just before branching. Finally, the increase of Γ/2 after branching suggests that a new branching event can occur if this quantity again reaches a critical value close to 2G c . We are not aware of previous works reporting such numerical results pointing towards an energetic criterion for the occurrence of branching. In particular, our results indicate that branching can be triggered without stress waves arriving at the crack tip as already indicated in [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF]. However, incoming waves at the crack tip can lead to a sudden increase of the dynamic energy release rate, which will exceed its critical value and form branches. This seems to be the case for the more traditional stress echelon problem [START_REF] Bobaru | Why do cracks branch? A peridynamic investigation of dynamic brittle fracture[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF].

When inspecting damage profiles in a direction perpendicular to the main crack propagation for ∆U = 0.04 mm (see Fig. 7), it can clearly be observed that, initially, the damaged zone is close to the one-dimensional analytic damage profile d(y) = (1-|y|/(2l 0 )) 2 (A), thus explaining why Γ ≈ G c at the early stage of propagation. A fully damaged zone (d = 1) with increasing width then appears during crack advance (B and C). It is to be noted that the shape of the transition zone from d = 0 to d = 1 remains similar, at all stages, to the theoretical 1D damage profile. The associated increase in Γ before branching is, thus, essentially due to this wide zone with d = 1. At the onset of branching, one can observe that the central region exhibits a damage level lower than 1, the maxima of the damage profile occurring at two symmetric positions away from the central path, as if two single crack tips were propagating side-by-side (D). With a further reduction of the central zone damage level, the two tips are well isolated and start to mutually screen each other, leading to a symmetric deviation of the two branches from the horizontal direction. After branching (E), the damage profile of a single branch is similar to the initial profile with a small fully damaged band.

Let us highlight that, for all simulations, the crack propagates from left to right with an already fully established damage band width, meaning that this wide band does not appear after the crack tip has further propagated. Hence, these two modes of propagation (crack tip advance and damage band widening) occur simultaneously and, once the crack tip further advances, no further evolution of the damage band is observed in the crack tail.

Let us briefly discuss here the effect of the mesh size on the different damage profiles and the computed value of Γ . In these simulations, a mesh size of h = 0.02 mm has been used near the middle plane of crack propagation corresponding to 5 elements over a distance of l 0 . In [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF], it has been mentioned that the quasistatic damage profile can be wider than the theoretical one by approximately δ h ≈ 2h (see Fig. 8-top) leading to the following overesimated fracture energy:

Γ qs est = G c 1 + 3(2h) 8l 0 = G c + Γ mesh (16)
In our case, the quasi-static fracture energy would be over-estimated by 15%. In practice, a smaller overestimation (around 10%) is observed, see for instance the measured value of Γ ≈ 1.1G c in Fig. 6 for ∆U = 0.035 mm at the beginning of crack propagation. In the dynamic case, we observe a similar behavior i.e. the two sharp transitions from the fully damaged zone to the parabolic profile can be over-estimated by a distance h (see Fig. 8-bottom). As a result, if we consider a damage profile consisting of a fully damaged zone of width 2s due to dynamic effects, two mesh-related overestimations of size δ h ≈ h and the two parabolic profiles, then the effective dynamic fracture energy can be estimated as:

Γ dyn est = G c      1 + 6s 8l 0 dynamic + 6h 8l 0 mesh size      = Γ dyn + Γ mesh (17)
where Γ mesh ≈ 0.15G c . The previous expression of Γ dyn holds only in the case of a single-crack propagation (profiles A to C). In the transition regime (profile D), the decrease of damage at y = 0 should be taken into account. Except in one case (∆U = 0.04 mm), no evident decrease of the crack velocity has been observed after branching. Although such results are not reported here, simulations in different configurations (loading, boundary conditions) exhibited branching at much lower velocities than those reported here (around 0.2c R ). These results suggest that the crack velocity is not a determining parameter for crack branching. However, a critical value of the energy release rate (close to 2G c ) seems to be a necessary condition for branching for this phasefield model. After branching, the damage rate drops to a value close to G c , which is associated with a smaller damaged zone, until the process eventually repeats itself (multiple branching events have indeed been obtained for higher loading values). This mechanism is very similar to the one suggested in [START_REF] Ravi-Chandar | Dynamic fracture[END_REF].

In view of the present results, we believe that an energetic criterion is a pertinent way to describe the branching process and would, therefore, require further investigations, either from numerical simulations or experiments. In particular, it would be very interesting to design experiments with patches of different materials with a known fracture energy: one can think of a situation in which a single crack propagates in a tougher material and branches when arriving in a weaker one.

Velocity-toughening results

For all loadings, the evolution of the damage dissipation rate Γ as a function of the instantaneous crack velocity v has been recorded during the single crack propagation i.e. slightly after the initiation phase and before the occurrence of macroscopic branching. Both quantities have been represented in Figure 9. A well-defined master-curve is obtained for all loading levels showing that Γ is slightly greater than the quasi-static fracture energy Γ (0) = G c at low velocities whereas a strong increase of Γ can be observed when approaching a limit speed around 0.7c R . Remarkably, this value is close to the limiting velocity of 0.75c R identified experimentally in [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] for the same geometry and loading conditions.

The L-shape aspect of the curve is associated with a competition between two behaviours already identified in dynamic crack propagation experiments: at low velocities, a single crack propagating according to Griffith theory whereas, at high velocities, a constant velocity regime with important increase of additional dissipation mechanisms at the crack front.

This observation can be linked to the experimentally observed transition between a mirror-like propagation and a mist then hackle propagation characterized by surface roughening and microbranching. Although the present 2D model is not able to resolve the small length scales associated with the onset of microbranches (and obviously its 3D nature) or surface roughening, it is possible for the damage field evolution to "choose" (as a by-product of energy minimization) between dissipating more energy by advancing the crack tip at a higher velocity or by increasing the width of the damaged zone by propagating a damage front on a small distance perpendicularly to the main crack propagation (this last possibility could be associated with some kind of crack tip blunting).

The degree of generality of Figure 9 has been addressed by changing the material parameters to those close to Homalite-100 with E = 4.55 GPa, ν = 0.35, ρ = 1230 kg/m 3 and G c = 38.5 J/m 2 and to soda-lime glass with E = 72 GPa, ν = 0.22, ρ = 2440 kg/m 3 and G c = 3.8 J/m 2 while keeping the same geometry. In these cases, the Rayleigh wave speed is c R = 1093 m/s for Homalite and c R = 3102 m/s for glass. We also varied the PMMA plate geometry by considering a square Fig. 9 The damage dissipation rate Γ is a well-defined increasing function of the crack speed for all initial loadings. The points correspond to instantaneous values of v and Γ during the single crack propagation phase i.e. after initiation phase and before branching for high loadings. The dotted line corresponds to the experimentally observed limiting velocity in [START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF].

geometry of 32 mm × 32 mm and a slender band of 32 mm × 8 mm with the same pre-notch of 4 mm.

The Γ (v) relation is represented in Figure 10 for these different cases along with the collective results of Figure 9. It can be observed that the same increasing trend is obtained for each case. Interestingly, the non-dimensional results corresponding to Homalite-100 and glass superimpose quite well to those obtained with PMMA for the same geometry. On the contrary, it seems that the Γ (v) relation is more dependent on the geometry, at least at high velocities.

The origin of this dependence on geometry is not yet clear, it may certainly be related to a different interaction of the crack with stress waves reflected at the boundary of the domain. In any case, it is a signature of a difference in the conversion process of kinetic and elastic energies into fracture energy. It is interesting that such a process seems independent of material properties when represented in terms of non-dimensional variables while elastic moduli, fracture energies and sound wave speeds differ by several factors.

Interestingly, without including any rate-dependency in the model, our results are able to reproduce a velocitytoughening mechanism at the macroscopic scale. Although this increase of fracture surface seems to be limited to a factor 2.5 whereas in [START_REF] Sharon | Energy dissipation in dynamic fracture[END_REF] a factor 6 has been observed due to the creation of microbranches, our simulations show that a limiting velocity consistent with experimental observations can be achieved. For instance, this was not the case with rate-independent cohesive models [START_REF] Zhou | A rate-dependent cohesive model for simulating dynamic crack propagation in brittle materials[END_REF] in which cracks could accelerate up to the Rayleigh wave speed for this particular configuration.

Crack propagation in a heterogeneous medium

Brittle fracture mechanics in heterogeneous materials is now receiving increasing attention and more specifically in the dynamic case for which numerous questions are still unanswered. The role of heterogeneities in dynamic fracture processes is studied across various length scales, ranging from kilometres for earthquakes [START_REF] Fisher | Statistics of earthquakes in simple models of heterogeneous faults[END_REF][START_REF] Lapusta | Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip[END_REF], to millimetres for composite plates [START_REF] Yu | Three-dimensional modeling of intersonic shear-crack growth in asymmetrically loaded unidirectional composite plates[END_REF] or heterogeneous thin films [START_REF] Vasoya | Finite size and geometrical non-linear effects during crack pinning by heterogeneities: An analytical and experimental study[END_REF] or even to nanometres in metallic glasses [START_REF] Murali | Atomic scale fluctuations govern brittle fracture and cavitation behavior in metallic glasses[END_REF]. Heterogeneities can have a complex influence on a crack propagation by nucleating daughter cracks or exhibit a transition from weak to strong pinning regime [START_REF] Patinet | Quantitative prediction of effective toughness at random heterogeneous interfaces[END_REF]. In particular, all these various dissipative mechanisms in presence of heterogeneities make it difficult to predict the effective toughness of a heterogeneous material.

In this section, we revisit the limiting speed, branching and velocity-toughening mechanism in presence of heterogeneities. As mentioned previously, experiments have also reported that crack velocities near c R can be obtained when constraining the crack path to a weaker interface or by drilling an array of holes ahead of the crack [START_REF] Washabaugh | A reconciliation of dynamic crack velocity and Rayleigh wave speed in isotropic brittle solids[END_REF]. We first aim at reproducing numerically these observations. The velocity-toughening mechanism is then reinterpreted for a crack propagating along such an interface and subjected to unsuccessful branching at-Fig. 10 The normalized damage dissipation rate Γ/G c as a function of the normalized crack velocity v/c R does not seem to depend on material properties whereas it is dependent on geometry. Fig. 11 Constrained propagation along a weak interface: the same problem as in Section 3 is considered, except that the crack is forced to propagate along the weak interface of fracture energy G i c , the surrounding material has the same elastic properties but is much tougher so that macroscopic branching is prevented.

tempts. Finally, we also investigate how a distant heterogeneity can influence the crack propagation path.

Constrained crack propagation along a weak interface

First, a band of material of width t = 2l 0 = 0.2 mm and fracture energy G i c = 300 J/m 2 surrounded by a tougher material (G b c = 100G i c ) is considered (see Fig. 11). Elastic properties are still homogeneous and correspond to those of Section 3 so that the initial prestressed state is the same as before. In this configuration, the crack never branched (Fig. 12). Contrary to Fig. 7, the phase field is almost fully saturated to d = 1 with a steep transition to d = 0 at the interface between both materials. For a given loading, the crack is accelerating to higher velocities than for unconstrained configurations. Final velocities for different initial loadings are reported in Table 1.

In particular, it has been possible to obtain velocities up to 0.98c R for the highest loading. It is to be noted that for ∆U > 0.15 mm, the whole interface starts to damage at the same time due to the initial stress state exceeding the elastic limit of the damage model.

Fig. 13 Constrained crack propagation in a perforated material: here the same homogeneous material as in Section 3 is considered except that an array of holes has been drilled in front of the crack path. This results in an apparently weaker interface so that propagation is also favoured along this plane.

Crack propagation through a perforated material

We now go back to the case of the homogeneous material of Section 3 except that an array of holes of diameter D, spaced by a distance S, is considered on the mid-plane of the sample ahead of the crack tip (Fig. 13). We considered 30 holes with D = 0.4 mm and S = 0.9 mm so that the holes are sufficiently small for inducing only a negligible change of the total initially stored elastic energy and sufficiently weaken the middle interface. Although this configuration induces stress concentrations at the hole boundaries, the stress levels remain sufficiently low so that no damage occurs between holes prior to the crack tip arrival.

Crack propagation is driven by a succession of the following events (Fig. 14): when arriving near a hole, the crack tip is first attracted to the boundary; after a short time without any further evolution of the phase field, a new crack nucleates on the middle plane at the opposite point of the hole boundary and starts to propagate towards the next hole.

As for the homogeneous case, progressive widening of the damage band can be observed with increasing crack advance (events A to C); for sufficiently available energy at the crack tip, this can translate into an attempt at forming two branched cracks at the nucleation point (events C to E). The attraction exerted by the next hole can then stop the branching process and reform one single crack from the two branches. As the crack further propagates, the length of the side branches increases. Similarly to the weakened interface, higher crack speeds can be obtained in this configuration. For instance, an initial loading of ∆U = 0.05 mm leads to averaged velocities up to 0.9c R (Fig. 15(a)).

It is interesting to note that this particular simulation shares some qualitative similarities with the micro-branching phenomenon when considering the role of nucleation, growth and coalescence of microcracks ahead of the main crack path [START_REF] Ravi-Chandar | On the role of microcracks in the dynamic fracture of brittle materials[END_REF][START_REF] Rice | Physical Aspects of Fracture, chap. Some Studies of Crack Dynamics[END_REF][START_REF] Dalmas | Damage mechanisms in the dynamic fracture of nominally brittle polymers[END_REF]. This process can be further understood when looking at the evolution of the damage dissipation rate during crack propagation (Fig. 15(b)). The five representative events are also reported and correspond approximately to the moment when the crack tip is located halfway between two holes.

At first, the apparent damage dissipation rate is less than the medium fracture energy due to the weakening presence of the holes. The apparent porosity along the main crack path is φ = D/S ≈ 0.44, a simple crude estimate of the weakened interface fracture energy using a rule of "mixture" would then be G c,weak = (1 -φ)G c = 0.56G c . This estimate is roughly consistent with the measured value of Γ at the beginning of the propagation (after the initiation phase) in Fig. 15(b), which is around 0.7G c (event A). It is also interesting to note that the increase of Γ is relatively regular in the early stages of propagation, corresponding to a straight propagation of the crack through the first holes (up to event B). However, when the crack further propagates, bursts of increasing amplitude can be observed in the evolution of Γ , which can be directly associated with unsuccessful attempts at forming macrobranches (events C to E).

From a macroscopic point of view, the evolution of Γ during crack propagation can be related to the observations of Fig. 6 in the homogeneous case. One important difference is that the increase of Γ is more important since it goes approximately from 0.7G c to 2.3G c . Referring to Fig. 14 from a macroscopic point of view, this stronger relative increase in Γ can also be associated with a widening of a macroscopic damage band. This widening is here essentially caused by unsuccessful branching attempts, with increasingly longer side branches as the crack propagates. It may seem surprising that a weakened interface would lead to a stronger increase in Γ . However, we might argue that, although the initial crack propagation phase can be associated with an effective fracture energy G c,weak , the attractive presence of the holes makes it more difficult for a crack to branch. Assuming the existence of an energetic branching criterion like Γ ≥ G heter c,branch , we may suppose that G heter c,branch is essentially the same as for the homogeneous case (because branching occurs inside the bulk between two holes) or is even higher since the holes attraction tend to delay branching. Thus, in this heterogeneous case, Γ would increase from G c,weak < G c to G heter c,branch ≥ G homog c,branch . The presence of heterogeneities along the crack path has, thus, an important influence Fig. 14 Crack propagation in a perforated material with ∆U = 0.05 mm (the same colormap as in Fig. 5 is used). Macrobranching is temporarily suppressed by crack tip attraction to the holes. Small attempts of branching can be observed at the end (events D and E). See also Movie 1 in supplementary material.

on the velocity-toughening mechanism.

When looking at the evolution of the different forms of energy (elastic, kinetic and dissipated through damage propagation), it can be observed that less energy is dissipated into fracture surface (Fig. 16 for ∆U = 0.05 mm) when the crack path is constrained either in a weak interface or an array of holes than in the homogeneous case where macroscopic branching is possible. The excess amount is then mostly transferred into kinetic energy which can be linked to higher crack velocities.

As mentioned before, when approaching a hole, the crack tip is attracted to its boundary. This results in an increase of its instantaneous velocity before being momentarily stopped by the hole. In order to have a closer look at the interaction between the crack and a hole, we considered a simulation with 10 holes (D = 0.4 mm) separated by a distance S = 2.55 mm, the local increase of the crack tip velocity when approaching each hole can clearly be observed in Fig. 17. Associated with this increase of crack tip velocities, strong bursts of kinetic energy can also be observed just before the crack encounters a hole (inset of Fig. 17).

Note also that, for this configuration, the time interval separating the moment the crack stops at the hole boundary and nucleates again at the opposite point is measured to be ∆t = 0.694 µs, which corresponds exactly to the time for a Rayleigh wave to travel half of the hole circumference: ∆t = πD/(2c R ).

Crack propagation in presence of distant heterogeneities

We investigate the influence of distant heterogeneities on the crack propagation path. The same geometry for the PMMA plate is considered here with different hole configurations. We present only preliminary results in the case of holes, although other kind of heterogeneities such as stiffer inclusions have also been considered. Investigating the mechanisms driving the interaction between a crack and various heterogeneities will require a more thorough study, which will be the purpose of another work.

First, we considered one single hole of diameter D = 0.4 mm, the center of which is situated 0.6 mm away from the middle plane of the sample and either 1 mm or 6 mm away from the pre-notch tip (∆U = 0.04 mm). Recalling that cracks accelerate progressively until reaching a limiting speed (Fig. 4), the crack will pass near the hole situated 1 mm from the pre-notch tip with a smaller velocity than when passing near the one situated 6 mm away from the notch. The comparison of the different crack paths between the two configurations is reported in Figure 18. It can be observed that for the closest hole, the crack interacts with the hole at a slower velocity and is only slightly deviated towards the hole but then continues its straight propagation. In the other case the crack arrives near the hole with a higher velocity, a microbranch appears and is attracted then stopped by the hole. The main part of the crack then continues its propagation with a smaller damage band width after the microbranching event. In both cases, a decrease of the instantaneous crack tip velocity can be observed when interacting with the hole. The crack accelerates again after passing the hole. Let us mention that the opposite is observed in the case of a stiffer inclusion, the crack accelerates when passing near the inclusion and then decelerates.

These results are consistent with known results regarding attraction by holes and deflection by stiff inclusions. However, one may have expected that a fast crack will have less time to interact with a distant heterogeneity than a slower one, contrary to what is observed.

Fig. 16 Evolution of the different contributions to the total energy for ∆U = 0.05 mm (normalized by the initially stored elastic energy): elastic (blue), kinetic (green) and dissipated (red) for the unconstrained case (solid) and weak interface (dashed) and array of holes layout (dotted). Dissipated energy is less important in the case of a constrained crack path, the difference being mostly compensated by an increase in kinetic energy and only slightly by a drop of elastic energy. Fig. 17 Increase of the instantaneous crack tip velocity when approaching holes (gray bands) for ∆U = 0.05 mm. The dashed parts are due to the smoothing of the velocity when the crack is stopped by a hole. Inset: Rate of kinetic energy change as a function of time. Strong bursts of kinetic energy at a macroscopic level are associated with the crack accelerations when approaching a hole.

On the other hand, the velocity-toughening mechanism suggests that a faster crack seeks additional dissipation mechanisms than just a straight propagation at higher velocities. As the damage band widens, the distance at which it can interact with a defect may also increase. However, we do not have yet a clear explanation of these observations but we believe that it is a striking example of interesting interactions between fast cracks and heterogeneities, which therefore calls for further numerical or experimental investigations.

Finally, we have considered a configuration of 15 small holes of diameter D = 0.4 mm, the centers of which are situated either 0.5 mm or 0.6 mm away from the middle plane of the sample. The spacing between two consecutive holes is S = 1.95 mm. Two different loadings of ∆U = 0.04 and 0.05 mm have been considered.

The results for these configurations are reported in Figure 19. It can first be observed that the value of the considered offset has an important influence on the crack path patterns although the two geometric configurations are quite similar. In the first case of a 0.5 mm offset (Fig. 19(a) and (b)), the crack localizes in the weak plane composed by the holes when arriving near the first one. The propagation, is then similar to the one already described in the previous subsection with attempts at branching and attraction of the different branches by the next holes.

On the contrary, in the second configuration with a 0.6 mm offset and ∆U = 0.04 mm (Fig. 19(c)), the crack is slightly attracted by the presence of the first hole but then continues its straight propagation while small microbranches emerge from the main path towards the next holes. For ∆U = 0.05 mm (Fig. 19(d)), this behavior is more pronounced and the emergence of microbranches tend to deflect the main crack away from the holed plane. Additional microbranches can emerge and are attracted by the holes even when the main crack is quite far from the holed plane. As a result, the toughest plate corresponds to the second case with a 0.6 mm offset since the crack does not propagate through the weakened plane.

It is therefore quite surprising that, despite two configurations being almost similar, such different crack patterns are observed. Either the main crack is attracted by the hole and propagate along a weaker plane or it is repulsed by it and only microbranches are attracted by the holes, resulting in the main crack deflection.

Finally, let us highlight once more that simulations have also been conducted in the presence of heterogeneities consisting of a stiffer material. While a crack is generally attracted by the presence of a hole, the presence of stiffer inclusions tends to repel the crack. Besides, the presence of a heterogeneity consisting of the same elastic properties but with a different fracture energy is not able to influence a crack unless it is lo-cated on its path, contrary to heterogeneities consisting of different elastic properties.

Conclusions and perspectives

This work has investigated the capacities of the phasefield approach to reproduce specific features of dynamic crack propagation in brittle media. We have focused on the physical aspects which can be reproduced by numerical simulations using such a method. We considered a pre-strained PMMA plate configuration in which the initially stored energy is well defined and which leads to a progressive acceleration of the crack before reaching a steady-state regime. By removing the effect of stress waves induced by a suddenly applied loading in other numerical benchmark, we are able to better understand different aspects of dynamic crack propagation and the onset of branching. More precisely, various key results can be retained from this work:

-Crack propagation is characterized by a progressive widening of the damaged band width. It does not correspond to a later evolution of the crack tail due to delayed diffusion but occurs simultaneously to the crack tip advance. Besides, this wider band does not correspond to a larger regularization length but to a fully damaged zone with d = 1 with increasing width, the transition zone from d = 1 to d = 0 remaining similar to the 1D solution profile. -The damage band widening is associated with an increase of the apparent fracture energy Γ . A welldefined master-curve relating Γ to the crack velocity in a single crack propagation phase is obtained. This relation does not depend on material properties but seems geometry-dependent. -Macroscopic branching is observed when Γ reaches a sufficiently high value, corresponding approximately to 2G c . This observation favours an energetic criterion for branching. To our knowledge, this has not been remarked in previous works and would, therefore, require further investigations, both from a numerical and experimental point of view. -A limiting velocity around 0.7c R is observed, which is in accordance with experimental results for this specific configuration. Computations in constrained propagation along a weakened interface enabled to reach velocities close to c R , which has also been observed in previous experiments. -The considered phase-field approach is able to naturally account for velocity-toughening both at a macroscale with damage band thickening and branching and at a smaller scale when interacting with heterogeneities through microbranching. -Computations in presence of heterogeneities clearly showed the influence of defects on the crack propagation dynamics. In particular, holes on the crack path can delay macroscopic branching and lead to a stronger increase of the damage dissipation rate. Besides, the interaction of a crack with distant heterogeneities show that complex configurations can be obtained depending on the crack velocity and defect locations.

Regarding the last point, these preliminary simulations in heterogeneous media reveal a complex interplay between crack dynamics, material heterogeneities, microbranching and increase of fracture surface. The prediction of an effective toughness in dynamics seems to be a challenging question as various dissipative process can lead to a value higher than a simple average toughness. This question will require further investigation to, eventually, pave the way to the design of more efficient materials regarding dynamic fracture.

We believe that this work enabled to show that the phase-field approach is a legitimate candidate to study complex mechanisms of dynamic fracture. However, numerous questions remain still open. First, the pertinence of an energetically-based branching criterion should be more precisely assessed. Secondly, the present work did not manage to reproduce crack velocities obtained from the pre-strained PMMA plate experiments with the corresponding loading levels. Other models showed the same deficiencies and suggest rate effects have to be taken into account for PMMA. Finally, it is necessary to investigate the influence of 3D effects and the role of random material heterogeneities to better reproduce such experimental results.

Fig. 1

 1 Fig. 1 Continuum regularization of a possibly branched crack C by a continuous phase field d(x) ∈ [0; 1] in a continuum Ω subjected to imposed displacement and tractions T at its boundary. The regularized crack corresponds to the domain C d for which d > 0, its typical width is related to l 0 .

Fig. 2

 2 Fig. 2 One dimensional localized solution of an infinite bar in traction for the different damage laws. The total dissipated energy for this test case corresponds exactly to G c . For w(d) ∝ d (12), the damage profile is d = (1 -|x|/(2l 0 )) 2 for |x| ≤ 2l 0 and d = 0 for |x| ≥ 2l 0 . For w(d) ∝ d 2 (11), the damage profile is an exponential with infinite support d = exp(-|x|/l 0 ).

Fig. 3

 3 Fig.3Pre-strained PMMA plate problem (left) and crack patterns (right) sketched from experimental observations[START_REF] Zhou | Study on the macroscopic behavior and the microscopic process of dynamic crack propagation[END_REF] for different imposed displacements (taken from[START_REF] Wolff | A non-local continuum damage approach to model dynamic crack branching[END_REF]). For low loading levels, single crack propagation is observed. For higher values, microbranches appear during propagation. For even higher loadings, macroscopic branches are formed.

Fig. 4

 4 Fig.4Evolution of crack tip velocity for different loadings. For all loadings, the crack accelerates up to a limiting velocity around 0.68c R . Branching events are indicated by the star-shaped symbol.

Fig. 7

 7 Fig. 7 Symmetric vertical damage profiles for y ≥ 0 at different positions from the left part of the domain (∆U = 0.04 mm): A (x = 5 mm), B (x = 12 mm), C (x = 18 mm), D (x = 20 mm), E (x = 22 mm). The profiles clearly indicate the existence of a fully damaged zone (d = 1) with increasing width before splitting into two distinct cracks. The black square line corresponds to the 1D parabolic profile.

Fig. 8

 8 Fig.8Potential effect of mesh size (solid line) on the damage profiles (dashed lines without mesh-size effect): quasi-static case (top) and dynamic case (bottom).

Fig. 12

 12 Fig. 12 Crack propagation in a constrained path configuration with ∆U = 0.05 mm. Branching is completely suppressed, dissipation takes place only inside the middle interface and the crack accelerates up to the Rayleigh wave speed.

Fig. 15

 15 Fig.[START_REF] Cazes | Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture[END_REF] Crack velocity (a) and damage dissipation rate (b) evolution during propagation through an array of holes (∆U = 0.05 mm). Events A to E of Fig.14are reported usdashed lines. The damage dissipation rate increases with crack propagation and large bursts can be associated with unsuccessful branching attempts (events D and E). Velocities close to the Rayleigh wave speed are obtained.

Fig. 18

 18 Fig.18Close-up view of the interaction of a crack with a single hole at 1 mm (a) and 6 mm (b) from the pre-notch tip with an offset of 0.6 mm from the middle plane (∆U = 0.04 mm).

  Fig. 19 Close-up view of the interaction of a crack with an array of holes located 0.5 mm ((a) and (b)) or 0.6 mm ((c) and (d)) away from the middle plane for two different loadings: ∆U = 0.04 mm ((a) and (c)) or ∆U = 0.05 mm ((b) and (d)). See also Movies 2-5 in supplementary material.

Table 1

 1 Crack velocities for the weak interface configuration. Velocities up to 0.98c R have been attained.

	∆U (mm) Stored energy (J/m 2 ) Crack velocity (c R )
	0.04	618	0.81
	0.05	966	0.87
	0.10	3,863	0.94
	0.15	8,691	0.98

Acknowledgements The authors would like to acknowledge Corrado Maurini and Li Tianyi for sharing FEniCSbased implementation of damage-gradient models.