Kevin Moreau
email: kevin.moreau@ec-nantes.fr

Nicolas Moes

Didier Picart

Laurent Stainier

N Moës

Explicit dynamics with a non-local damage model using the thick level set approach

Keywords: thick level set, damage mechanics, explicit dynamics, impact

come

Explicit dynamics with a non-local damage model using the thick level set approach

INTRODUCTION

This paper focuses on the dynamical response of a material body modeled in continuum mechanics and subjected to impact loadings. Impact loadings are understood here as intense and brutal loadings that may lead to a rupture of the material body. Such loading conditions are encountered during crashes, explosions or any energetic and short-time phenomena. Numerically, models that predict the evolution of a material body up to the rupture are efficiently solved by explicit time integration schemes. These schemes proved to be robust and sufficiently accurate for engineering applications.

A lot of research efforts concentrate on modeling the rupture, and two major trends can be summarized: fracture mechanics and damage mechanics. On the one hand, fracture mechanics models propagation of large-scale defects (cracks) that are points (1D), lines (2D), or surfaces (3D) where unilateral discontinuity of kinematic fields occurs. But it fails to initiate new cracks. On the other hand, damage mechanics models an homogeneization of small-scale defects (micro-cracks) and handles both cracking and pulverisation (micro-cracking on a segment (1D), a surface (2D), or a volume (3D)). The continuity of a broken material is most of the time ambiguous. But it easily initiates new defects. This work takes place in damage mechanics.

Experimentally, the transition from a sane material body towards a broken material body is often limited to particular areas that experience intense strain. This phenomenon is called strain localization. It is shown that strain localization occurs whenever material exhibits strain softening, that is, a decline of uniaxial stress at increasing strain as defined by [START_REF] Bažant | Continuum theory for strain-softening[END_REF]. For a large number of materials such as steel, rock, soil, and concrete, localized areas present a characteristic length, which corresponds to the scale at which the localization phenomenon occurs.

For the class of local and rate-independent constitutive models presenting strain-softening (that is the consequence of a damage model or more generally the consequence of a plastic model with negative hardening) strain localization cannot be accurately modeled because there is no information about the localization scale. It is commonly observed numerically that for this class of models, localization tends (while refining the space discretization) to occur on a zone of null measure. A theoretical result may be found in [START_REF] Bažant | Wave propagation in a strain-softening bar: exact solution[END_REF] on a uni-dimensional material body in which localization occurs on a point. A consequence is that there is no energetic dissipation associated to the rupture. Numerically in the context of FEM discretization, a pathological mesh sensitivity is symptomatic of this lack of characteristic length.

Nevertheless, it seems that the class of local and rate-dependent constitutive models does give a more satisfactory solution [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF]. The characteristic time involved in such models induces a characteristic length by a space/time coupling due to the equation of dynamics. However, the characteristic time involved in rate-dependent model has nothing in common with the localization problem and leads to physically non acceptable localization areas. With the exception of the delay-damage model, where the characteristic time is intentionally introduced for this purpose [START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF]. But in this last model, an estimation of the resulting characteristic length given by [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF] shows that it depends on the state of stress inside the body and likely vanish when stress wave has small magnitude.

