R Desmorat
email: desmorat@lmt.ens-cachan.fr

F Gatuingt

M Jirásek

Nonlocal models with damage-dependent interactions motivated by internal time

Keywords: nonlocal continuum, localization, regularization, anisotropic damage, eikonal equation, metric

There are open questions concerning isotropic/anisotropic nonlocal models, as the possible evolution of internal length, as special treatments near boundaries and cracks. First, one makes nonlocal weight kernel function of information/wave propagation time normalized by an internal time and obtains strong localization with a non-spreading damage zone. The limit value of full damage is attained at a single point in 1D. Second, WKB approximation for wave propagation in 3D damaged media define interaction distances as solutions of an eikonal equation. This leads to the interpretation that damage, possibly anisotropic, curves the space in which the interaction distances are calculated.

Nomenclature

Notation Meaning x = max(x, 0) positive part of scalar x T positive part of symmetric second-order tensor T x norm of vector x

Introduction

The effect and the formulation of boundary conditions -such as free edges, notches and initial cracks -remain an open question for nonlocal models. The main drawback of the standard nonlocal (NL) integral theory [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] consists in the nonphysical interaction, through the nonlocal averaging process, of points across a crack or a hole. The definition of natural boundary conditions of vanishing normal strain gradient at a free edge is still under discussion for implicit gradient formulations [START_REF] Aifantis | The physics of plastic deformation[END_REF][START_REF] Peerlings | Gradient-enhanced damage model for quasi-brittle materials[END_REF]. The continuous nucleation of a crack of zero thickness is not so simple as the thicknesses of the localization band and also of the band in which damage tends to 1 are not zero. They are proportional to the internal length introduced in the nonlocal approach. Local behavior along free edges -i.e. with a vanishing internal length -has been obtained by some authors [START_REF] Pijaudier-Cabot | Comments on boundary effects in non local damage based models[END_REF][START_REF] Krayani | Boundary effect on weight function in nonlocal damage model[END_REF][START_REF] Pijaudier-Cabot | Non local damage model: boundary and evolving boundary effects[END_REF][START_REF] Bažant | Nonlocal boundary layer (nbl) model: overcoming boundary condition problems in strength statistics and fracture analysis of quasibrittle materials[END_REF]. The consideration of an internal length evolving with either the damage, the strain or the stress [START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF][START_REF] Pijaudier-Cabot | Nonlocal Damage Model with Evolving Internal Length[END_REF][START_REF] Simone | From continuous to discontinuous failure in a gradient enhanced continuum damage model[END_REF][START_REF] Nguyen | A damage model with evolving nonlocal interactions[END_REF][START_REF] Giry | Stress-based nonlocal damage model[END_REF][START_REF] Saroukhani | A simplified implementation of a gradient-enhanced model with transient lenght scale effects[END_REF] seems a way to properly bridge Continuum Damage Mechanics and Fracture Mechanics as the internal length may then vanish for large values of damage. Such a feature is nevertheless quite difficult to enforce in nonlocal integral theories, for instance by keeping the nonlocal connectivity matrix symmetric. See [START_REF] Comi | Numerical aspects of nonlocal damage analyses[END_REF][START_REF] Borino | A symmetric nonlocal damage theory[END_REF] for the study of the symmetry of the nonlocal weight.

One attempts here to propose a solution -bringing also questions -to these main difficulties, first within the framework of the integral theory. The idea [START_REF] Desmorat | Introduction of an internal time in nonlocal integral theories[END_REF][START_REF] Desmorat | Introduction of an internal time in nonlocal integral theories[END_REF] is to keep the nonlocal integral averaging process but to quantify the distance between points as an effective distance, i.e. as a distance function for instance of the geometry and the matter encountered between interacting points. One proposes to define such an effective distance with respect to a dynamic process: how information (wave) propagates between interacting points. This is made through the introduction of an internal time τ c , constant, instead of an internal length l c , measured as evolving. Note that in a similar manner local thermal expansions can be used instead of dynamic impulses [START_REF] Rojas-Solano | Interaction-based non-local damage model for failure in quasi-brittle materials[END_REF]. Waves propagations give us information on interacting / non interacting points for the definition of the nonlocal averaging. Such a point of view is consistent with the fact that dynamics is important to define a link between a characteristic time and a characteristic length, either when viscosity is introduced [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Allix | Delay-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF] or when the physical defects obscuration phenomenon encountered in high speed dynamics and multi-fragmentation is taken into account [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF]. Some authors even introduce the standard nonlocal theory by comparing the characteristic wavelength of the deformation field to an intrinsic length of the material [START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF][START_REF] Jirásek | Comparison of integral-type nonlocal plasticity models for strain-softening materials[END_REF], still a dynamics vocabulary.

The question of the numerical implementation with damage/softening models naturally arises and is addressed in 1D in section 4, in which the properties of the modified (wave) nonlocal approach are pointed out. The one-dimensional bar test studied has shown its usefulness to explore the localization characteristics of damage models with evolving internal length [START_REF] Pijaudier-Cabot | Nonlocal Damage Model with Evolving Internal Length[END_REF][START_REF] Nguyen | A damage model with evolving nonlocal interactions[END_REF][START_REF] Pham | From the onset of damage to rupture: construction of responses with damage localization for a general class of gradient damage models[END_REF][START_REF] Saroukhani | A simplified implementation of a gradient-enhanced model with transient lenght scale effects[END_REF]]. An important result of the 1D analysis performed with the nonlocal integral theory with internal time will be that mesh independency can be obtained with damage approaching 1 at a single point (and not in a spreading zone). The uniaxial pure geometric point of view that interacting points are "pushed apart" by the damage field present in between is finally extended to 3D by means of the differential geometry interpretation that damage curves the space.

The pupose here is to propose nonlocal appproaches that can apply to both isotropic and anisotropic damage.

Nonlocal integral theory with internal time

Softening constitutive equations classically lead to spurious dissipation modes and to mesh dependency. The need of the definition -and the introduction -of an internal length in the models is now established. But in which form? In a gradient form? In an integral form? From an internal viscosity (delay-damage) combined then with dynamics? The main idea of such regularizations is to average the variable -i.e. the thermodynamics force, denoted next V in the general case or Y or ε for damage -responsible for the strain localization. The procedure to define a nonlocal variable V nl from its local expression V introduces a characteristic length l c considered as a material parameter.

Nonlocal integral theories -Boundary effect

The standard nonlocal theory [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] uses the integral

V nl (x) = 1 V r (x) V α 0 x -ξ l c V(ξ)dξ = V α x -ξ l c V(ξ)dξ (1)
V r (x)

= V α 0 x -ξ l c dξ (2)
which is refered to as the standard averaging over the whole domain V in order to define a nonlocal quantity V nl -built from the local variable V -to be used in thermodynamics state or evolution laws. In Eq. (1) α 0 is the nonlocal weight function, positive and decreasing with respect to the distance xξ between the considered point x and all the structure points ξ. The normalizing factor is denoted V r (x) and l c is an internal or characteristic length. Two frequently used weight functions are the Gaussian function

α 0 (ζ) = e -1 2 ζ 2 (3)
and the bell-shaped polynomial function [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: Survey of progress[END_REF] α

0 (ζ) = 1 -ζ 2 2 (4)
where the brackets . . . denote the positive part, defined as x = max(x, 0). For the bell-shaped polynomial function it is common to call the characteristic length the interaction radius (because points at distance larger than l c do not interact at all) and to denote it as R instead of l c . The expressions (1) then define the same contribution to V nl of points x and ξ -across a crack as for any points x and ξ separated by the same distance ℓ = xξ (for instance path x 1 -x 6 versus path x 2 -x 5 of figure 1); this flaw has been pointed out and studied e.g. by Jirásek & Marfia [START_REF] Jirásek | Nonlocal damage models: displacement-based formulations[END_REF] and Pijaudier-Cabot et al [START_REF] Pijaudier-Cabot | Comments on boundary effects in non local damage based models[END_REF]. A numerical averaging adaptation neglecting the communication between some integration points can be found in [START_REF] De Vree | Comparison of nonlocal approaches in Continuum Damage Mechanics[END_REF];

-across a damaged zone: for example the interaction between points x 3 -x 4 is not affected by the zone at a damage level D 0 , i.e. at a Young's E(1 -D 0 ) much lower than the Young's modulus of the virgin material; much less studies address this point.

This illustrates the need to deal with such a boudary condition effect but also to make equivalent a real crack and a completely damaged zone.

Time instead of distance

In order to solve the problem, one proposes in this work to keep the nonlocal integral framework and to consider the information time propagation τ xξ between points x and ξ instead of the classical distance ℓ = xξ .

A nonlocal theory with internal time is then simply defined by replacing Eq. ( 1) by:

V nl (x) = 1 V r (x) V α 0 τ xξ τ c V(ξ)dξ = V α τ xξ τ c V(ξ)dξ (5) V r (x) = V α 0 τ xξ τ c dξ (6)
with τ xξ the information propagation time [START_REF] Desmorat | Non standard thermodynamics framework for robust computations with induced anisotropic damage[END_REF] taken next as the time for an elastic wave to propagate from point x to point ξ and τ c a material parameter. Such a nonlocal averaging process may apply to physical laws of different nature. Concerning continuum mechanics and elasticity, plasticity and/or damage, the variables V are often equivalent strains, accumulated plastic strains or strain energy densities.

Example: isotropic and anisotropic nonlocal damage models

For example, local isotropic damage theories for quasi-brittle materials define -the damage as a loss of stiffness D = 1 -Ẽ/E with Ẽ (resp. E) being the damaged (resp. initial) Young modulus, equation also rewritten in 3D as the elasticity law coupled with isotropic damage

σ σ σ = (1 -D)E E E : ǫ ǫ ǫ (7)
with σ σ σ, ǫ ǫ ǫ and E E E respectively the stress, the strain and the Hooke tensors;

-the damage evolution as a function of a local variable V either equal to the thermodynamic force Y = 1 2 ǫ ǫ ǫ : E E E : ǫ ǫ ǫ associated with damage (Marigo model) or to an equivalent strain as Mazars strain ε = √ ǫ ǫ ǫ : ǫ ǫ ǫ ,

D = g(max t V) (8)
with g(max t V) a nonlinear non-decreasing function and . . . the positive part of a tensor (obtained by applying the scalar positive part operator on the principal values). The maximum max t at current time t enforces a non-negative damage rate for any loading, Ḋ ≥ 0.

The nonlocal damage law is simply written as

D = g(max t V nl ) (9)
instead of Eq. ( 8) with the nonlocal averaging process [START_REF] Belytschko | Crack Propagation by Element-Free Galerkin Methods[END_REF]. An example using Mazars equivalent strain is the damage evolution law

D = g(κ) = 1 - ǫ 0 κ exp - κ -ǫ 0 ǫ f -ǫ 0 κ = max t εnl (10)
with parameters ǫ 0 (damage threshold) and ǫ f . Setting f = εnlκ, κ = g -1 (D), the Kuhn-Tucker loading/unloading conditions f ≤ 0, κ ≥ 0, κ f = 0 apply. For concrete, the microcracks due to tension are mainly orthogonal to the loading direction, when the microcracks due to compression are mainly parallel to the loading direction. The damage state has then to be represented by a tensorial variable D either a fourth rank tensor or a second rank tensor [START_REF] Murakami | A constitutive equation of creep damage in polycrystalline metals[END_REF][START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF][START_REF] Murakami | Mechanical modeling of material damage[END_REF][START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF][START_REF] Lemaitre | Mécanique des matériaux solides[END_REF]. The use of a second-order damage tensor is more convenient for practical applications (as well as for the material parameters identification) and this is the choice which is made. The damage anisotropy induced by either tension or compression is simply modeled by the consideration of damage evolution laws ensuring a damage rate proportional to the positive part of the strain tensor, i.e. a damage governed by the principal extensions [START_REF] Mazars | The unilateral behavior of damage concrete[END_REF].

The full set of constitutive equations for the local anisotropic damage model reads [START_REF] Desmorat | Modèle d'endommagement anisotrope avec forte dissymétrie traction/compression, 5e journées du Regroupement Francophone pour la Recherche et la Formation sur le Béton (RF2B)[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF][START_REF] Souid | Pseudodynamic testing and nonlinear substructuring of damaging structures under earthquake loading[END_REF] -Elasticity,

ǫ ǫ ǫ = 1 + ν E σ σ σ - ν E tr σ σ σ 1 1 1 (11)
with E the Young modulus and ν the Poisson ratio.

-Effective stress,

σ σ σ = (1 1 1 -D) -1/2 σ σ σ D (1 1 1 -D) -1/2 D + 1 3        tr σ σ σ 1 -1 3 η tr D --tr σ σ σ        1 1 1 (12)
where (.) D = (.) -1 3 tr(.) 1 1 1 stands for the deviatoric part and with η ≥ 1 the hydrostatic sensitivity parameter [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF].

-Damage criterion,

f ( ε, D) = ε -κ(tr D) (13)
with ε local Mazars equivalent strain and where

κ(tr D) = a • tan tr D aA + arctan ǫ 0 a ( 14
)
-Induced damage anisotropy governed by the positive extensions,

D = λ ǫ ǫ ǫ 2 (15)
In the rate independent formulation, the damage multiplier λ is determined from Kuhn-Tucker loading-unloading conditions f ≤ 0, λ ≥ 0, λ f = 0.

In total, 5 material parameters are needed: E, ν for elasticity, ǫ 0 as damage threshold and A and a as damage evolution parameters.

The model is simply made nonlocal, either from the standard integral theory or from the new integral nonlocal with internal time theory, by replacing Mazars equivalent strain ε by its nonlocal form in the damage criterion function, becoming

f ( εnl , D) = εnl -κ(tr D) (16)
instead of Eq. ( 13) with the nonlocal averaging process [START_REF] Belytschko | Crack Propagation by Element-Free Galerkin Methods[END_REF],

εnl (x) = V α τ xξ τ c ε(ξ)dξ (17)

Effective or "dynamic" distance -Vanishing internal length

In a plain and uncracked medium the internal length and internal time concepts are equivalent as

ℓ ≡ x -ξ = c 0 τ xξ ( 18
)
with c 0 the information celerity taken as the wave speed in the undamaged medium (a constant). Then, if the internal time is related to the internal length as l c = c 0 τ c the weight functions are equal,

α 0 x -ξ l c = α 0 τ xξ τ c (19)
In a homogeneously damaged medium at D = D 0 , the wave speed is proportional to the square root of the damaged Young's modulus and depends on the damage level as c0 = c 0 1 -D 0 [START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF] One has in this case

τ xξ = x -ξ c0 (21)
and

τ xξ τ c = x -ξ c0 τ c = x -ξ l c √ 1 -D 0 = l l c = x -ξ lc > x -ξ l c ( 22
)
which shows that the effective or "dynamic" distance l

= x -ξ √ 1 -D 0 ( 23
)
between two points increases, as expected, with damage. Eq. ( 22) defines in an equivalent manner an effective evolving internal length lc = l c √ 1 -D 0 (with then the definition for the distances unchanged) which tends to zero when D 0 tends to unity in accordance with Pijaudier-Cabot et al results [START_REF] Pijaudier-Cabot | Comments on boundary effects in non local damage based models[END_REF] of a material behavior becoming local on free edges (at least in the direction normal to the edge).

In 1D and in the non-uniform case, the effective distance is evaluated by integrating over the interval between x and ξ as

l = max(x,ξ) min(x,ξ) d l = max(x,ξ) min(x,ξ) dx ′ √ 1 -D(x ′ ) ( 24
)
and in 3D as l = c 0 τ xξ .

Equivalence between a crack and a damaged zone

In order to illustrate the formulation ability, consider the specimens of figure 2: (a) a square plate, (b) a notched specimen and (c) a specimen with a damaged zone at D = D 0 . Vertical z-acceleration fields are also plotted at the same instant for the 3 specimens. It is clearly noticed that the information initiated as an impulse at x = x 1 takes longer to reach point x 6 for both notched and damaged specimens (the wave generated in the example has to turn around the notch). Note the waves reflexion on upper notch and on upper damaged zone sides. When the damage D 0 becomes large (0.99 in the figure), notched and damaged specimens are found equivalent.

x p -x q l c x 1 x 2 x 3 x 4 x 5 x 6 x 1 0 l l c
Table 1: Connectivity table built from the classical distances xξ = x px q ( √ 2 and √ 5 replaced by their numerical values 1.414 and 2.236 to make easier further comparisons).

Connectivity matrices

To quantify the approach, consider six points x p in these specimens. The distances between these points, used in the standard nonlocal theory, are calculated in Table 1 for plate (a), with l denoting a fixed distance l = 0.5 m.

The standard connectivity table is the same (and is symmetric) for the 3 structures (a), (b) and (c). It does not take into account the presence of notches nor the occurence of damage as component 16

= component 25 = component 34 = l/l c .
In order to determine the connectivity table with the internal time formulation, one proceeds as follows:

-an impulse force according to the out-of-plane z-axis (close to a Dirac) is applied at each point x = x p in a dynamic finite element analysis of a 3D model of a plate (thickness of 0.1 m, 720000 nodes and 360000 elements, free boundary conditions),

-the times τ xξ at which the first peak of the z-acceleration information arrives at point ξ = x q are recorded. Divided by τ c they are put in the form of connectivity tables, -the expressions for τ xξ /τ c are synthesized by use of the relation l c = c 0 τ c and compared to Table 1 results for the different specimens.

For the concrete square plate (E = 35000 MPa, ν = 0.2, ρ = 2400 kg/m 3 ) and notched specimen one ends up to connectivity tables 2 and 3 where the notch presence is taken into account (16 and 61 components), leading for the notched specimen to τ (b) 16

/τ c = 1.494 l/l c instead of τ (a)
16 /τ c = l/l c for the plate with no notch. For the specimen with the damaged zone at D = D 0 = 0.99, the connectivity table corresponding to the nonlocal internal time analysis is Table 4. Again only the components 16 and 61 of the connectivity matrices are changed and found close to the value 1.494 l/l c obtained with reals notches: the damaged zone behaves as a notch, damage and notch being both taken into account by the proposed nonlocal with internal time analysis.

τ xξ τ c = c 0 τ xξ l c x 1 x 2 x 3 x 4 x 5 x 6 x 1 0 l l c 2.
τ xξ τ c = c 0 τ xξ l c x 1 x 2 x 3 x 4 x 5 x 6 x 1 0 l l c 2.008 l l c 2.233 l l c 1.420 l l c 1.494 l l c x 2 l l c 0 l l c 1.420 l l c l l c 1.420 l l c x 3 2.008 l l c l l c 0 l l c 1.420 l l c 2.233 l l c x 4 2.233 l l c 1.420 l l c l l c 0 l l c 2.008 l l c x 5 1.420 l l c l l c 1.420 l l c l l c 0 l l c x 6 1.494 l l c 1.420 l l c 2.233 l l c 2.008 l l c l l c 0
τ xξ τ c = c 0 τ xξ l c x 1 x 2 x 3 x 4 x 5 x 6 x 1 0 l l c 2.008 l l c 2.233 l l c 1.420 l l c 1.500 l l c x 2 l l c 0 l l c 1.420 l l c l l c 1.420 l l c x 3 2.008 l l c l l c 0 l l c 1.420 l l c 2.233 l l c x 4 2.233 l l c 1.420 l l c l l c 0 l l c 2.008 l l c x 5 1.420 l l c l l c 1.420 l l c l l c 0 l l c x 6 1.500 l l c 1.420 l l c 2.233 l l c 2.008 l l c l l c 0 Table 4:
Connectivity table built from times τ xξ for the specimen (c) with a damaged zone at D = D 0 = 0.99.

Straight 1D wave propagation

For a better understanding, let us focus on times τ 16 , τ 25 , τ 34 of wave propagation in specimen (c) from x 1 to x 6 , from x 2 to x 5 , from x 3 to x 4 . The corresponding distances are equal in the 3 cases (equal to l). Due to the presence of the notch the time τ 16 is of course larger than the time τ 25 corresponding to a straight path with no notch (and τ 25 /τ c = l/l c ).

A wave propagating along the path from x 1 to x 6 in specimen (c) slows down to speed c0 when meeting the damaged zone. If one only considers the straight path from x 1 to x 6 for the wave propagation, the propagation time is given by

τ 16 = l c 0 + e c 0 • 1 - √ 1 -D 0 √ 1 -D 0 ( 25
)
or

τ 16 τ c = l l c 1 + e l • 1 - √ 1 -D 0 √ 1 -D 0 ( 26
)
where e is the thickness of the damaged zone at D = D 0 (see Fig. 2). Fig. 3 shows the normalized increase of time τ 34 = τ (c) 34 with respect to the damage value D 0 and for different ratios e/l. The ratio (τ 34 /τ c )/(l/l c ) in case (c) is equal (resp. close) to unity for a zero (resp. small) damage, the nonlocal theory with internal time recovering then the standard nonlocal theory. The very large increase obtained for large values of the damage, enhanced by a large damaged zone thickness, shows that both a real crack and a highly damaged zone are equivalent in the internal time nonlocal framework.

Weight functions and nonlocal strain for piecewise constant strain and damage

As an example emphasizing how a damaged zone is taken into account, consider a bar in tension at a stress level σ. The length of the bar is 2L, a damaged zone at large D = D 0 corresponds to x ∈ [-l, l]. Damage is assumed piecewise constant in the bar, equal either to 0 or to D 0 . The nonlocal analyses evaluate the integrals (1) over the whole bar, 2 ) is considered. The normalizing factors V r (x) are compared for both analyses in figure 4a where D 0 = 0.99, l = 2.5 cm, L = 50 cm and where the characteristic length is taken as l c = 10 cm (twice the size 2l of the damaged zone). The standard normalization does not "see" the damaged zone and averages across it when the new nonlocal with internal time approach behaves for the undamaged domains almost as for two independent bars, as expected. The nonlocal strains obtained with both approaches are compared in figure 4b. The local strain is piecewise constant, the standard nonlocal strain εnl 0 = ε nl 0 corresponds to the dot line, the modified (wave) nonlocal strain εnl = ε nl new to the solid thick line. In this particular piecewise constant strain field, the formulation with internal time gives as expected a nonlocal strain field closer to the local strain field than the standard nonlocal approach for which too much importance is given to points over the damaged zone when performing the integral [START_REF] Borderie | Phénomènes Unilatéraux dans un Matériau Endommageable : Modélisation et Application à l'Analyse de Structures en Béton[END_REF]. Note that in case of structural failure with damage growth, strain localization leads to non homogeneous fields, as studied next.