To regularize local constitutive models, the concept of non-locality is often used. There are several ways to introduce non-locality in the formulation of models going from a complete rewrite of governing equations as in microstructured continuum theory and related particular cases like Cosserat theory, second gradient theory, and Cosserat second gradient theory, see [START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies[END_REF], to a minimal rewrite of governing equations as in non-local integral models [START_REF] Bažant | Continuum theory for strain-softening[END_REF][START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] or gradient-enhanced models [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]. A comparison of non-local models can be found in [START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF] and concludes that few formulations become mathematically well-posed after non-local regularization. It also concludes that for mathematically well-posed formulations, such models lack of thermodynamic grounds. It seems, however, that this does not affect the effectivity of the regularization. Another comparison between non-local integral model and delay-damage model is proposed in a dynamic context in [START_REF] Desmorat | Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading[END_REF]. Classifications of non-local approaches can be found in [START_REF] Chambon | Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies[END_REF][START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF][START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF].

We can mention as well phase field approaches [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Miehe | A phase-field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] emanating from fracture mechanics, which regularize the discontinuity introduced by the crack with a smooth auxiliary phase field. Then, the fracture energy is written in terms of the phase field and a smooth loss of some material properties with respect to the phase field is added. These approaches lead to a formulation close to a continuum damage mechanics formulation for which damage variable would be smoothed. Extension to dynamic brittle fracture can be found in [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF].

We now focus on the Thick Level Set (TLS) approach, a young approach to introduce non-locality in the formulation. It has been first introduced in a quasi-static context [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF][START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF][START_REF] Stolz | A new model of damage : a moving thick layer approach[END_REF] and compared with phase field approaches [START_REF] Cazes | Comparison of a phase-field model and of a thick level set model for brittle and quasi-brittle fracture[END_REF] and is presented here in a dynamic context. This approach has the particularity to provide the non-locality only when a first defect emerges in the material body. More precisely, the computation is purely local until damage initiates, for instance, when the local energy release rate (the thermodynamic force conjugate to damage rate) reaches a threshold value. Then, non-locality starts acting on the formulation by means of a non-local energy release rate and a non-local damage rate, which are responsible of the damage growth. This newly appeared nonlocality is also limited in space to damaged areas. When damage grows enough, the material body presents fully damaged areas and the non-locality stops acting on it. This fully damaged areas are dumped from the computation, and a discontinuity appears in kinematic fields. To sum up, the nonlocality appears, grows, and stops with damage. The two non-local fields are the average of their local counterparts but in a different way than that for non-local integral models. All this features are handled by an auxiliary field: a level set function field. This level set function is chosen as a signed distance function so that we have a way to introduce a characteristic length. The iso-zero curve of this level set function is also coincident with the damage front so that we are able to easily discretize the interface between the sane material and the damaged material. The damage growth problem is transformed on a level set function propagation problem, and as a consequence, all the work of the non-locality is to provide an appropriate level set function rate for the level set function to be updated. In this paper, we propose a new way to compute non-local fields based on the construction of an approximation space that has all required properties to compute the average of local counterparts.

The paper is organized as follows: in the second section, we recall governing equations of classical continuum damage mechanics with a constitutive model written by means of a free-energy potential and a dissipation potential. We introduce the TLS approach and we define two non-local fields, the non-local energy release rate, and the non-local damage rate. Once every governing equations are given, we give in the third section a time discretization and we write three variational problems that must be solved at each time step; the third one is the weak form of the balance of momentum equations, and the two others are related to the TLS approach. They permit to compute non-local fields and are limited to the localization region. Then, in the fourth section, we give a space discretization and we present the new way to compute non-local fields. The fifth section is dedicated to numerical applications. Finally, the sixth section proposes a discussion of the TLS approach exposed in this paper.

GOVERNING EQUATIONS

Mechanical problem

We assume a material body as a continuum of homogeneous density whose material points belong to the domain R 2 and study its evolution during the time interval T D OE0I t f R C ,where t f is the final time, in terms of displacement u and damage d . Both fields are supposed continuous in space and time. We place ourselves under the small strain assumption. We provide a constitutive model, formally denoted by functionals F D F .".u/; d / and F d D F d .".u/; d /,where" is the small strain tensor. We impose three initial conditions: the initial displacement u 0 , the initial velocity P u 0 , and the initial damage d 0 . We also impose two boundary conditions: the Dirichlet boundary condition u d and the Neumann boundary condition T d , respectively, on parts @ D and @ N of the domain frontier such that @ D [ @ N D @ and @ D \ @ N D; . It leads to the study of the following initial boundary value problem:

R u Dr D > in T D F .".u/; d / P d D F d .".u/; d / (1) of initial conditions u D u 0 P u D P u 0 in ¹0º d D d 0 (2)
and boundary conditions

u D u d on @ D T n D T d on @ N T (3)
where is the Cauchy stress tensor and n is the outward normal. The aims of the TLS is to provide a constitutive model that is consistent with regard to the localization problem. The constitutive model is built by defining both a free-energy potential ' and a dissipation potential , and applying on them normality rules to derive, respectively, state and evolution laws. We will see in the following that in the TLS the state laws remain unchanged compared with classical thermodynamic framework whereas the evolution law introduces non-locality in the formulation.

Free-energy potential, local state laws

We define the small strain tensor

" D r s u D 1 2 r u C r u > in T( 4 )
and the local energy release rate Y . The free-energy potential ' D ' .";d/ provides relations between ."; / and .d; Y / thermodynamic pairs. Thermodynamics stability requires that the freeenergy is a lower semi-continuous and convex functional of small strain and damage [START_REF] Germain | Continuum thermodynamics[END_REF]. Then, state laws may be derived by calling normality rules:

. ; Y / 2 @' .";d/ in T [START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF] where @' stands for subdifferential. For a differentiable free-energy potential, it simplifies to

D @' .";d/ @" in T( 6 ) Y D @' .";d/ @d in T( 7 )
From state laws, we derive the energetic dissipation by using Clausius-Duhem inequality:

W P " P ' D Y P d > 0 in T( 8 )
The fulfilment of this inequality is ensured by introducing a dissipation potential. Nevertheless, it is known that if we provide a local evolution model to relate the damage rate to the energy release rate, the obtained formulation will fail to model accurately the localization problem. A solution is to introduce a characteristic length in the formulation to add non-locality.

Dissipation potential, non-local evolution law

In order to introduce a characteristic length, we use the TLS approach first described in [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF]. We assume a continuous level set function that takes value zero on the damage front so that we can define it as

.t/ D¹x 2 j .x;t/ D 0º (9)
This front is a closed interface or is in intersection with domain boundary. It splits domain into two regions: one composed of sane material and the other one composed of damaged material. We further assume that the level set function verifies the property:

kr k D 1 in T ( 10
)
This equation is a particular case of the stationary eikonal equation and the associated boundary value problem admits many distinct solutions as highlighted in [START_REF] Crandall | Viscosity solutions of Hamilton-Jacobi equations[END_REF] for the Dirichlet boundary value problem. The level set function that we are looking for is given by the fast marching methods [START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF] and often called the signed distance function to the damage front. For simple geometries, the signed distance function is given by

.x;t/ D˙min y2.t/ kx yk (11)
In the following numerical applications, Equation ( 11) is used instead of the fast marching method. The use of the fast marching method is a perspective of this work. We still call the level set function in the following. With this function at hand, we map the damage field on it in the following way:

d./ D 0 6 0 d 0 ./ > 0 0 < < l c (12) d./ D 1 > l c
where l c is the characteristic length and the prime operator 0 denotes the derivative with respect to the level set function. We call this function the damage profile. Note that the introduction of a characteristic length is possible because the level set function has the unit of a length. The damage growth problem is thus transformed on a level set function propagation problem:

² P d D F d .".u/; d / in T d D d 0 in ¹0º ! ² P D F .".u/; / in T D 0 in ¹t i º (13)
where we denote formally F D F .".u/; / as an evolution model for the level set function and where t i 2 T is the time at which material body begins to suffer damage. Let us remark that when the initial damage d 0 is zero, there is no damage front and therefore it is somewhat hard to define the initial level set function 0 . In fact, as the damage is not initiated, non-locality is not acting and the computation is purely local. As soon as damage initiates, we place a damage front by considering a level set function 0 . In practice, we place a circular iso-0 with radius l i <l c and centered on the material point that suffers damage. This implies a small perturbation of the damage field but is not inconvenient in practice. This is the way the TLS approach handles transition from a local formulation to a non-local formulation.

With the representation (12) of the damage field, the damage rate is given by

P d D d 0 ./ P in T ( 14
)
An important consequence is that to keep the property [START_REF] Desmorat | Delay-active damage versus non-local enhancement for anisotropic damage dynamics computations with alternated loading[END_REF], the level set function rate P verifies

d dt kr k D r kr k r P D 0 in T (15)
A simple manner to manipulate fields that verify such a constraint is to define a new coordinate system: x D .x; y/ ! ´D .´; s/ attached to the damage front where ´follows the level set function gradients and s follows the iso-curves. The iso-0 curve (the damage front) defines ´D 0. And s D 0 is taken arbitrary along the iso-0. The jacobian of the change of coordinate is

J.´;t/ D 1 .s; t / (16)
where is the radius of curvature of the iso-0 curve. The jacobian is valid when the damage front is regular. The damage front is supposed regular with the exception of the first application where a specific comment is addressed. Note that the coordinate system itself is depending on time variable because it is attached to a moving front. Note as well that in practice, we never construct such a coordinate system. Then, property [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF] sure that within a one-element thick fully damaged area, the discontinuity is naturally (understood as no need to remesh) handled by the approach. It provides a transition from damage to fracture.

In the TLS approach, there is no straightforward relation between the damage rate and the local energy release rate. Instead, Equation ( 17) is used, and an evolution model is provided for the level set function. To this purpose, we introduce a pair of thermodynamically conjugate non-local variables N Y; N P d called, respectively, the non-local energy release rate and the non-local damage rate, Figure 2. A damage band, the curvilinear coordinate system attached to it, the length l.s; t/ (the upper bound of the interval we integrate over in [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF] and ( 24)) and the local and non-local energy release rate Y and N Y .

which are defined on C only and verify property [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF], that is, N Y D N Y.s;t/and N P d D N P d.s; t/ and we want the pair to verify the additional property:

Z C N Y N P d dV D Z C Y P d dV in T (22)
meaning that the global dissipation due to non-local thermodynamic pair is the same as the local one. The Clausius-Duhem equation is therefore verified at least globally. The expression of the nonlocal energy release rate has already been given in [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF], it is resulting from the stationarity of the potential energy functional (that makes sense in quasi-statics) with respect to a variation of the level set function: These non-local fields appear to be some average of their local counterpart in the direction of the level set function gradient. Note that definition of the non-local energy release rate introduces an additional weighting function through the first derivative of the damage profile. An evolution model is written in terms of non-local quantities by calling the normality rule on a dissipation potential D N Y . To preserve thermodynamics stability, this potential must be a positive, lower semicontinuous, and convex functional of non-local energy release rate, null for null non-local energy release rate [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF]. Then, the evolution model is given as

N Y.s;t/D R l
N P d 2 @ N Y in C T ( 25
)
For a differentiable dissipation potential, it simplifies in

N P d D @ @ N Y N Y in C T (26)
Equation ( 17) is introduced inside Equation ( 24), and we use notation

N d 0 .s; t / D R l.s;t/ 0 d 0 ..´; t //J.´;t/dŔ l.s;t/ 0 J.´;t/d´( 27
)
to obtain an expression of the level set function rate

P .s; t/ ˇC D N P d.s; t/ N d 0 .s; t / (28)
So far, the damage evolution model is not prescribed. The only constraints on the level set velocity are given by Equation ( 15) and the fulfilment of Clausius-Duhem equation at the global stage with an appropriate choice of the dissipation potential. The level set function rate is known inside C and has to be extended to the other side of the damage front. In level set methods, see [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF], the level set function must be updated at least on a narrow band defined around the iso-0. We introduce a narrow band N for which @ N and we apply the following advection boundary problem:

r kr k r P j N D 0 in N T P j N D P j C on T (29)
which completes the governing equations needed to compute the formal problem given in [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF]. We now discretize in time governing equations and give variational problems that have to be solved at each time step. Note that from ( 28) and ( 29), the level set function can only be updated in C [ N .

As commonly carried out in level set methods, a reinitialization step is added to update the level set function on all the domain . In numerical simulations, we use [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] for that purpose.

TIME DISCRETIZATION

The time discretization is built by considering a growing sequence of time values

t k 2 T, k 2 N.
We define time steps as t D t kC1 t k and we simplify notations from f x;t k to f k for clarity. The main unknowns are the displacement field u.x;t/and the level set function .x;t/.W e have to integrate in time two partial differential equations: the second-order balance of momentum equations, given in Equation ( 1), and the first-order level set propagation equation, given in Equation [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF]. We apply the second-order accurate central difference time integration scheme to the first one and the first-order accurate forward Euler method to the second one, leading to

u kC1 D u k C t P u k C 1 2 t 2 R u k (30) kC1 D k C tF " u kC1 ; k (31)
R u kC1 D 1 r F " u kC1 ;d kC1 (32)
P u kC1 D P u k C 1 2 t R u k C R u kC1 (33)
This is the most straightforward way to discretize in time the TLS approach. Both schemes are explicit, they remain stable if time steps are upper bounded by a critical value t c .I ti sg i v e ni n the next section after a space discretization. Note that the obtained scheme is fully explicit, because the evolution model itself is computed from k . Knowing the displacement field u kC1 and the level set function k , we deduce the strain field " kC1 from (4) and the damage field d k from ( 12) and derive the local energy release rate Y kC1 by applying state law [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]. We then compute non-local fields N Y kC1 from ( 23) and N P d kC1 from ( 26) and deduce the level set function rate field P kC1 from (28) only in domain C k . We then extend the level set function rate to the other side of the damage front in the narrow band N by applying [START_REF] Frey | Mesh Generation: Application to Finite Elements[END_REF]. We deduce kC1 by updating the damage front and applying [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF]. Computations of the non-local energy release rate, and the level set function rate are not straightforward because they are defined on a curvilinear coordinate system attached to the damage front. One way to compute them is to establish variational formulations, as proposed in [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF], that we discuss in the next section.

Variational formulations

Following [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF], we introduce the constrained space

N S D ´y regular on C k ˇr k r k r y D 0 μ (34)
Any function lying on this function space is constant along level set function gradients. In other words, these functions depend only on coordinate s of the curvilinear coordinate system attached to the damage front. Then, we apply the two following variational problems:

Problem 1. knowing k and Y kC1 both regular on C k ,find N Y kC1 2 N S such that Z C k N Y kC1 Y kC1 d 0 k dV D 0 8 2 N S ( 35
)
Problem 2. knowing k regular on C k and N P d kC1 2 N S,find P kC1 2 N S such that

Z C k d 0 k P kC1 N P d kC1 dV D 0 8 2 N S ( 36
)
The third variational problem is given by

Problem 3. knowing u kC1 2 U , kC1 , T kC1 d , u kC1 d and R u kC1 d ,find R u kC1 2 R U such that Z R u kC1 u dV C Z F " u kC1 ;d kC1 W r s u dV Z @ N T kC1 d u dS D 0 8u 2 U 0 (37)
in spaces

U D °u regular in j u D u kC1 d on @ D ± (38)
R U D °R u regular in j R u D R u kC1 d on @ D ± (39)
U 0 D ¹u regular in j u D 0 on @ D º (40)
The difficulty lies in the computation of variational formulations in the space N S. In the next section, we build an approximation space that weakly verifies the constraint.

SPACE DISCRETIZATION

This section is dedicated to the computation of Problem 1, Problem 2, and Problem 3. Problem 3 is standard in the FEM, the only difference is the numerical integration that takes into account the damage front. For Problem 1 and Problem 2, the difficulty lies in the constrained function space N S. A solution adopted in [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF] was to simplify the space by removing its constraint and add more difficulty to the variational equation. This was carried out by using lagrangian techniques. Here, we propose a new approach in which we build a finite dimensional function basis where every function weakly verifies the constraint. The new approach can be summarized as follow: first, we build some approximation functions on the damage front that we call modes. These are obtained by well chosen linear combinations of FEM approximation functions intersecting the damage front. Then, we advect these modes to the damage band by solving several times an advection equation where the advective velocity is the level set function gradient and modes are given as Dirichlet boundary conditions. We end up with extended modes and we use them as basis functions of a finite dimensional constrained function space N S h . Note that these advection problems are solved in an efficient way, compatible with an explicit dynamics computation. We first explain how to build modes on the damage front and then explain how to extend it to the damage band. Then, we generalize such an extension to a narrow band on the other side of the damage front to be able to update the level set function.

Kinematic fields

We consider a space discretization h of domain and we approximate kinematic fields using FEM piecewise linear approximation functions defined on three-noded triangular element (we omit the notation h for clarity):

u.x/ D X i 2S D u i N i .x/ (41)
P u.x/ D X i 2S D P u i N i .x/ (42) R u.x/ D X i 2S D R u i N i .x/ (43)
where a DOF and its associated approximation function are numbered with a unique number (ID) belonging to the set S D .Therearen D IDs in S D . Applying the space discretization to the variational formulation given in Problem 3, we write the problem using matrix notations

OEM 2 M n D , ¹F º2 R n D , ¹Rº2R n D that is: find ¹ R U º2R n D such that OEM ¹ R U ºC¹F ºD¹Rº (44)
where we have

M ij D Z N i N j dV (45)
F i D Z W r s N i dV (46) R i D Z @ N T d N i dS (47)
With the discretization at hand, we can express the critical time step of the whole scheme used in numerical simulations. For the central difference scheme alone, we use the lower estimate critical time step provided by linear elastodynamics, see [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF],

t u c D h c l (48)
where c l is the longitudinal elastic wave velocity and h is the smallest element characteristic length.

For an element undergoing damage, the wave velocity is altered, leading to a greater critical time step. For the forward Euler method alone, we first consider the level set propagation equation

P v kr k D 0 in ( 49
)
where v is the physical velocity imposed to the level set function. From ( 10) and ( 13), we express this velocity as

v D P D F .".u/; / (50)
The critical time step is given as, see [START_REF] Osher | Level Set Methods and Dynamic Implicit Surfaces[END_REF],

t c D h max x2 .jv j/ (51)
Finally, the critical time step for the whole scheme is chosen as

t c D min t u c ;t c (52)
Without going into too much detail, the numerical integration takes into account the damage front. In a cut element, damage is strictly positive on the C side, and zero on the other side. The mass matrix is lumped as usual in explicit dynamics.

4.2.

Non-local energy release rate, level set function rate 4.2.1. Modes on damage front. We use the previous space-discretization h to build a nonconforming space discretization C h of domain C and a space discretization h of damage front . Space discretization h is built by finding intersections between the damage front and some C h edges. We define on each vertices belonging to C h a DOF. All IDs are collected in the set S d and we define hat functions N i ;i 2 S d . We then use algorithm described in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] to select some vertices (called vital vertices) belonging to h . For each vital vertex, we define a DOF. All IDs are collected in another set S and we define basis functions ' ;˛2 S , called modes. The DOF of a vital vertex is linked to DOFs of some vertices belonging to C h . All the corresponding IDs belong to a last set S l S d . See Figure 3 for an illustration of all different sets. Then, modes are defined as linear combinations of hat functions N i ;i 2 S l and restricted to h (we omit the notation h for clarity):