ǫ nl (x) = 1 V r L -L α(x, ξ) ε(ξ)dξ V r (x) = L -L α 0 dξ (27) with strain ε(x) = σ/E if |x| > l and ε(x) = σ/E(1 -D 0 )

Strain-damage localization in a bar

A bar of a constant cross section subjected to uniaxial tension is simulated with a one-dimensional nonlocal damage model based on the stress-strain law σ = E(1 -D)ε and damage evolution law [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF], with equivalent strain equal to the positive part of strain, ε = ε .

Standard nonlocal approach

The simulation is first done using the standard nonlocal averaging approach. Due to boundary effects, damage would tend to localize at one of the boundaries. Since we are interested in the evolution of a process zone that forms inside the body, the solution localized at one boundary is artificially suppressed by enforcing symmetry. The finite element mesh consisting of 327 linear elements is refined near the axis of symmetry, so that the shape of the localized strain profile can be resolved with sufficient accuracy. The total length of the bar, 2L = 0.1 m, is five times the interaction radius R = 0.02 m used for evaluation of dimensionless variable ζ = |x-ξ|/R in nonlocal weight function [START_REF] Bažant | Nonlocal boundary layer (nbl) model: overcoming boundary condition problems in strength statistics and fracture analysis of quasibrittle materials[END_REF]. Parameters of the constitutive law are set to E = 30 GPa, ε 0 = 10 -4 and ε f = 10 -3 . This means that the tensile strength is f t = Eε 0 = 3 MPa and, if the strain remained uniform, the initial post-peak slope of the load-displacement curve would be ten times smaller than the slope of the elastic branch. Due to localization, the actual softening branch becomes steeper, but it does not exhibit snapback. Localization is triggered by reducing the sectional area of one element in the middle of the bar by 1‰ (the length of this element is 20 µm, i.e., R/1000).

The load-displacement diagram is plotted by the dash-dotted curve in Fig. 5a, in terms of the normalized stress (or load), σ/ f t , against the normalized displacement, u/u peak = u/2Lε 0 . Fig. 6a-c shows the evolution of the (local) strain, ε, damage-driving variable (maximum nonlocal strain), κ = max t εnl , and damage variable, D. Strain-like variables are normalized by the elastic limit strain, ε 0 , and the spatial variable x on the horizontal axis is normalized by the interaction radius R. It is obvious that the strain profile is becoming more narrow during the localization process, but the damage growth takes place in a relatively wide band. This is documented in more detail in Fig. 6d, which shows the process zone size L d normalized by the interaction radius R as a function of the maximum damage at the center of the zone, D max . Physically, the active part of the damage process zone should be understood as the set of points at which the damage is growing (at the current stage of loading). The corresponding curve is plotted by the solid line. It is seen that as D max approaches its limit value 1, the size of the active damage zone approaches 2R. For comparison, Fig. 6d also contains a dashed curve that corresponds to the size of the zone of growing local strain, which is smaller than the active damage zone (characterized by growing nonlocal strain).

Modified (wave) approach

As described in section 1 the strength of nonlocal interaction between two points now depends not just on their physical distance ℓ but rather on the time needed by an elastic wave to propagate from one point to the other. Since the wave speed is proportional to √ 1 -D, the same effect is obtained if the actual physical distance with differential dx is replaced by an effective distance with differential d l = dx/ √ 1 -D. A pure geometric approach can therefore be developed in the one-dimensional setting, with no need of dynamics computations, as the effective distance between points x and ξ is simply obtained as

l(x, ξ) = max(x,ξ) min(x,ξ) dx ′ √ 1 -D(x ′ ) (28)
and standard nonlocal averaging is modified as follows, according to the introduction of the internal time in nonlocal theory:

α(x, ξ) = α 0 l(x,ξ) R V r (x) V r (x) = L -L α 0 l(x, ξ) R dξ (29)
For uniform damage D, this would be equivalent to a reduction of the interaction radius R by the factor √ 1 -D. As D tends to 1, the modified interaction radius tends to zero, and the model response becomes local. In the general case of nonuniform damage, the reduction of the interaction radius is strong in the highly damaged region. Such a modification is expected to eliminate certain pathological effects, such as the growth of damage to its limit value 1 in a wide band of thickness 2R.

The intuitive expectations are confirmed by numerical simulations. Fig. 7 shows the evolution of the strain profile. Shrinking of the zone in which strain keeps growing is much more pronounced than for the standard nonlocal approach; see also the dashed curve in Fig. 8. Fig. 7a shows the early stages of localization, during which the width of the strain profile is comparable to the (initial) interaction radius R. Later on, the zone of growing strain becomes much thinner; see parts (b) and (c) of the figure. Note that the scales on the horizontal axis in different parts of the figure are different, and the width of the profile in the bottom part is two orders of magnitude smaller than R. The zone in which the nonlocal strain is growing is also progressively shrinking; see Fig. 9. The evolution of the damage profile is depicted in Fig. 10. Part (a) of that figure provides the global picture, while parts (b) and (c) are close-ups focusing on the late stages of the softening process. It is apparent that damage tends to 1 only at the center of the process zone, which is confirmed by the solid curve in Fig. 8.

As the process zone gets thinner, the structural response becomes more brittle. Consequently, the load-displacement diagram exhibits a snapback; see the solid curve in Fig. 5a. Convergence upon mesh refinement from element size h = 8 mm to 0.5 mm is documented in Fig. 5b. Two sensitivity analyses are performed, first with respect to damage parameter ε f , which is varied from 10 -3 to 4 × 10 -3 while keeping the interaction radius R = 20 mm, second with respect to R, varied from 20 mm to 60 mm while keeping ε f = 10 -3 . The corresponding graphs show that the dissipated energy can be controlled by variation of ε f and R but the snapback is always present, in contrast to the dash-dotted curve with a long tail, corresponding to the standard nonlocal approach. This indicates that the reduction of effective characteristic length at late stages of damage evolution is closely linked to the increased brittleness. It seems next to impossible to obtain simultaneously a properly localized damage profile and a long tail of the load-displacement curve. Similar phenomena have been observed in [START_REF] Jirásek | Localization study of a regularized energetic damage model[END_REF] for a family of energetic gradient damage models. Note nevertheless that (i) the shape of quasi-static post-peak response of concrete obtained from direct tension experiments is still under discussion [START_REF] Mazars | The unilateral behavior of damage concrete[END_REF][START_REF] Nouailletas | Comportement d'une discontinuité dans un géomatériau sous sollicitations chemo-mécanique-expérimentations et modélisations[END_REF], and (ii) in the present and many other approaches the permanent strains due to damage are not modeled. Such a modeling may prove pertinent, it is nevertheless left to further work (see for instance references [START_REF] Borderie | Phénomènes Unilatéraux dans un Matériau Endommageable : Modélisation et Application à l'Analyse de Structures en Béton[END_REF][START_REF] Ragueneau | Damage model for concrete-like materials coupling cracking and friction, contribution towards structural damping : first uniaxial applications[END_REF][START_REF] Lebon | Analyse de l'endommagement des structures de génie civil: techniques de sous-structuration hybrides couplées à un modèle d'endommagement anisotrope[END_REF][START_REF] Richard | Continuum damage mechanics based model for quasi brittle materials subjected to cyclic loadings: Formulation, numerical implementation and applications[END_REF] for permanent strain representation in the pure Continuum Damage Mechanics framework).