' .x/ D X j 2S l a j N j .x/j ˛2 S ;a j 2 R (53)
These modes are well chosen to impose stiff interface conditions within the eXtended Finite Element Method (X-FEM), see [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF]. But this is not necessarily a good choice for a front propagation. Indeed, depending on the location of the damage front with respect to the space discretization C h , the 'size' of modes may present a large disparity and some mode can be very thin. To keep a front propagation smooth enough, we propose to remove modes that are considered too thin by combining them with neighboring ones. To this purpose, we consider the front length

l D Z dS (54)
Because the mode basis is a partition of unity, see once more [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF], we have The length is obtained as the sum of terms l ˛of a lumped mass matrix with unit density assembled on the damage front. Matrix entries are good quantities to evaluate mode sizes. Practically, we define a threshold value l red under which the mode is eliminated and combined to its neighbors. Starting from n modes ' ;˛2 S , we end up with n red filtered modes Q ' ;˛2 S red ,w h e r e S red S is a set of DOFs IDs. We add a tilde notation on filtered modes. To build filtered modes, we start from modes and recombine them based on Algorithm 1. Removing a mode ˛affects mode ˛ 1 and mode ˛C 1

l D X ˛2S Z ' dS " ƒ' … l ˛ (
Q ' ˛ 1 Q ' ˛ 1 C 1 2 Q ' (56) Q ' ˛C1 Q ' ˛C1 C 1 2 Q ' (57)
Finally, we get

Q ' .x/ D X j 2S l Q a j N j .x/j ˛2 S red ; Q a j 2 R ( 58
)
Mode filtering is basically a rearrangement of a j ;˛2 S coefficients to Q a j ;˛2 S red coefficients. remove mode by combining it to neighbors using (56) and (57), set smallest term of mass matrix in l, 10: end while 4.2.2. Modes advection to damage band. Now that several modes are defined on the damage front, we advect them in the damage band by applying n red times the following advection boundary value problem

r kr k r Q ' ˛D 0 in C ˛2 S red (59) Q ' ˛DQ ' on 2 S red ( 60
)
where the advection velocity is the level set function gradient and modes Q ' are imposed as Dirichlet boundary conditions. We recognize in this problem the constraint [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF] on the level set function rate and non-local fields. It means that these extended modes are good candidates to build function space N S. One can note that on , the inequality r kr k n <0is verified, because n D r kr k , as required for such a purely advective boundary value problem; see [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF]. This problem can be satisfactory solved by applying a Streamline-Upwind Petrov-Galerkin (SUPG) weak formulation as established in [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], that is, with

Problem 4. for ˛2 S red ,find Q ' ˛2 S ˛such that Z C r kr k r Q ' ˛ ' dV C Z C N r kr k r Q ' ˛ r kr k r ' dV D 0 8' 2 S 0 (61)
S ˛D ® ' 2 H 1 C ˇ' DQ ' on ¯(62) S 0 D ® ' 2 H 1 C ˇ' D 0 on ¯(63)
where N is a diffusion term. Using the TLS approach, the damage front is in general not conforming to the mesh and the spatial integration must take into account only parts of finite element that are inside the domain C ; see Figure 4 for an illustration. More importantly, Dirichlet boundary condition are harder to prescribe on a damage front that is not mesh conforming. A solution is to combine the SUPG formulation to a weak imposition of Dirichlet boundary condition by means of a Lagrange multiplier defined on a stable Lagrange multiplier space; see [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] for a precise definition and motivation on the use of such a space. The weak formulation becomes Problem 5. for ˛2 S red ,find Q ' ˛2 H 1 C and ˛2 H 1=2 ./ 0 such that

8 < : Z C r kr k r Q ' ˛ ' dV C Z C N r kr k r Q ' ˛ r kr k r ' dV C Z ˛' dS D 0 8' 2 H 1 C Z Q ' ˛ Q ' dS D 0 8 2 H 1=2 ./ 0 (64)
Note that the Lagrange multiplier space must be defined as proposed in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF], otherwise, the problem fails to pass the inf-sup condition. To update the damage front position, the level set function rate must be defined on both side of it, as given by Equation [START_REF] Frey | Mesh Generation: Application to Finite Elements[END_REF]. In practice, the level set function rate is only extended to a narrow band N defined on the other side of the damage front with respect to the damage band, as depicted in Figure 4 and problem (29) similar to mode extension is solved on the narrow band. Extended modes are approximated using FEM piecewise linear approximation functions defined on three-noded triangular element and Lagrange multipliers are approximated using non-filtered modes

Q ' ˛.x/ D X j 2S d b j N j .x/2 S red ;b j 2 R (65) ˛.x/ D X j 2S c j˛'j .x/2 S red ;c j˛2 R (66)
Variational formulations given in Problem 1 and Problem 2 are then discretized by using the extended modes as basis functions, that is, trial and test fields belong to the following finite dimensional function space:

N S D 8
n d ® n ® n d ' …" ƒ A C B n '…"ƒ C > C0 n red ' …" ƒ D 1 D n red E 1 E n red D n red ' …" ƒ 0 0 D 1 D n red ( 68
)
where

OEA 2 M n d , OEB 2 M n d , OEC 2 M n n d and ® D ¯2 R n with A ij D Z C N i r kr k r N j dV (69) B ij D Z C N r kr k r N i r kr k r N j dV (70)
C i j D Z ' i N j dS (71) D ˛i D Z Q ' ' i dS (72)
Now, we can build the extended modes from (65) and the finite dimensional function space N S from (67). Variational formulations given in Problem 1 and Problem 2 reduces to find

® N Y ¯R 2 R n red and ® P ¯R 2 R n red such that OEM R ® N Y ¯R D¹Y º R (73) OEM R ® P ¯R D °N P d ± R ( 74
)
where

OEM R 2 M n red , ¹Y º R 2 R n red and °N P d ± R 2 R n red with M ˛ˇD Z C Q ' ˛Q ' ˇd 0 ./ dV (75) Y ˛D Z C Y Q ' ˛d 0 ./ dV (76) N P d ˛D Z C N P d Q ' ˛dV (77)
We introduce matrix OEP D b j 2 M n d n red containing all DOF values of all extended modes resulting from ¹D ˛º vectors concatenation. Therefore, we can rephrase Problem 2 and Problem 3 using classical FEM basis functions

OEM R D OEP > OEM OEP (78) ¹Y º R D OEP > ¹Y º (79) °N P d ± R D OEP > °N P d ± (80) ® N Y ¯D OEP ® N Y ¯R (81) ® P ¯D OEP ® P ¯R ( 82
)
where

OEM 2 M n d , ¹Y º2R n d , °N P d ± 2 R n d , ® N Y ¯2 R n d and ® P ¯2 R n d with M ij D Z C N i N j d 0 ./ dV (83) Y i D Z C YN i d 0 ./ dV (84) N P d i D Z C N P dN i dV (85)
So, we finally solve the two problems in the following way: find ® N Y ¯and ® P ¯such that

OEP > OEM OEP ® N Y ¯R D OEP > ¹Y º ® N Y ¯D OEP ® N Y ¯R (86) OEP > OEM OEP ® P ¯R D OEP > °N P d ± ® P ¯D OEP ® P ¯R (87)
The mass matrix can be lumped to improve efficiency, we define ¹M º2R n d and ¹M º R 2 R n red with

M i D X j 2S d M ij (88) M ˛D X ˇ2S red M ˛ˇ(89)
We have the relation

¹M º R D OEP > ¹M º (90)
So that the two problems become the following: find