To illustrate the evolution of the nonlocal weight function during the damage process, Fig. 11 shows α(x, ξ) as a function of ξ for several selected values of x at the initial undamaged state (part (a)) and at an early and late stage of damage evolution (parts (b) and (c)). The initial and modified weight functions are compared in Fig. 12a. Finally, Fig. 12b shows the evolution of the scaling factor V r (x), which is initially constant in the central part of the bar and variable only near the boundaries, but later decreases dramatically in the damaged zone. This indirectly reflects the reduction of the effective interaction radius.

For one-dimensional problems, the wave approach requires only a relatively simple modification of the standard nonlocal damage model implemented in a finite element code such as OOFEM [START_REF] Patzák | Design of object oriented finite element code[END_REF][START_REF] Patzák | OOFEM -an object-oriented simulation tool for advanced modeling of materials and structures[END_REF]. In a typical implementation, each Gauss integration point x p creates and stores a list of pointers to all the "neighboring" Gauss points x q whose distance is smaller than the interaction radius R. Simultaneously, a list of corresponding nonlocal interaction weights a pq is precomputed and stored. The interaction weights depend on the value of the nonlocal weight function but also on the weight attributed to point x q in the Gauss quadrature formula and on the Jacobian of the isoparametric mapping. It is a good practice to compute first the non-normalized weights

a * pq = w q J q α 0 x p -x q R ( 30
)
which would be used in a numerical approximation of the integral

V α 0 x p -ξ R V(ξ)dξ ≈ q α 0 x p -x q R V(x q )w q J q = a * pq V(x q ) ( 31
)
and then normalize them to

a pq = a * pq r a * pr ( 32
)
This procedure leads to interaction weights that exactly satisfy the condition q a pq = 1 for each p, and so the numerical approximation of nonlocal averaging by a finite sum does not modify a constant local field.

For the standard nonlocal model, the interaction weights a pq are precomputed for all interacting pairs of Gauss points before the start of the actual incremental-iterative nonlinear analysis. They are stored in the form of lists and reused in every evaluation of the nonlocal variable. For the modified wave approach, these lists are also precomputed and stored, but during the nonlinear analysis it is necessary to adjust the weights, taking into account the current damage field. In 1D, it is sufficient to sort the Gauss points q in the list of the nonlocal neighbors of Gauss point p in the ascending order according to their spatial coordinate x q and mark the position of the receiver point p (which always belongs among its own nonlocal neighbors). Then the list is processed first from point p to the end of the list ("to the right") and then from point p to the beginning of the list ("to the left"), and the effective distances are computed by accumulating the contributions of individual segments between the nearest neighbors. To describe that formally, we can assume for simplicity that the global numbering of Gauss points is done in the ascending order according to the spatial coordinate. The evaluation of effective distances lpq ≡ l(x p , x q ) can then be described by the recursive scheme lpp = 0 (33)

lp,q+1 = lpq + x q+1 -x q 2        1 1 -D q + 1 1 -D q+1        , q = p, p + 1, p + 2, . . . e p -1 (34)
lp,q-1 = lpq + x qx q-1 2

       1 1 -D q + 1 1 -D q-1        , q = p, p -1, p -2, . . . b p + 1 ( 35
)
where b p and e p are respectively the numbers of the Gauss points at the beginning and at the end of the list of neighbors of Gauss point p. In practice, the evaluation of distance in the direction to the right (or to the left) is terminated when lp,q+1 (or lp,q-1 ) exceeds the interaction radius R. Since 1/ √ 1 -D ≥ 1, the domain that contributes to the nonlocal average at point p is always contained in the interval (x p -R, x p + R) and the initial list of neighbors never needs to be extended. Once the effective distances are known, the interaction weights are updated as

ã * pq = w q J q α 0       lpq R       (36) ãpq = ã * pq r ã * pr ( 37
)
Of course, an extension of this procedure to multiple spatial dimensions is not straightforward, and it turns out to be very helpful to consider an alternative interpretation of the wave approach, to be presented next.

Interpretation in terms of differential geometry and extension to anisotropic damage

In the previous so-called modified (wave) approach, the points of a structures have been "pushed farther apart" by the damage field present between them. This has been formulated in 1D based on the concept of internal time for nonlocal theories, as d l = (1 -D(x)) -1/2 dx. In the uniaxial setting, there is in fact no need for dynamics computations, since a purely geometric equivalent formulation can be used. If one needs to model as many wave propagations as the number of Gauss points in a structure, the computational cost may become prohibitive in 3D! The wave type has also to be chosen.

Let us show that the 1D purely geometric approach can be extended to 3D if isotropic damage is considered, dropping the reference to the concept of internal time.

Eikonal approximation in the 3D case with isotropic damage

Let us consider the unstrained medium with a non-homogeneous isotropic damage field D(x), constant in time, previously used to calculate the information propagation time τ xξ by a dynamic simulation. Setting G = G(1 -D) = effective shear modulus and K = K(1 -D) = effective bulk modulus, the displacement u corresponding to an elastic wave satisfies the equation

G∇ 2 u + K + 1 3 G ∇ (∇ • u) + (∇ • u)∇ K - 2 3 G + ∇ G • ∇u + ∇u • ∇ G -ρ ∂ 2 u ∂t 2 = 0 ( 38
)
which is the wave equation for an isotropic but non-homogeneous medium, with non-homogeneity caused by the damage gradient ∇D. Of course, Eq. ( 38) in general does not have a closed-form solution.

If we are interested in how a wave with a high frequency ω and a small wavelength λ 0 propagates at velocity c 0 , i.e. at a high wave number k 0 = ω/c 0 = 2π/λ 0 , we can seek the complex-valued displacement field in the polar form of WKB approximation [START_REF] Wentzel | Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik[END_REF][START_REF] Kramers | Wellenmechanik und halbzählige Quantisierung[END_REF][START_REF] Brillouin | La mécanique ondulatoire de Schrödinger: une méthode générale de résolution par approximations successives[END_REF],

u(x, t) = u 0 (x) e i(k 0 S (x)-ωt) ( 39
)
where u 0 is the position-dependent wave amplitude and S (x) is the so-called eikonal function, selected such that the initial wave front is its zero level set. It will now be demonstrated that S (x) corresponds to the effective distance l(ξ, x) between the initial wave front at which S (ξ) = 0 and a given (arbitrary) point x. Consequently, the eikonal function can be interpreted as a level set function. Differentiation of the displacement field (39) leads to

∇u = (ik 0 u 0 ⊗ ∇S + ∇u 0 ) e i(k 0 S -ωt) ≈ ik 0 u 0 ⊗ ∇S e i(k 0 S -ωt) (40)
∇ • u = (ik 0 ∇S • u 0 + ∇ • u 0 ) e i(k 0 S -ωt) ≈ ik 0 ∇S • u 0 e i(k 0 S -ωt) (41) ∇∇u = -k 2 0 u 0 ⊗ ∇S ⊗ ∇S + ik 0 ∇u 0 ⊗ ∇S + ∇(ik 0 u 0 ⊗ ∇S + ∇u 0 ) e i(k 0 S -ωt) ≈ -k 2 0 u 0 ⊗ ∇S ⊗ ∇S e i(k 0 S -ωt) (42) ∇ (∇ • u) = -k 2 0 ∇S (∇S • u 0 ) + ik 0 ∇S (∇ • u 0 )∇ (ik 0 ∇S • u 0 + ∇ • u 0 ) e i(k 0 S -ωt) ≈ -k 2 0 ∇S (∇S • u 0 ) e i(k 0 S -ωt) (43)
∇ 2 u = ∇ • ∇u = -k 2 0 (∇S • ∇S )u 0 + ik 0 ∇S • ∇u 0 + ∇ • (ik 0 u 0 ⊗ ∇S + ∇u 0 ) e i(k 0 S -ωt) ≈ -k 2 0 (∇S • ∇S )u 0 e i(k 0 S -ωt) (44)
∂ 2 u ∂t 2 = -ω 2 u 0 e i(k 0 S -ωt) (45)
The approximation signs above mean "at high wave number k 0 ", i.e. at high frequency. WKB approximation also uses assumptions on the gradients of material properties, which are here proportional to the damage gradient, since ∇ K = -K∇D and ∇ G = -G∇D.

Assuming that the wave length is much smaller than the inverse value of the damage gradient, λ 0 ≪ 1/ ∇D , we can neglect all terms in the wave equation [START_REF] Nouailletas | Comportement d'une discontinuité dans un géomatériau sous sollicitations chemo-mécanique-expérimentations et modélisations[END_REF] which contain the product of the displacement gradient (or divergence) with the gradient of reduced (damaged) elastic moduli, because they are much smaller than the terms of order k 2 0 that are considered by the WKB approximation. After these simplifications, Eq. ( 38) is reduced to

-Gk 2 0 (∇S • ∇S )u 0 -K + 1 3 G k 2 0 ∇S (∇S • u 0 ) + ρω 2 u 0 = 0 (46)
This is a set of homogeneous equations for u 0 , which has a nontrivial solution if and only if

det K + 1 3 G k 2 0 ∇S ⊗ ∇S -ρω 2 -Gk 2 0 ∇S 2 1 1 1 = 0 ( 47
)
It is easy to show that ∇S ⊗ ∇S has eigenvalues ∇S 2 and 0, the latter being a double eigenvalue. Therefore, the condition of zero determinant ( 47) is satisfied if either

ρω 2 -Gk 2 0 ∇S 2 = K + 1 3 G k 2 0 ∇S 2 (48)
or ρω 2 -Gk 2 0 ∇S 2 = 0 (49) The first case corresponds to a P-wave (∇ × u = 0) and gives

k 2 0 ∇S 2 = ρω 2 K + 4 3 G (50)
while the second case corresponds to an S-wave (∇ • u = 0) and gives

k 2 0 ∇S 2 = ρω 2 G ( 51
)
Both cases can be presented in the form

k 2 0 ∇S 2 = ω 2 c 2 0 (1 -D) ( 52
)
where c 0 is the wave speed in the undamaged medium, given by c 0 = (K + 4G/3)/ρ for the P-wave and c 0 = G/ρ for the S-wave.

The ambiguity in the definition of k 0 and S (only their product matters) can be removed by setting k 0 = ω/c 0 . Eq. ( 52) then reduces to the (isotropic) eikonal equation

∇S = 1 √ 1 -D (53)
When the WKB approximation is combined with the condition that the initial wave front is the zero level set of function S , it turns out that S represents the signed distance l from the initial wave front ξ. For a uniformly damaged medium, the usual distance ℓ = xξ is scaled by the factor 1/ √ 1 -D; see Eq. ( 23) in 1D. This is extended to a three-dimensional non-uniformly damaged medium by the choice of high frequency elastic waves, therefore by defining effective distances S ≡ l from eikonal equation [START_REF] Sethian | Level set methods and Fast Marching Methods: evolving interfaces in computational geometry[END_REF]. The result derived here for isotropic damage is found independent of the type of wave, since the same eikonal equation applies to both Pwaves and S-waves. For the undamaged medium with D = 0, we obtain the condition ∇S = 1, which means that S ≡ l = ℓ = x-ξ .