® N Y ¯and ® P ¯such that OEP > ¹M º ® N Y ¯R D OEP > ¹Y º ® N Y ¯D OEP ® N Y ¯R (91) OEP > OEM ® P ¯R D OEP > °N P d ± ® P ¯D OEP ® P ¯R ( 92
)
where operator denotes the element-wise product. We apply a linear system similar to linear system (68) to compute matrix OEP N that is the same operator as OEP but on the narrow band. We deduce the level set function rate on the narrow band simply by computing

® P ¯N D OEP N ® P ¯R (93)
All the computation process is resumed in the Algorithm 2 for the case of a lumped mass matrix. We consider this process in numerical applications.

Algorithm 2 Computation of the level set function rate in domains C and N . = denotes the element-wise division. 1: compute OEP , 2: assemble lumped FEM mass matrix ¹M º, 3: assemble right and side ¹Y º,

4: compute ® N Y ¯! OEP OEP > ¹Y º = OEP > ¹M º , 5: compute °N P d ± ! @ ® N Y ¯ =@ N Y , 6: compute ® P ¯R ! OEP > °N P d ± = OEP > ¹M º , 7: compute ® P ¯! OEP ® P ¯R, 8: compute OEP N , 9: compute ® P ¯N ! OEP N ® P ¯R,
By solving (91), ( 92), (93), and (44) while using the two explicit time integration schemes, one may solve the following problems: Kalthoff and Winkler experiments and the single-edge notched tension test that are proposed in the following section.

Mesh adaptation

In this work, there are two length scales: the characteristic length of the structure and the characteristic length of the material. When the last one is much smaller than the first one, a uniform mesh often leads to mesh considered as too thin. While accuracy is increased, so do computational costs.

A natural solution is to use an adaptive mesh to preserve accuracy in areas of interest while decreasing computational costs by coarsening the mesh everywhere else. For this purpose, we use a time adaptive mesh refined around the damage band.

It is based on a balanced quadtree data structure with hanging nodes (incompatibly placed nodes) [START_REF] Frey | Mesh Generation: Application to Finite Elements[END_REF][START_REF] Krysl | Natural hierarchical refinement for finite element methods[END_REF][START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF] combined with a straightforward projection of displacement, velocity, and level set function fields. The mesh is obtained from the quadtree data structure by assigning two three-noded triangular elements per cell. This division step is alternated for two neighboring cells. A refinement level is assigned to cells and elements to identify their sharpness. The unique cell of refinement level 0 is the bounding box of domain . Others are obtained by recursive and conditional subdivisions of a cell into four cells of one higher refinement level. Domain boundary @ conforms to the mesh so that elements outside domain are simply filtered out. Damage front is given by the iso-0 of the level set function and is treated within the X-FEM [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finiteelement method[END_REF]. Hanging nodes specific treatment is detailed in [START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF].

The mesh is adapted in time using refinement only, that is, the number of DOFs only increases. Initial refinement is specific to numerical applications. Next, refinements are then based on three criteria. The first criterion concerns the sign of the level set function. If the level set function is positive, then the highest refinement level (chosen as a numerical parameter) is used. If the level set function is negative, then the second criterion is used. This second criterion accounts for the absolute value of the level set function (the distance to the front). If this value is lower than a certain length l m , the highest refinement level is used, otherwise, the refinement level decreases as a function of jj l m . It takes into account the balance of the quadtree data structure so that a single hanging node per edge is possible. The last criterion consist in refining elements for which the refinement level is too low to a minimum refinement level (also chosen as a numerical parameter). Based on the two first criteria, the adaptive mesh is closely related to the level set function field and values on new nodes are simply interpolated.

For the displacement and velocity fields, the projection is realised by solving the two following problems:

Problem 6. knowing 1 u 2 V,find 2 u 2 V such that Z 2 u 2 u dV D Z 1 u 2 u dV 8 2 u 2 V (94) Problem 7. knowing 1 P u 2 V,find 2 P u 2 V such that Z 2 P u 2 P u dV D Z 1 P u 2 P u dV 8 2 P u 2 V ( 95
)
where

V D ¹u regular in º (96)
and we used left superscript 1 and 2 to specify the mesh a field is discretized on: 1 for old mesh, 2 for new mesh.

The mesh is adapted when the damage band advances on a distance ˛lm ;˛2 R so that it always stays inside the highest refinement level area. Note that this solution is preferred over a fully dynamic mesh adaptation as in [START_REF] Sochnikov | Level set calculations of the evolution of boundaries on a dynamically adaptive grid[END_REF] because mass matrix assembly and projection step are not fully optimized from an implementation point of view. This projection step takes place after Equation (33).

NUMERICAL APPLICATIONS

V-shaped damage band

In this section, we illustrate the non-local computation based on modes on a V-shaped damage band by computing the non-local energy release rate N Y . We define a domain R 2 as a unit length square and a damage band as depicted in Figure 5. We assume an energy release rate field given by Y. (99)

The prescribed damage front is not regular and the jacobian given in Equation ( 16) is singular for x D .1I 1=2/. Nevertheless, the computation of Equation ( 23) on a small variation of the damage front shape to get a circular tip of radius " is possible and leads to

N Y D 1 l c R " 0 .Y C ´/ 1 " d1 l c R " 0 1 " d´( 100) D Y C 1 3 " (101)
where D cos arctan 1 2 , by taking the limit when " tends toward 0, the non-local energy release rate at the singular tip is found equal to the energy release rate. The solution is discontinuous, and the discontinuity is due to that of the level set function gradient; see Figure 2. We are interested in the capacity of the method to approximate such a solution. We therefore study its convergence with respect to space discretization. To do so, computation is carried on five cartesian meshes going from 16 16 to 256 256 and made of three-noded triangular and alternated elements for which the element characteristic length are h D 6:25010 2 I 3:12510 2 I 1:56310 2 I 7:81310 3 I 3:906 10 3 m. The method possesses two parameters . N ; l red / that we decline in six values for the first parameter N D h 10 I h 4 I h 2 I hI 3h 2 I 2h and ten values for the second parameter l red D hI 2hI :::I 10h. We establish the convergence towards the analytical solution by considering the relative error with respect to the analytical solution

relative error D R C N Y N Y h 2 dV R C N Y 2 dV ( 102
)
Convergence curves are depicted in Figure 6. For l red D 3h, the smallest relative errors are obtained for h 4 6 N 6 h 2 ; see Figures 6(i) and 7. The relative error rapidly increases when N decreases as we loose the SUPG stabilization term in Problem 5. The relative error slowly increases when N increases, so greater values seem acceptable in practice. On the other hand, for N D h 2 , the smallest relative error is obtained when l red D h; see Figures 6(i) and 8. The relative error increases when l red increases. It increases even faster as the mesh is coarser, this is probably due to the small number of modes; see Figure 8. In practice, we will not choose a value l red D h, but we keep in mind that the value should not be too big and even more on a coarse mesh. The obtained non-local energy release rate is depicted in Figure 9 and compared with the lagrangian-based computation [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF]. We observe on the damage band tip that the discontinuity of the non-local energy release rate is well described by the computation based on modes. Iso-values are always orthogonal to the level set function isovalues. Oscillations are localized around the discontinuity of the solution, this is inherent to the SUPG but it is acceptable because it is localized in a band of a few elements thickness.

Application model

The constitutive model assumed in following numerical applications is a non-linear elastic damageable model with asymmetric traction/compression behavior and a delayed damage growth. It is summarized with two potentials. The first is the free-energy potential; see [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] '

.";d/ D .1 d/ 1 2 htr"i 2 C C h"i C W h"i C C 1 2 htr"i 2 C h"i W h"i (103)
where and are homogeneous Lamé parameters. The second is the dissipation potential

N Y D Y c c * N Y Y c 1 + C 1 a 1 exp a * N Y Y c 1 + C !!! (104)
where Y c is the critical energy release rate over which damage initiates and propagates, c is the characteristic time and a is an homogeneous parameter that influence the rate at which the non-local damage rate tends toward the upper bound value 1 c . Operators h i C and h i , respectively, stand for the positive and negative part of a quantity:

(1) for a scalar

hxi C D 1 2 .x C jxj/ (105) hxi D 1 2 .x jxj/ (106)
(2) for a second-order tensor

h"i C D 3 X iD1 h" i i C v i ˝vi (107) h"i D 3 X iD1 h" i i v i ˝vi (108)
where " i are eigenvalues and v i are eigenvectors of the tensor. It reads

D .1 d / htr"i C I C 2 h"i C C . htr"i I C 2 h"i / (109) Y D 1 2 htr"i 2 C C h"i C W h"i C (110)
N P d D 1 c 1 exp a * N Y Y c 1 + C !! (111)
where ( 111) is a delay-damage evolution model [START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF] derived from a positive, continuous, and convex dissipation potential that vanishes for a vanishing non-local energy release rate. The positiveness of the non-local damage rate implies the positiveness of the level set function rate. Moreover, this model is chosen because it upper bounds the non-local damage rate:

N P d< 1 c (112)
By considering [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], we write an upper bound on the level set function rate:

P < 1 c 1 N d 0 (113)
For a well chosen damage profile function, such as d./ D =l c and d./ D 2.=l c / .= l c / 2 for 0<<l c , see Appendix B for details, one can prove that

P < l c c (114)
The right hand side of (114) is taken as

l c c D c R (115)
where c R is the Rayleigh wave speed and is the crack velocity upper bound in brittle fracture. Then, it is ensured that

P <c R (116)
With (115), the critical time step for the forward Euler method is t c D h c R , noticing that c R <c l and using (52) the estimate of the critical time step is t c D h c l .

Kalthoff and Winkler experiments

In this section, we are simulating the Kalthoff and Winkler experiments [START_REF] Kalthoff | Failure mode transition at high rates of shear loading[END_REF]. A large amount of work has already been carried out on the simulation of such experiments, we summarize some of them: in [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF], a damage model is used until strain localization occurs, then a cohesive crack model is introduced; in [START_REF] Réthoré | An energy-conserving scheme for dynamic crack growth using the eXtended finite element method[END_REF][START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method (X-FEM)[END_REF][START_REF] Elguedj | An explicit dynamics extended finite element method. Part 1: Mass lumping for arbitrary enrichment functions[END_REF], brittle materials are considered, modeling is not focused on strain localization aspects but on efficient and accurate strategy to use the X-FEM in explicit dynamics. It focuses on schemes, enrichment, and lumping strategy;

in [START_REF] Menouillard | Time dependent crack tip enrichment for dynamic crack propagation[END_REF], authors focuse on a time-dependent enrichment to accurately treat crack propagation within the X-FEM; in [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], a phase field approach is used to regularize the crack into a smooth field and to model the loss of material properties around the crack, it leads to a model similar to continuum damage mechanics with an additional length scale. A projectile is launched at velocity v 0 on an edge of a rectangular l 2l steel plate as depicted in Figure 10. The plate is notched on this edge with two notches of length l=2 and the projectile impacts in between them. The plate is modeled under the plane strain assumption. The problem is symmetric, and therefore only half of the plate is modeled and we impose the boundary condition u e 2 D 0.The projectile is assumed to have the same impedance as the plate, and therefore the impacted part of the edge as velocity P u e 1 D v 0 2 . The material is a maraging steel modeled by a constitutive model given by potentials ( 103 The two meshes are referred respectively as coarse mesh and fine mesh in the following. The time step is chosen as t D 0:8h=c l . We focus our interest on the damage distribution, in Figures 11 and 14, on the front velocity, in Figure 12 and the dissipated energy E d , in Figure 13. The front velocity is estimated by computing the time derivative of a damage band length estimate:

L D 1 2l c Z C \D i dVi D 1; 2 (117)
where D 1 and D 2 are two domains illustrated in Figure 10. This way, we can estimate the front velocity for the two cracks. The dissipated energy is obtained by computing

E d D Z t f 0 Z C N Y N P d dV dt (118)
it is a good quantity to evaluate mesh sensitivity. We observe that the damage distribution is very close to the one obtained in [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF]. Two cracks appear and propagate into the domain. The presence of the second crack, although not reported in Kalthoff and Winkler experiments, exists in continuum damage mechanics-based approaches as in [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF]. In these studies, a damage area is initiated, but material instabilities are not lost. Note that in [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF], the maximum damage obtained in this area (close to d D 1) has a similar shape as our second crack. The approach presented in this paper is in fact limited to brittle failure, diffuse micro-cracking are not modeled within this paper. We are going further on this subject in the discussion section. There are almost no pathological mesh dependencies, neither on the dissipated energy, the front velocity, and the damage distribution. It exists in a mesh dependency on the size of the fully damaged domain c visible in Figure 14. The width of the domain is about h. This is due to the lack of enrichment to disconnect both sides of domain c . In [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF], a ramp Heaviside enrichment is introduced to remove this mesh dependency. The use of this enrichment in dynamics is not straightforward and remains a perspective in the approach. The angle between the main crack and the horizontal axis is about 65 ı (evaluated from the initiation location to a point where a bifurcation is observed, otherwise from initiation location to crack tip at instant t D 85:11s, the angle is about 71 ı ). Kalthoff and Winkler reported in their experimental studies an angle of about 70 ı . The average front velocity is below the Rayleigh wave speed. In Figure 22, we observe a first impulse in all curves, going away from plots bounding boxes. These impulses correspond to the initiation of damage where we place a circular iso-0 of radius l i . We observe as well small negative values, probably due to level set function reinitialization step.

Single-edge notched tension test

In this section, we are interested by the branching capacity of the TLS approach and the influence of the parameter a on the crack velocity. Experimental crack branching with a complex crack pattern can be found in [START_REF] Ravi-Chandar | An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching[END_REF]. Numerically, the branching problem of a single-edge notched tension test has already been treated by many authors [START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF][START_REF] Linder | Finite elements with embedded branching[END_REF]. We consider the domain as a L l rectangle with a notched edge. The length of the notch is e. It is modeled by a predefined damaged band C of thickness 2l c . The problem is illustrated in Figure 15. We apply a constant stress of magnitude T on both side of edges parallel to the notch. The problem is modeled under the plane strain assumption using the constitutive model given by potentials ( 103 The geometry and material parameters are chosen to be close to the dynamic crack branching test considered in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF] where we estimate our parameters product Y c l c equal to their parameter G c following the relationship with fracture mechanics given in [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] for which h is taken as 0. We focus our interest on the branching capacity by looking at the damage distribution in Figure 17. The influence of parameter a is seen on the maximum front velocity curve given in Figure 16 as well as in Figure 17 where the crack branching is delayed.

We observe in Figure 17 that the predefined damaged band opens until the fully damaged area has a one element thickness. The main crack propagates and branches into two cracks. The symmetry and the smoothness of the crack pattern seems better for smaller values of a. Crack branches sooner for bigger values of a even if we observe a longer main crack length. No other initiation occurs. In Figure 16, we observe that the maximum level set function rate is significantly slower than the upper bound c R . Parameter a affects the average velocity of the crack without significantly changing variations of curves. Oscillations amplitude is reduced for smaller a so that a does affect the stability of the propagation. The final increase in maximum front velocity corresponds to the emergence of the through-crack. When the damage front cut the domain boundary, a specific treatment is needed for the computation of modes and the level set function reinitialization into a signed distance function. This is not treated in this article and remain a perspective of the approach. The obtained results remain close to results obtained in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF] for a D 5. Compared with [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], the kinematic field discontinuity is modeled. Compared with [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF], no localization criterion is used to model transition from a continuous formulation to a discontinuous formulation.

DISCUSSION

Physical features

The foundation of the TLS approach is continuum damage mechanics. Nevertheless, the TLS approach tries to add into the formulation the physical benefit of modeling discontinuities that is the key idea of fracture mechanics. This choice is different from the phase-field approach that is built on fracture mechanics grounds but regularize the discontinuity for a computational benefit. This capacity of TLS to introduce kinematic fields discontinuities takes all its importance for problems with large sliding of cracks lips inducing mesh distortions.

The TLS approach, thanks to level set methods, has a great potential in the branching and merging of cracks. An interesting physical feature of the TLS approach concerns the branching capacity obtained without any numerical help and any additional physical criteria.

The parameters c and a introduced by the delay-damage evolution model serve two purposes: c imposes an upper bound to the crack velocity and a controls the average velocity. Parameter c is easy to determine from the maximum crack velocity. The choice of parameter a could then be made on the basis of an experimental average crack velocity.

Computational efficiency

Extra computational costs to handle non-locality is restricted to small spatial and temporal zones compared to other non-local approaches where the extra computational costs are widespread in space and time. Indeed, in the TLS approach, damage initiation is local leading to purely local computation before localization emerges somewhere. Afterwards, computational costs are increased, but in a smaller proportion than other approaches because the non-local computation is only needed in small areas that suffer damage. Computational costs increase with damage expansion.

The proposed extension to dynamics has almost no linear system to solve, it only appears when computing matrices OEP and OEP N with DOFs from C h and N h . Although all the benefits have not been implemented yet, level set methods used in the TLS approach are efficient to track the damage front. Computation of level set function rate is limited to C [ N . In this paper, the level set function is defined on all the domain, updated on C [ N and reinitialized on using a tree-based implementation of the distance function for the reinitialization step. This narrow band technique is well known in the community of level set methods; see [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF].

Limitations and perspectives

The approach presented in this paper is limited to brittle fracture under the small strain assumption as in previous paper dedicated to quasi static [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF][START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF]. It is extended to dynamic but leaves out the ramp Heaviside enrichment introduced in [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: efficient and accurate discretization for quasi-static loadings[END_REF]. Therefore, the fully broken area has a one-element thickness and is still mesh dependent. A perspective of this work is to determine a lumping strategy for the mass matrix as well as an efficient implementation with growing DOFs due to enrichment.

As first commented in the Kalthoff and Winkler experiments section, the presented approach is limited to brittle fracture in the sense that diffuse damage is not accounted for. A broader presentation is given in [START_REF] Moës | Coupling local and non-local damage evolution with the Thick Level Set model[END_REF] where the key idea is to relax the constraint imposed on the level set function kr k D 1 to kr k 6 1 to get a spatial damage distribution that satisfies at all time kr d k 6 f.d/ on . Without going into details, it means that the TLS approach can be extended to treat localized as well as diffuse fracturing.

So far, the level set reinitialization is based on a straightforward distance computation based on Equation [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] or a slightly improved one as used in Kalthoff and Winkler experiments. This computation does not correspond to the level set function needed in the TLS approach in presence of holes, notches. In Kalthoff and Winkler experiments, this low-quality reinitialization step affects non-local fields computation and level set function rate. A better suited method to reinitialize the level set function is the fast marching method. We believe that this method would improve simulation of the Kalthoff and Winkler experiments.

CONCLUSION

In this paper, we presented the TLS approach in the dynamical context of a material body subjected to impact loadings. We recalled that the TLS handles both the transition from a purely local computation to a non-local computation when damage initiates and the transition from a partially damaged area to a fully damaged area when damage propagates. We saw that the fully damaged area is not taken into account anymore so that it naturally brings the discontinuity in the displacement field (up to now, for a fully damaged band broader than one element).

The approach introduces a characteristic length in the formulation with the help of an auxiliary field: the level set function. To this end, we saw that the level set function must be a signed distance function. We wrote an evolution model for the level set function based on a non-local computation by introducing a thermodynamically conjugate pair of non-local variables: the non-local energy release rate and the non-local damage rate. We highlighted properties that must be verified by such fields so that the level set function remains a signed distance function during its propagation.

We gave a purely explicit time discretization based on the central difference time integration scheme for balance of momentum equations and the forward Euler time integration scheme for the level set function update.

We insisted on the fact that efficiency is preserved because computation of non-locality is limited to the damaged area.

We gave a new way to compute non-locality that is based on specific basis functions that weakly impose desired properties on non-local fields and level set function rate. We proposed an efficient computer implementation for which the construction of those basis function is not required.

Further improvements would be to introduce enrichment in the formulation so that discontinuity is treated even when the thickness of the fully damaged area is finer than one element. This work has already been treated in a quasi static context and is under development in a dynamics context. This would lead to a proper transition from damage to fracture.

APPENDIX A: NON-LOCAL DAMAGE RATE COMPUTATION

We derive the non-local damage rate expression considering [START_REF] Sethian | Fast methods for the Eikonal and related Hamilton-Jacobi equations on unstructured meshes[END_REF]:

Figure 1 .

 1 Figure 1. Domains and interfaces used in the thick level set approach.

55 )Figure 3 .

 553 Figure 3. Sets involved in modes definition.

  l red or more to corresponding mass matrix entry, 9:

Figure 4 .

 4 Figure 4. Damage band and narrow band after discretization.

Figure 5 . 8 <: 2x C y p 5 x 2 C and y< 1 2 r 4x 2 C y 2 4xy 16x C 8y C 16 5 x 2 C

 582222 Figure 5. V-shaped damage band, to illustrate the computation of the non-local energy release rate field.

Figure 6 .

 6 Figure 6. Relative error versus mesh size for (i) different diffusion coefficient N and l red D 3h; (ii) different mode support length l red and N D h=2.

Figure 7 .

 7 Figure 7. Relative error versus diffusion coefficient for (i) different mesh size h and l red D 3h; (ii) different mode support length l red and h D 3:906 10 3 m.

Figure 8 .

 8 Figure 8. Relative error versus mode support length for (i) different mesh size h and N D h=2; (ii) different diffusion coefficient N and h D 3:906 10 3 m.

Figure 9 .

 9 Figure 9. Non-local energy release rate computed with the mode-based method (left) and the lagrangianbased method (right) from [17]. Parameters are h D 7:813 10 3 m, N D h=2,a n dl red D 3h. Errors concentrate on the discontinuity of the solution. The mode-based is able to represent the discontinuity of the solution.

Figure 10 .

 10 Figure 10. Geometry for the Kalthoff and Winkler plate.

Figure 11 .

 11 Figure 11. Damage fields for Kalthoff and Winkler experiments.

  ) and (104). Material properties are E D 190 GPa, D 0:3, D 8000 kg/m 3 , Y c D 4:434 MPa, l c D 5 mm, c D l c =c R , a D 4 and d./ D 2.=l c / .= l c / 2 for 0<<l c . It leads to longitudinal wave speed c l D 5654 m/s, transversal wave speed c t D 3022 m/s, and Rayleigh wave speed c R D 2799 m/s. The length of the plate is l D 0:1 m. The velocity of the projectile is v 0 D 33 m/s. Numerical parameters are l i D 2 mm, N D h=2 and l red D 2h. Numerical computations are achieved on two 128 128 and 256 256 quadtree adaptive meshes made of three-noded triangular and alternated elements, the minimal level of refinement is 3 and the maximum level of refinement is, respectively, 7 and 8. It leads to h D 7:8 10 4 mandh D 3:9 10 4 m. The mesh adaptation parameters are l m D 5 mm and ˛D 1=5. The mesh is initially refined around the notch on a band of l m thickness.

Figure

  Figure 12. Front velocity versus time.

  Figure 12. Front velocity versus time.

Figure 13 .

 13 Figure 13. Dissipated energy versus time.

Figure 14 .

 14 Figure 14. Damage fields at instant t D 71:85s for Kalthoff and Winkler experiments and corresponding meshes.

  ) and (104). Material properties are E D 32 GPa, D 0:2, D 2450 kg/m 3 , Y c D 1:2 kPa, l c D 2:5 mm, c D l c =c R , a D 0:5I 1I 2I 5 (declined in four values) and d./ D 2.=l c / .= l c / 2 for 0<<l c .

Figure 15 .

 15 Figure 15. Geometry for the single-edge notched tension test.

Figure 16 .

 16 Figure 16. Maximum front velocity versus time for different values of parameter a.

Figure 17 .

 17 Figure 17. Damage fields for single-edge notched tension test at instant t D 78:75 µs.

  The domain C is called the damage band. All domains and interfaces are illustrated on Figure1. As said before, when initial damage is zero, we have D . Then, when damage initiates, we have D [ C ,t h a ti s , c D; . Finally, when damage propagates enough, we have D [ C [ c so that it exists a fully damaged area where d D 1. A tricky question is what should we do in this fully damaged area, keeping in mind that it does not mean that the stress tensor is zero, for instance, by considering asymmetric traction/compression constitutive model. In this work, we decide to force the constitutive model to provide .x;t/ D 0; x 2 c .t /.Thisway,weare

	C .t / D¹x 2 j 0 6 .x;t/ 6 l c º	(18)
	.t / D¹x 2 j .x;t/ < 0º	(19)
	c .t / D¹x 2 j .x;t/ > l c º	(20)
	c .t / D¹x 2 j .x;t/ D l c º	(21)
	on the level set function rate can be written P apply the change of coordinate to Equation (14):	.s; t/.W e D P
	P d.´;t/ D d 0 ..´; t // P .s; t/	(17)
	which means that the damage evolution model is linked to the level set function evolution that is
	always normal to the damage front. The evolution model for the level set function must provide such
	a behavior. We define the following domains and interface:

  is the length of the developed damage band along level set function gradients; see Figure2. It is the minimum between l c and value ´for which a discontinuous level set function gradient is found.We have R

	0	.s;t/ l.s;t/ Y.´;t/d 0 ..´;t// J.´;t/dŔ 0 d 0 ..´; t //J.´;t/d´( 23)
	where l.s; t/ N P d.s; t/ D	R l.s;t/ 0 0 l.s;t/ P d.´;t/J.´;t/dŔ	J.´;t/d´(	24)

C .t / dV D R .t/ R l.s;t/ 0 J .´;t/ d´ds. We therefore deduce from (22); see Appendix A for details:

  We now write Problem 5 using matrix notation. In order to solve the n red boundary value problems efficiently, a global matrix is assembled only once and Dirichlet boundary conditions are assembled in a matrix. It leads to several linear systems that can be written in a compact form: for ˛2 S red ,find¹D ˛º D b j 2 R n d and ¹E ˛º D c j˛ 2 R n such that

			9
	< : y.x/ D	˛2S red X	; = y ˛Q ' ˛.x/; y ˛2 R	(67)
	4.2.3. Efficient computer implementation.

we obtain [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] by identification of (A.1) and (A.2).

APPENDIX B: LEVEL SET FUNCTION RATE UPPER BOUND

For the two following damage profile functions (defined for 0<<l c ), the upper bound (113) reads

where f.l; / D

For these two damage profile functions, we obtain (114).

ACKNOWLEDGEMENT

We gratefully acknowledge the Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) Le Ripault. We would like to thank professor Antonio Huerta from the Polytechnic University of Catalonia for useful remarks on the SUPG and on the TLS in general.


Maria Sayu Yamamoto
email: maria-sayu.yamamoto@universite-paris-saclay.frsciencesconf.org

Fabien Lotte

Florian Yger

Sylvain Chevallier

User-specific frequency band and time segment selection with high class distinctiveness for Riemannian BCIs

Keywords: Riemannian geometry, oscillatory activity based BCIs, frequency band and time segment selection

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Abstract

User-specific settings are known to enhance Brain-Computer Interface (BCI) performances [1]. In particular, the optimal frequency band and time segment for oscillatory activity classification are highly user-dependent and many selection methods have been developed in the past two decades. However, it has not been studied well whether those conventional methods can provide optimal settings for the Riemannian BCIs, a recent family of BCI systems that utilize different data representation, based on covariance matrices, compared to conventional BCI pipelines. In this work, we proposed a novel frequency band and time segment selection method considering class distinctiveness on the Riemannian manifold. The class distinctiveness of each combination of frequency band and time segment is quantified based on inter-class distance and intra-class variance on the Riemannian manifold [2]. An advantage of this method is the user-specific settings can be adjusted without computationally heavy optimization steps. To the best of our knowledge, this is the first optimization method for selecting both the frequency band and the time segment on the Riemannian manifold. We evaluated the contributions of the 3 different selection models (frequency band, time or frequency band+time), comparing classification accuracy with a baseline (a fixed frequency band of 8-30 Hz and a fixed 2s time window) and a conventional popular method for non-Riemannian BCIs [1], on the BCI competition IV dataset 2a (2-class motor imagery). Our method showed higher average accuracy than baseline and a conventional method in all three models, and especially the frequency band selection model showed the highest performance. This preliminary result suggests the importance of developing new selection algorithms considering the properties of the manifold, rather than directly applying methods developed prior to the rise of Riemannian BCIs as they are. Blankertz B & al., IEEE Sig Proc Mag, 2007;[2] Lotte F & Jeunet., JNE, 2018

1

  User-specific frequency band and time segment selection with high class distinctiveness for Riemannian BCIs Maria Sayu Yamamoto * †1 , Fabien Lotte 2 , Florian Yger 3 , and Sylvain Chevallier 4 Laboratoire d'Ingénierie des Systèmes de Versailles (LISV, Université Paris-Saclay) -Laboratoire d'Ingénierie des Systèmes de Versailles (LISV), University paris saclay -10-12 Avenue de l'Europe, France 2 Inria Bordeaux Sud-Ouest / LaBRI (Inria Bordeaux Sud-Ouest / LaBRI (Univ. Bordeaux / CNRS / Bordeaux INP)) -Inria Bordeaux Sud-Ouest, LaBRI -200 Av. de la Vieille Tour, Talence, France 3 LAMSADE, PSL-Université Paris-Dauphine (LAMSADE) -LAMSADE, Université Paris-Dauphine, PSL -Place du Maréchal de Lattre de Tassigny, Paris, France 4 Laboratoire d'Ingénierie des Systèmes de Versailles, Université de Versailles Saint-Quentin-en-Yvelines, Université Paris-Saclay (LISV, UVSQ, Université Paris-Saclay) -Laboratoire d'Ingénierie des Systèmes de Versailles (LISV), Université de Versailles Saint-Quentin-en-Yvelines (UVSQ), University Paris-Saclay -10-12 avenue de l'Europe, France