Note that the eikonal equation is derived as ∇S = n in geometric optics (with n the index of the medium), and allows then for the determination of a light ray path without solving the wave equations (electrodynamic Maxwell's equations in optics). Efficient algorithms to solve the eikonal equation are available, they include, e.g., Dijkstra-like algorithms and Fast Marching Methods [START_REF] Dijkstra | A note on two problems in connexion with graphs[END_REF][START_REF] Tsitsiklis | Efficient Algorithms for Globally Optimal Trajectories[END_REF][START_REF] Sethian | Level set methods and Fast Marching Methods: evolving interfaces in computational geometry[END_REF].

Anisotropic eikonal equation in anisotropic damage case

In the case of induced anisotropic damage, represented by a symmetric second-order tensor D with components D i j , the coupling of initially isotropic elasticity tensor with anisotropic damage Ẽ E E = Ẽ E E(D) is model dependent (a fourth-order tensor with components D i jkl can similarly be considered). For the example given in section 1.3, the damaged stiffness tensor is given by

Ẽ E E = 2G (1 1 1 -D) 1/2 ⊗(1 1 1 -D) 1/2 - (1 1 1 -D) ⊗ (1 1 1 -D) 3 -tr D + K 1 - 1 3 η tr D H (tr ǫ) 1 1 1 ⊗ 1 1 1 ( 54
)
where H(x) is the Heaviside function, and tensorial product ⊗ is defined by the identity (A⊗B) : C = A • C • B for all symmetric tensors A, B, C. The elastic waves in the unstrained medium are solutions of

∇ • Ẽ E E : ∇u -ρ ∂ 2 u ∂t 2 = 0 (55)
The WKB approximation at high frequency leads to

Ẽ E E . . . ∇(∇u) -ρ ∂ 2 u ∂t 2 = 0 (56)
The same developments as in the previous isotropic damage case now lead to an anisotropic eikonal equation for effective distance S = l,

∇ l •       u 0 u 0 • Ẽ E E ρc 2 0 • u 0 u 0       • ∇ l = 1 (57)
The choice c 0 = (K + 4 3 G)/ρ where k 0 = ω/c 0 can be made in order to recover the limiting case of isotropic damage with D = D 1 1 1 and η = 1.

Both isotropic and anisotropic damage cases are unified within a single eikonal equation

∇ l • g -1 • ∇ l = 1 (58)
with g -1 denoting a symmetric second-order tensor, dependent on damage. In the anisotropic case g -1 also depends on the position, on the wave type (which is even evolving when damage anisotropy and intensity increase) and on the damage model itself.

An alternative is to define the information propagation times in a damaged medium from a straightforward extension of eikonal equation [START_REF] Sethian | Level set methods and Fast Marching Methods: evolving interfaces in computational geometry[END_REF] to anisotropy by setting g -1 = 1 1 1 -D and rewriting [START_REF] De Vree | Comparison of nonlocal approaches in Continuum Damage Mechanics[END_REF] as

∇ l • (1 1 1 -D) • ∇ l = 1 (59)

Space curved by damage

The Hamilton-Jacobi standard theory allows to show that the eikonal equation in fact defines a geodesic of differential geometry [START_REF] Smirnov | A Course on Higher Mathematics[END_REF]: the information time τ xξ defined from Eq. ( 53) or (58) multiplied by the wave velocity c 0 is the shortest path in a Riemannian space curved by damage.

In order to derive this result in the present damage mechanics analysis, let us consider the tensor g = (1 1 1 -D) -1 with covariant components g i j = g ji = [(1 1 1 -D) -1 ] i j as the (symmetric) metric of a differential geometry manifold, and let us define the increment of effective distance l by

d l2 = dx • g • dx = dx i g i j dx j (60)
In such a curved space contravariant and covariant coordinates differ. One has the usual rule to lower or raise indexes with the metric, a i = g ik a k , A i j = g ik g jl A kl , g ik g k j = δ i j . The case of isotropic damage is the case of spherical metric

g = 1 1 1/(1 -D), g i j = δ i j /(1 -D), with d l2 = dx i dx i /(1 -D) = dx 2 /(1 -D).
Note that other contributions where geometry has been used to describe damage can be found in [START_REF] Lagoudas | A gauge theory of damage[END_REF][START_REF] Lagoudas | Finite element implementation of the gauge theory of damage[END_REF][START_REF] Ganghoffer | A new framework in nonlocal mechanics[END_REF][START_REF] Ganghoffer | New concepts in nonlocal continuum mechanics and new materials obeying a generalised continuum behaviour[END_REF].

The effective distance l between points x 1 and x 2 is a measure of curvilinear abscissa along a geodesic, l =

x 2 x 1 ẋ • (1 1 1 -D(x)) -1 • ẋ dt = x 2 x 1 L(x, ẋ) dt (61)
If extremity positions x 1 and x 2 are fixed, the Euler-Lagrange extremality condition of path length l is

∂L ∂x - d dt ∂L ∂ẋ = 0, with Lagrangian L(x, ẋ) = ẋ • (1 1 1 -D(x)) -1 • ẋ.
This gives the standard geodesic equation as

d 2 x k d l2 + Γ k i j dx i d l dx j d l = 0 where Γ k i j = 1 2 (δ kl -D kl ) [(1 1 1 -D) -1 ] li, j + [(1 1 1 -D) -1 ] l j,i -[(1 1 1 -D) -1 ] i j,l (62)
The Levi-Civita connection is used here. The Christoffel symbols Γ k i j are damage-dependent, they vanish in a homogeneously damaged medium (in which geodesics are straight lines).

The Hamilton-Jacobi theory can provide additional information [START_REF] Smirnov | A Course on Higher Mathematics[END_REF], because the change δ l of the effective distance between points ξ and x resulting from a change of x by δx (at fixed ξ) can be expressed as

δ l = ∂L ∂ẋ • δx + x ξ ξ ξ ∂L ∂x - d dt ∂L ∂ẋ • δx dt (63)
Since the integration is done along a geodesic, the integral in Eq. ( 63) vanishes. The gradient of the effective distance can thus be expressed as

∇ l ≡ ∂ l ∂x = ∂L ∂ẋ = (1 1 1 -D) -1 • ẋ ẋ • (1 1 1 -D) -1 • ẋ ( 64
)
and it is easy to show that it satisfies the anisotropic eikonal equation

∇ l • (1 1 1 -D) • ∇ l = 1 (65)
In the case of isotropic damage, this reduces to ∇ l = (1 -D) -1/2 . The eikonal function (the wave front level set function l) stands for the effective distance, algebraic in the present variational calculations, defined as based on the propagation along a geodesic. When l is divided by the wave velocity c 0 in an undamaged medium, it can be interpreted in the isotropic damage case as the travel time τ ξx of a wave propagating from point ξ to point x. As the geodesic is then the path of the shortest travel time of a wave propagating from ξ to x, the minimal path principle can be seen as the shortest travel time principle. In addition, the eikonal equation, therefore the numerically efficient Fast Marching Methods, can be used to find the shortest wave path and hence the effective length l introduced into the nonlocal integral theory.

Due to non-homogeneity of the damage field, the space in which the effective distances lxξ = | l(x, ξ)| between points ξ and x are measured (and then used as an argument of the nonlocal interaction function α( l/l c )) is not a flat space anymore: it is curved by the damage field. The shortest distance between two points is not a straight line but a geodesic curve, either of the eikonal equation (65), or of the more classical geodesic equation (62).

Eikonal NonLocal (ENL) integral damage models use the effective distance lxξ , instead of the physical distance xξ , within the averaging weight function, setting then for the initial internal time (IT) approach

V nl (x) = 1 V r (x) V α 0       lxξ l c       V(ξ)dξ, V r (x) = V α 0       lxξ l c       dξ ENL-IT approach (66)
A potential alternative would be to use the mapping x → x, ξ → ξ, from metric g to perform the integration in the "effective" (curved) space, with V(ξ) = Ṽ( ξ), V r (x) = Ṽr (x) and

V nl (x) = 1 V r (x) V α 0       lxξ l c       Ṽ( ξ)d ξ, V r (x) = V α 0       lxξ l c       d ξ ENL-E approach (67)

1D and 3D diagonal cases

In one dimension, the geometric interpretation gives

d l = dx √ 1 -D (68)
In three dimensions, if the coordinate axes are aligned with the principal axes of damage, the damage tensor is represented by a diagonal matrix

D =         D 1 0 0 0 D 2 0 0 0 D 3         (69)
The effect of damage anisotropy (different values of D 1 , D 2 and D 3 ) is illustrated by the calculation of the effective distances lx , ly , lz -found to be different -from the origin ξ = (0, 0, 0) to a point x at the same distance ℓ but respectively taken as x = (ℓ, 0, 0), x = (0, ℓ, 0), x = (0, 0, ℓ):

lx = ℓ 0 dx √ 1 -D 1 ly = ℓ 0 dy √ 1 -D 2 lz = ℓ 0 dz √ 1 -D 3 (70)
In full 3D, the knowledge of the metric g, directly obtained from the calculation of the anisotropic damage, allows then to calculate lengths, i.e. effective distances, from any point, at any time step, in a continuum damage mechanics computation.

Remark on gradient-enhanced nonlocal models

Let us make a final remark concerning the gradient-enhanced formulation with the nonlocal equivalent strain εnl (see section 1.3) or any other quantity V nl determined from the differential equation [START_REF] Peerlings | Gradient-enhanced damage model for quasi-brittle materials[END_REF] εnl -

1 2 l 2 c ∇ 2 εnl = ε (71)
with l c representing a (constant) characteristic length.

A straightforward transformation of such a partial differential equation by the use of metric g gives εnl -

1 2 l 2 c 1 det g ∇ • det g g -1 • ∇ εnl = ε (72)
which defines the nonlocal equivalent strain in what can be called the Eikonal NonLocal gradient-enhanced damage models. Eq. ( 72) can be derived from the "effective" eikonal nonlocal integral approach ENL-E; for details see the Appendix.

For isotropic damage D(x, t), possibly evolving, g = 1 1 1/(1 -D), det g = 1/(1 -D) 3 , and Eq. (72) reads εnl -

1 2 l 2 c (1 -D) 3/2 ∇ • (1 -D) -1/2 ∇ εnl = εnl - 1 2 l 2 c (1 -D)∇ 2 εnl - 1 4 l 2 c ∇D • ∇ǫ nl = ε (73)
At unloading points with decreasing εnl , damage and its gradient remain constant and Eq. ( 73) can be solved for the nonlocal equivalent strain. In the active part of the process zone where damage is growing, one has f = εnlκ(D) = 0, D = κ -1 ( εnl ) (Eq. ( 10) for example) and ∇D = dκ -1 d εnl ∇ εnl . In such a loading case the nonlocal equivalent strain εnl is the solution of εnl -

1 2 l 2 c 1 -κ -1 ( εnl ) ∇ 2 εnl - 1 4 l 2 c dκ -1 d εnl ∇ǫ nl 2 = ε (74)
This indirectly introduces an evolving internal length, damage-dependent as in [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Pham | From the onset of damage to rupture: construction of responses with damage localization for a general class of gradient damage models[END_REF]. Let us point out that there is a division by 73) which makes the recovery of a local strain at D = 1 somewhat tricky. It is nevertheless automatic if function κ -1 becomes constant (equal to 1) with a continuously vanishing derivative with respect to εnl , so the last term on the left-hand side of (74) vanishes at D = 1 -. This is the case for damage law [START_REF] Denoual | A damage model for the dynamic fragmentation of brittle solids[END_REF].

√ 1 -D in term 1 2 l 2 c (1 -D) 3/2 ∇ • ((1 -D) -1/2 ∇ εnl ) of Eq. (
The extension to second-order anisotropic damage is straightforward from the geometric point of view with metric g = (1 1 1 -D) -1 . An interesting particular case is the tensile damage case (e.g. tension along direction 1), for which the damage tensor is diagonal with D 11 = D 1 as the only non-zero component, and then D = diag[D 1 , 0, 0] in basis (e x , e y , e z ) and there is no equivalent strain gradient along directions 2 and 3 (e y • ∇ εnl = e z • ∇ εnl = 0). Then det g = 1/(1 -D 1 ) and the governing equation can be written as εnl -

1 2 l 2 c 1 -D 1 ∂ ∂x 1 -D 1 ∂ εnl ∂x = ε (75)
In the plane of full damage (D 1 = 1), i.e., the crack plane, the nonlocal equivalent strain εnl tends to the local value ε.

Conclusion

An internal time, instead of an internal length, has been introduced into the integral nonlocal theory, leading to the definition of an effective or "dynamic" distance and to an internal length evolving with the damage field present between interacting points. The important feature is that the distance between these points is not the argument of the nonlocal averaging weight function anymore. It is replaced by the information on wave propagation time between them. Pre-computations of wave propagation in the considered structure are needed to determine the corresponding nonlocal connectivity matrix, with the open questions of their updating, of the wave type and of the computational cost. The standard nonlocal theory is recovered far from nonconvex boundaries and for low damage. The presence of cracks and notches is naturally taken into account by the wave propagation within the proposed nonlocal theory with internal time. Such a modified (wave) approach leads to an equivalence between a crack and a highly damaged zone, since points shielded by a notch have a reduced contribution to the nonlocal averaging. An important point is that no assumption on the medium isotropy is made. The proposed nonlocal framework includes anisotropy, either initial or induced.

Convergence upon mesh refinement has been demonstrated in a uniaxial example of a bar loaded in tension, with the full damage (D = 1) reached at a single point (and not in a spreading zone). At late stages of the localization process, when the damage concentrates in a band much smaller than the initial one, the global response in terms of the load-displacement diagram becomes very brittle. However, the total dissipation per unit sectional area remains finite and can be controlled by adjusting the nonlocal interaction radius or the parameter that controls the damage evolution.

The pure geometric point of view that the effective distance between interacting points is increased by the damage field in between has finally been extended to 3D by means of the differential geometry interpretation that damage, possibly anisotropic, curves the space. WKB approximation for high-frequency wave propagation in a damaged medium has allowed us to define effective (interaction) distances as the solution of an eikonal equation, either isotropic or anisotropic, depending on the nature, scalar or tensorial, of the damage variable. This builds a pure geometric approach with no need anymore of many additional dynamics computations. Such a result will simplify the numerical implementation in 2 or 3D, as efficient algorithms for solving the eikonal equation (e.g. Fast Marching Methods) are available for both the isotropic and the anisotropic cases.

The pure geometric approach applies to nonlocal integral as well as gradient-enhanced theories, making the internal lengths evolve and defining Eikonal NonLocal (ENL) damage models. A comparative study of different nonlocal averaging schemes and of the corresponding gradient enhancements is left to future work. as over the Euclidian unit sphere Ω e r ⊗ e r dΩ = 4π 3 1 1 1 and where ∇2 Ṽ = ∇i ∇i Ṽ = pi j ∇ j ∇i Ṽ, defined as the trace of the inverse of metric tensor times the Hessian, is the Laplacian (in curved space) of the local equivalent strain Ṽ at the considered point x(x). Following [START_REF] Peerlings | Gradient-enhanced damage model for quasi-brittle materials[END_REF], one has then for the Gaussian weight α 0 ( l/l c ) = e -1 2 ( l/l c ) 2 , detailing terms up to order 2 only,

Ṽnl = Ṽ + 1 2 l 2 c ∇2 Ṽ + . . . Ṽnl - 1 2 l 2 c ∇2 Ṽnl = Ṽ + . . .
At the second order we recover the implicit Peerlings gradient nonlocal enhancement Ṽnl -1 2 l 2 c ∇2 Ṽnl = Ṽ, coupled with damage as in the physical coordinate system

V nl = V nl (x) = Ṽnl (x), V = V(x) = Ṽ(x), as then nonlocal variable V nl is solution of V nl - 1 2 l 2 c 1 det g ∇ • det g g -1 • ∇V nl = V Specimen (a) Specimen (b) x 1 x 2 x 3 x 6 x 5 x 4 x 1 x 6 x 1 x 2 x 3 x 6 x 5 x 4 4
Specimen (c)

1 X 2 X 3 X 6 X 5 X 4 Specimen (a) X 1 X 2 X 3 X 6 X 5 X 4 Specimen (b) X 1 X 2 X 3 X 6 X 5 X 4
Specimen (c)

  otherwise, and with α(x, ξ) = α 0 |x-ξ| l c for the standard nonlocal analysis, α(x, ξ) = α 0 τ xξ τ c for the nonlocal analysis with internal time. The Gaussian weight function α 0 (ζ) = exp(-1 2 ζ

x 1 x 2 x 3 x 6 x 5 x 4 D=D 0 Figure 1 :

 401 Figure 1: Notched and damaged specimens.

X

Figure 2 :Figure 3 :Figure 4 :RFigure 5 :

 2345 Figure 2: Square plate, notched and damaged specimens (notch thickness e = 0.1m). Geometry and pictures of wave propagation at the same instant.

Figure 6 :

 6 Figure 6: Standard nonlocal approach: evolution of (a) normalized strain ε/ε 0 , (b) damage-driving variable κ = max t εnl , (c) damage variable D, and (d) normalized process zone size, L d /R.

Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 789101112 Figure 7: Evolution of normalized local strain, ε/ε 0 , for the modified (wave) approach

Table 2 :

 2

Connectivity table built from times τ xξ for the concrete square plate (a) (no notch).

Table 3 :

 3 Connectivity table built from times τ xξ for the notched specimen (b).

Acknowledgements

The third author gratefully acknowledges financial support of the Czech Science Foundation under project 108/11/1243. He would also like to thank Ecole Normale Supérieure de Cachan for hosting him as a visiting professor.

Appendix : Gradient-enhanced model based on "effective" nonlocal integral averaging (ENL-E)

Let us consider the "effective" eikonal nonlocal averaging scheme ENL-E (Eq. ( 67)) in an infinite medium, in which the integration is done over a curved space equipped by metric g.

Let us use polar geodesic coordinates in the space curved by damage, so that d l2 = d ξ • p • d ξ = dξ • g • dξ defines a polar metric p, det p = 1. One sets x = ξ(x) = 0 0 0 as the origin, so that the first component of polar position ξ is ξ1 = lxξ = l (point at distance l from the origin point x = 0 0 0 along a geodesic). The volume measures are d ξ = l2 d ld Ω = det g dξ with d Ω denoting the solid angle measure in the curved space. One has then

Taylor expansion of local equivalent strain in curved space gives Ṽnl

where integrals of scalar terms α 0 l l c l2+n with respect to l and integrals of spherical terms

l over the solid angle are independent in an infinite medium. Due to normalisation by V r (x) the term multiplying the local equivalent strain Ṽ = V is classically found equal to 1. The middle term, contaning the integral Ω ξ ξ ξ l d Ω vanishes due to spherical symmetry. The integrations over l -from 0 to ∞ -are performed along geodesics of the curved space. The integration with respect to the solid angle can be performed at any distance l, and so it is evaluated in the neighborhood of the origin, in Euclidian tangent frame at negligible curvature, with Riemann change of variable r =