Estelle Berthier

Vincent Démery

Laurent Ponson

Damage spreading in quasi-brittle disordered solids: I. Localization and failure

We propose a novel approach inspired from non-local damage continuum mechanics to describe damage evolution in heterogeneous quasi-brittle materials. Material heterogeneities are introduced at a mesoscopic continuous scale through spatial variations of the resistance to damage. The damage field evolution is computed from irreversible thermodynamics principles by assuming that the elastic energy released during loading is dissipated into failure. The onsets of damage localization and catastrophic failure in the material are studied as a fonction of the strength of the heterogeneities and the interaction function involved in the non-local formulation of the model. The predictions obtained numerically are explained theoretically for weak heterogeneities using a linear stability analysis and confirmed through a complementary approach based on a global energy minimization. Two distinct quasi-brittle failure behaviors are identified: For interaction functions that impose a reloading of the material points after the occurence of a damage event, the damage grows rather uniformly in the material until catastrophic failure takes place. On the contrary, when damage events trigger reloading, but also a sufficiently strong unloading in some material regions, catastrophic failure is preceded by a stable regime of damage localization characterized by a lengh scale emerging from the structure of the load redistribution. Our study reveals the cooperative nature of the damage localization process, showing that quasi-brittle failure emerges from the interaction between the elements constituting the material. It also highlights the central role played by the mechanism of load redistribution that is shown to control the failure behavior of quasi-brittle solids.

Macroscopic failure properties of quasi-brittle materials emerge from the complex evolution of a large number of microcracks in interaction [START_REF] Kachanov | Elastic solids with many cracks -a simple method of analysis[END_REF][START_REF] Lockner | Quasi-static fault growth and shear fracture energy in granite[END_REF][START_REF] Bazant | Nonlocal damage theory based on micromechanics of crack interactions[END_REF], making their prediction a challenging task. Yet, a remarkable property of these materials is that, despite this complexity and the variety of microscopic damage mechanisms involved like e.g. crack face friction [START_REF] Lawn | Nonlinear stress-strain curves for solids containing closed cracks with friction[END_REF][START_REF] Poon | Damage accumulation and hysteretic behavior of max phase materials[END_REF], crack kinking [START_REF] Horii | Compression-induced microcrack growth in brittle solids: axial splitting and shear fracture[END_REF][START_REF] Ravichandran | A micromechanical model for high strain behavior of ceramics[END_REF] or pore collapse [START_REF] Fortin | Acoustic emission and velocities associated with the formation of compaction bands in sandstone[END_REF][START_REF] Baro | Similarity between the compression of a porous material and earthquake[END_REF], their failure behavior shares many similarities. They generally display an initial regime where damage events are homogeneously distributed within the material, followed by a localization regime during which failure events concentrate into bands that, subsequently, lead to the catastrophic failure of the specimen. Another striking feature of the failure behavior of quasi-brittle solids is the intermitent evolution of the damage that grows through bursts of activity separated by silent periods. These giant fluctuations revealed from the acoustic signals emitted during fracture [START_REF] Lockner | The role of acoustic emission in the study of rock[END_REF][START_REF] Fortin | Acoustic emission and velocities associated with the formation of compaction bands in sandstone[END_REF] display characteristic power law statistics that are still poorly understood [START_REF] Hirata | Omori's power law aftershock sequences of microfracturing in rock fracture experiment[END_REF][START_REF] Petri | Experimental evidence for critical dynamics in microfracturing processes[END_REF][START_REF] Guarino | An experimental test of the critical behavior of fracture precursors[END_REF][START_REF] Alava | Statistical models of fracture[END_REF][START_REF] Davidsen | Scaling and universality in rock fracture[END_REF][START_REF] Girard | Failure as a critical phenomenon in a progressive damage model[END_REF]. Last but not least, quasi-brittle fracture leaves scale-free failure patterns that remain undeciphered [START_REF] Weiss | Fracture and fragmentation of ice: a fractal analysis of scale invariance[END_REF][START_REF] Baud | Compaction localization in porous sandstones: A spatial evolution of damage and acoustic emission activity[END_REF][START_REF] Carpinteri | Size effect on tensile fracture properties: a unified explanation based on disorder and fractality of concrete microstructure[END_REF][START_REF] Renard | Constant dimensionality of fault roughness from the scale of micro-fractures to the scale of continents[END_REF]. The observation of a failure behavior common to a large range of quasi-brittle materials motivates this work that aims at (i) identifying the physical ingredients underlying it and (ii) propose a unified theoretical description that captures it.

To describe the damage spreading preceding the failure of quasi-brittle materials, two distinct approaches have been proposed. From a continuum mechanics perspective, failure processes are described at a coarse scale and described through a continuous damage field. This field provides the amount of damage accumulated in a material point of the specimen that describes an elementary volume including a large number of microcracks. Microcracking affects the local mechanical behavior of the material through the degradation of its elastic properties, leading to a local softening of the material [START_REF] Castaneda | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF][START_REF] Heap | The evolution of elastic moduli with increasing crack damage during cyclic stressing of a basal from mr. etna volcano[END_REF]. Within this so-called damage mechanics theory, the evolution of the damage field can be inferred from a damage driving force, often defined as the rate of elastic energy released as damage increases.

The material damages when this quantity reaches a critical value that characterizes the material resistance. This approach has been shown to capture various aspects of the failure behavior of quasi-brittle solids, like damage localization [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Lasry | Localization limiters in transient problems[END_REF][START_REF] Lemaitre | A course on damage mechanics[END_REF][START_REF] Aifantis | On the role of gradients in the localization of deformation and fracture[END_REF][START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]. However, it relies on the introduction by hand of non-local effects through e.g. the dependance of the damage driving force to the damage level in an extended region of the material or a gradient of the damage field in the expression of the damage driving force [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in the uniaxial tests with gradient damage models[END_REF]. This amounts to introduce an internal length that controls damage spreading, the physical origin of which is however poorly understood. In addition, these approaches are limited to homogeneous materials so, among many aspects of quasi-brittle failure, the intermittent evolution of the damage field and the scale-free failure patterns cannot be captured. Here, we will not address the origin of non-local effects and the emergence of a length scale characterizing the damage growth process that will be treated in subsequent articles (Berthier and Ponson; Démery et al.). Instead, we will focus on the combined effect of the material heterogeneities and the non-local interactions that will be shown to play a centrale role on damage localization and failure.

To capture the effect of material heterogeneities on damage spreading, alternative approaches based on a discrete description of failure processes have been proposed. In these models, solids are described by an array of elements with randomly distributed properties, as e.g. their failure thresholds. An element is broken when the load applied to it reaches a critical value. The load is then redistributed in the network either to satisfy equilibrium conditions or according to a redistribution law imposed by the model [START_REF] Herrmann | Statistical Models for the Fracture of Disordered Media[END_REF][START_REF] Zapperi | Plasticity and avalanches behavior in microfracturing phenomena[END_REF][START_REF] Alava | Statistical models of fracture[END_REF]. One of the most successfull model is the fiber bundle model that allows for the exploration of the interplay between disorder and elastic interactions as the bundle is progressively loaded up to catastrophic failure [START_REF] Pradhan | Crossover behavior in failure avalanches[END_REF]. These models predict not only the macroscopic failure behavior of the network and how microstructural parameters may affect it [START_REF] Pradhan | Dynamic critical behavior of failure and plastic deformation in the random fiber bundle model[END_REF][START_REF] Manzato | Fracture strength of disordered media: Universality, interactions, and tail asymptotics[END_REF], but they also provide a detailed description of the damage spatio-temporal evolution and its fluctuations [START_REF] Delaplace | Damage cascade in a softening interface[END_REF][START_REF] Lennartz-Sassinek | Emergent pattens of localized damage as a precursor to catastrophic failure in a random fuse network[END_REF][START_REF] Manzato | Damage accumulation in quasibrittle fracture[END_REF]. Their main drawback however is that they are not derived from thermodynamical principles so they do not ensure energy conservation during damage growth. As a result, they only provides a qualitative description of quasi-brittle failure.

A major challenge is to conciliate these both antagonist approaches and describe the transfer of mechanical energy into fracture energy during the damage of solids as described in continuous damage mechanics, but taking into account the effect of material heterogeneities and microcracks interactions that are succesfully captured by discrete models. [START_REF] Delaplace | Progressive damage in discrete models and consequences on continuum modeling[END_REF] have used lattice models to provide insights on continuous damage theory. They showed that a length scale characterizing damage spreading emerges from the interaction between microcracks, paving the way for non-local models with an internal length that may evolve with time [START_REF] Pijaudier-Cabot | A review of non local continuum damage: Modeling of failure[END_REF]. [START_REF] Van Mier | Fracture mechanisms in particle composites: statistical aspects in latteice type analysis[END_REF] and [START_REF] Rinaldi | Rational damage model of 2d disordered brittle lattices under uniaxial loadings[END_REF] have taken lattice models to the next level, by applying them to simulate damage spreading within realistic microstructures. More recently, the new generation of discrete approaches in quasi-brittle failure, the so-called discrete element models, have opened many perspectives for the simulation of quasi-brittle failure [START_REF] Kun | Approach to failure in porous granular materials under compression[END_REF][START_REF] Scholtès | A dem model for soft and hard rocks: Role of grain interlocking on strength[END_REF]. Another successfull numerical approach is the one proposed by [START_REF] Tang | Numerical simulation of progressive rock failure and associated seismicity[END_REF] and [START_REF] Amitrano | From diffuse to localised damage through elastic interaction[END_REF], as it takes into account material microstructural disorder while describing damage processes at a continuous scale. Their approach relies on the finite element computation of the stress field after each damage event, allowing for an incremental prediction of the damage spreading process even beyond localization. These works have greatly improved our ability to describe accurately failure of quasi-brittle solids in numerical simulations.

However, a deep understanding of the basic mechanisms underlying damage spreading, localization and failure is still missing, as well as a unified theoretical framework built on them. We believe that it can be found at the interface between discrete and continuous models, taking advantage of the strength of both approaches. Here, we narrow the gap between these two points of view by proposing a novel model that rigorously derives from energy conservation, but takes into account material heterogeneities. In our approach, heterogeneities represent a mesoscopic volume that contains a large number of microcracks so that a damage variable that goes continuously from zero for an intact element to unity for a fully broken one can be used. The local stiffness of the material is affected by the damage field through a weighted integral made over the whole specimen. We will show that the introduction of such nonlocal behavior allows the description of the load redistributions mechanism taking place after failure events while ensuring an energetically consistent model formulation.

Finally, material evolution is inferred from energy conservation through the introduction of a thermodynamic driving force for damage that describes the transfer of elastic energy stored in the loaded specimen into fracture energy.

A similar approach was proposed by [START_REF] Lyakhovsky | Distributed damage, faulting, and friction[END_REF], even though the elastic energy released is instead absorbed into entropy generation. Our description is sufficiently general so that two important ingredients can be varied and hence studied: (i) The strength of the microstructural disorder that is introduced at the mesoscale through a statistical distribution of the damage energy; (ii) the load redistribution mechanism, controlled through the nonlocal interactions. In the following, we focus on non-local interactions with internal length, resulting in a finite-range redistribution process. This differs from the redistribution function underlying damage processes in 3D elasto-damageable that shows power law behavior (Démery et al.). In addition, the redistribution computed from 3D elasticity displays a quadrupolar symmetry that ressembles to the Eshelby solution for the mechanical fields around a soft inclusion embedded in an infinite elastic medium [START_REF] Eshelby | The determination of the elastic field in an ellipsoidal inclusion, and related problems[END_REF] while in our model, the sign of the redistribution varies with the distance to the event. Despite these differences, our approach captures qualitatively well the main features of the redistribution process, namely the decay of the reloading with the distance to the fracture event and the variations of the sign of the redistribution implying that some regions of the material are reloaded while other one are unloaded.

Practically speaking, finite-range redistribution regularizes the localization process that will be shown to be characterized by a single localization mode. As a result, localization and catastrophic failure can be predicted analytically from a linear stability analysis of the homogeneously damaged state. We show that our predictions capture well the material behavior computed numerically, even for a sigificant level of material disorder. This theoretical analysis is complemented by an independant procedure based on global energy minimization that refines the previous predictions and sheds light on the nature of these instabilities. Overall, our study reveals the central role of the damage driving force redistribution in the mechanical behavior of quasi-brittle materials and the importance of considering material heterogeneities to describe their evolution towards localization and failure. A description of quasi-brittle failure as a four steps process also emerges: An initial elastic regime is followed by a progressive increase of damage with a rather homogeneous spatial distribution. Beyond some localization threshold, the damage grows heterogeneously to accumulate in some regions of the material only. This localization process ultimately leads to catastrophic failure of the whole specimen at some critical loading. Interestingly, some redistribution functions lead to unstable failure without prior localization. Our extensive numerical study shows not only the good agreement with the analytical predictions but also that our approach captures most generic features observed during quasi-brittle fracture experiments, namely localization and catastrophic failure, but also intermittency of the microcracking process and scale-free fracture pattern -the characterization of the precursors to failure and the spatial structure of the damage field will be described in a companion paper (Berthier and Ponson). Our approach is here applied to one and two dimensional systems for some specific types of interaction. However, the methodology proposed is amenable to the exploration of more complex situations with other interaction functions. Hence, we believe that this work paves the way towards the development of a thermodynamically consistent description of quasi-brittle failure that captures the co-action of microcrack interactions and material heterogeneities.

The paper is organized as follows: Section 2 introduces the assumptions of the model, defines the loading conditions and the material properties, and describes how damage evolution is derived from it. The section 3 is devoted to the failure behavior of homogeneous materials while heterogeneous materials are investigated numerically in section 4. Our approach brings out two major steps as the material is driven towards failure, namely damage localisation and catastrophic failure. In section 5.1, these transitions are explored analytically through a linear perturbation analysis of the homogeneous material problem, and the predictions are used to interpret the simulation results. Finally, we propose in section 5.1 an alternative approach to the stability analysis based on a global energy minimization that sheds light on the nature of damage localization and failure in quasi-brittle heterogeneous materials.

An energy-based damage model for heterogeneous materials with tunable interactions

Our approach relies on the description of the material at a mesoscopic scale, intermediary between the microscale at which microfracturing processes take place and the macroscale at which the material is homogeneous. As a result, it allows for the introduction of material heterogeneities which will be shown to play a central role in the damage evolution, while describing the mechanical quantities like stress, strain but also damage through continuous fields.

The interactions within the materials are introduced through an interaction function that describes how the driving force for damage is redistributed after a damage event.

Evolution law of the damage field

In the following, we consider a material constituted of N D heterogeneities that are distributed equidistantly in space and indexed by their position x in the specimen Σ of dimension D and of size L D , where D equals 1 or 2, with periodic boundary conditions. The size of the heterogeneities is noted ξ, so that the specimen size is effectively N D ξ D . The material is clamped between two rigid plates, as schematically represented in Fig. 1. The bottom plate is maintained fixed whereas a uniform macroscopic displacement ∆, perpendicular to the x-axis in 1D and to the (x, y)-plane in 2D, is applied to the upper plate. Due to the clamping conditions used in our model, the microscopic displacement is homogeneous, equal to the macroscopic one, such that the macroscopic force writes as

F = Σ σ( x)d x = ∆ Σ k[d( x)]d x (1)
where σ is the stress field and k[d( x)] = σ( x)/∆ the local stiffness which depends on the level of damage of the heterogeneity. However, k does not depend directly on the local damage parameter d( x) that goes from zero when the element is intact to unity when it is fully broken. Instead, we introduce a dependency on the nonlocal parameter d defined as

d( x) = (α * d)( x) = Σ α( x -x )d( x )d x (2)
where α is a weight function that will be specified below, that is used in the expression of the stiffness

k[d( x)] = k d( x) . (3)
The introduction of a non-local variable allows to introduce interactions in the material, as originally proposed by Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]. However, contrary to this work and various studies inspired from it, the interactions are introduced in our model at the level of the damage field, and not on the driving force. In other words, the stiffness in x depends on the damage level in an extended region defined by the range of the weight function. This description was chosen to ensure a proper energy balance during the transfer of mechanical energy into fracture energy, as shown in Appendix A. What are the physical implications of this non-local approach? As shown subsequently in Sec. 5.1, using a non-local damage variable produces a redistribution of the damage driving force after failure that is fully controlled by the shape of the weight function. A typical example of redistribution function is shown schematically in Fig. 1(a). For elements close to the heterogeneity located at the center that has just damaged, the driving force will increase (red areas) while elements located further away might be unloaded (blue areas), depending on the actual shape and sign of the weight function. The introduction of such a non-local failure behavior results in a complex temporal and spatial organization of damage illustrated on the damage spatio-temporal evolution of Fig. 1(b) that resembles qualitatively to the experimental observations. It also leads to a non-monotonic macroscopic failure response that will be the central point of this study. To predict damage evolution, the total energy E tot The damage driving force that applies on the central element results from a spatial average of the driving force that applies on the other elements weighted by the interaction function α represented here. This interaction function also interprets as the spatial structure of the load redistribution that takes place in the material after the central heterogeneity is damaged. (b) Spatio-temporal evolution of the damage field for a 1D specimen of N = 2000 elements. Initially, damage is zero everywhere (d = 0 for a driving ∆ = 0, upper line), until catastrophic failure of the specimen for ∆ = ∆ c (d = 1, lower line).

of the system that comprises the elasto-damageable material and the loading device is expressed, under some imposed loading ∆, as the sum of the elastic energy E el , the fracture energy E d and the work of the external force W

E tot = E el + E d -W = Σ 1 2 ∆ 2 k d( x) d x + Σ d( x) 0 Y c ( x, d)d d d x - ∆ 0 F ( ∆)d ∆.
(4)

where Y c ( x, d) refers to the damage energy of the medium in x at a damage level d. It is defined so that an increase δd of damage in a region δ x of the material results in an energy dissipation δE d = Y c ( x, d) δd δ x.

Focusing on the quasi-static response of the material, we introduce the thermodynamic driving force for damage as

F( x) = - δE tot δd = Y ( x) -Y c ( x). (5)
Under quasi-static loading conditions, i.e. assuming that damage mechanisms take place at a much faster rate than the external driving, the damage increases in the material under fixed applied displacement ∆. This implies that the work of the external force does not contribute to the damage energy release rate Y which writes as

Y ( x) = - δE el δd = (α * Y )( x). (6)
Interestingly, with the introduction of the non-local damage parameter in the expression (4) of the elastic energy, the energy release rate Y takes a nonlocal form and writes as the convolution of the weight function α with the local energy release rate defined as

Y ( x) = - ∆ 2 2 dk dd . (7)
The mechanical equilibrium that derives from the condition δE tot δd = 0 during damage growth is thus reached when

Y ( x) = - δE el δd = δE d δd = Y c ( x). (8)
In Appendix A, we show that this damage evolution law can be derived equivalently from the variation of the total energy of the system with respect to the external driving ∆.

Material parameters

Fracture energy and material heterogeneity

To account for the heterogeneities of the material, we introduce a damage energy field dependent on the position x in the medium through a random quenched noise g c drawn from a constant distribution of zero average value and variance σ. In addition, the material resistance is chosen to depend on the damage level at a points x as

Y c ( x, d) = Y c0 [1 + g c ( x, d) + ηd( x)] (9)
where η > 0 is a hardening parameter. The increase of the average damage energy

Y c ( x) x = Y c0 + η d( x)
x with d is often assumed in damage models (see e.g. [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in the uniaxial tests with gradient damage models[END_REF]). It qualitatively reflects microscopic toughening mechanisms taking place during material failure like friction at crack faces [START_REF] Lawn | Nonlinear stress-strain curves for solids containing closed cracks with friction[END_REF][START_REF] Poon | Damage accumulation and hysteretic behavior of max phase materials[END_REF], wing cracks [START_REF] Ashby | The damage mechanics of brittle solids in compression[END_REF][START_REF] Ravichandran | A micromechanical model for high strain behavior of ceramics[END_REF] or crack bridging [START_REF] Bower | A 3-dimensional analysis of crack trapping and bridging by tough particles[END_REF].

It also allows to explore material behaviors with an extended post-peak response and a softening. Indeed, without hardening (η = 0), the load redistribution as described in Fig. 1(a) induces systematically an increase of the driving force of the damaging elements even at constant applied load, resulting in an unstable failure of this element, and subsequently of the whole material right after the linear elastic regime.1 This is indication that hardening is a relevant physical mechanism involved in quasi-brittle materials displaying an extended softening regime.

Material stiffness

To account for the degradation of the elastic properties of the material during damage, we describe the stiffness decay of a single element with the damage parameter through a polynomial law

k(d) = k 0 [ad 3 -(a + 1)d + 1]. ( 10
)
k 0 is the stiffness of the intact material, taken homogeneous in the material, and a is a constant verifying -1 < a < 1/2 to ensure that the stiffness is a decaying function of d. This polynomial law allows for the exploration of a large range of damage behaviors for a single elements, as illustrated in Fig. 2(a)

where the stiffness variations with d are shown for three values of a. Indeed, the sign of a gives the curvature of k(d), and so the stability of the individual elements since ∂Y ∂d ∆ = -∆2 2 d 2 k dd 2 = -3∆ 2 ad as derived from Eq. ( 7). For negative values of a and no hardening, a point that starts to damage experiences a growing driving force under fixed displacement condition, resulting in an abrupt failure. On the contrary, positive values of a lead to a stable damage process, so that the external displacement has to be increased to damage further the fiber.

The case a = 0 corresponds to the largely used damage law k(d) = k 0 (1d) [START_REF] Lemaitre | A course on damage mechanics[END_REF] and leads to a damage-independent driving force. Since interactions are mediated through the non-local damage variable d, this would amount to consider a local damage model: Heterogeneities would evolve independently from each other, and the material behavior would follow that of a single element. 2 The order of the polynomial, if too low, migh also result in a singular behavior. If the stiffness evolves linearly with d, we recover the particular case a = 0 that prevents interactions between elements. If a second-order polynomial is chosen, the element stability remains constant during damage, so that the acceleration of the damage growth observed experimentally is not captured [START_REF] Guarino | An experimental test of the critical behavior of fracture precursors[END_REF][START_REF] Fortin | Acoustic emission and velocities associated with the formation of compaction bands in sandstone[END_REF][START_REF] Baro | Similarity between the compression of a porous material and earthquake[END_REF]. The general case is obtained for polynomials of order three or higher. Damage spreading is then an unstationnary process, and the macroscopic material response is characterized by common failure behaviors, namely damage localization and sudden failure.

Note that the choice of a polynomial variation of the stiffness with d does not limit the generality of our model. In the following, we will keep the function k(d) and its derivatives in the theoretical expressions so that they can be specified to other behaviors. As evidenced here for a single element, the failure behavior of the material will be shown to depend on the sign of the second derivative k (d), but not on the actual function k(d) used to describe the stiffness degradation.

Failure behavior prediction for homogeneous materials without interactions

Before investigating the failure response of heterogeneous materials, we first explore the homogeneous situation. To impose constant failure properties through the material, the noise g c (x, d) introduced in the expression Eq. ( 9) of the damage energy is set to zero, leading to

Y c ( x, d) = Y c0 [1 + ηd( x)]. (11)
Starting from homogeneous initial conditions with an intact material d( x) = 0, the damage field remains uniform all along the failure process with d( x) = d 0 . Indeed, the driving force F(d 0 ) remains homogeneous, so that damage evolution is identical in each fiber. This has several important consequences:

• Considering an homogeneous material amounts to consider a local and interactionfree model of damage where the evolution of a heterogeneity is independent of the response of the other ones. • Each material point evolves identically. Therefore, the problem of damage evolution in the specimen is reduced to a single element problem easily tractable analytically. • The damage evolution is independent of the dimension D of the specimen, so that the predicted force-displacement response is valid for any dimension.

Using the previous remarks, the damage evolution law (8) reduces to its local form

Y (d 0 ) = Y c (d 0 ). ( 12
)
Combining the above equation with the expressions ( 7) and ( 11) of the energy release rate and the damage energy, the macroscopic displacement and force as a function of the level of damage follow

       ∆ = -2Y c0 (1 + ηd 0 ) k (d 0 ) F = ∆L D k(d 0 ) (13)
where the force evolution has been derived from its general expression given in Eq. ( 1). The force-displacement response of the homogeneously damaged material is here a parametric curve where the damage variable plays the role of the parameter increasing from d = 0 until the first fiber starts to damage up to d = 1 at full failure when the material cannot sustain any loading. The normalized macroscopic mechanical response is represented in Fig. 2(b) for non-hardening (η = 0) materials for the three values of a also used in Fig. 2(a).

Normalization is here achieved by dividing the displacement by ∆ br 0 and force by F br 0 which correspond, respectively, to the force and displacement at the end of the elastic regime for a brittle material (a = 0), hence ∆ br 0 = 2Y c0 /k 0 and F br 0 = 2k 0 Y c0 . As shown previously, the parameter a controls the fiber stability so that a > 0 results in an extended softening regime characterized by a slowdown of the damage dynamics as complete fracture is approached.

d 0 0 0.2 0.4 0.6 0.8 1 k/k 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 a = -0.3 a = 0 a = 0.3 (a) ∆ 0 /∆ br 0 0 0.5 1 1.5 2 F 0 /F br 0 0 0.2 0.4 0.6 0.8 1 1.2 a = 0.3 a = 0 a = -0.3 (b)
Figure 2. Effect of the stiffness parameter a introduced in Eq. ( 10): (a) Variations of the stiffness with the level of damage; (b) Normalized Force-Displacement responses of the homogeneous non-hardening (η = 0) material. The normalization constants ∆ br 0 and F br 0 correspond, respectively, to the displacement and force at the end of the elastic regime when a = 0.

In presence of hardening, the final failure can be delayed, as shown in Fig. 3 where the curves are normalized by the load and displacements at the end of the elastic regime, ∆ 0 = 2Y c0 k 0 (a + 1)

and

F 0 = L D 2k 0 Y c0 (a + 1)
. The snapback, which position is indicated by a dot, is observable for all the values of η explored here, even though less discernible as the hardening parameter is large. By tuning the value of η, its position can be shifted towards higher critical displacements and pushed away from the peak position. For simulating damage evolution in heterogeneous materials, we will specify our model to a = -0.3 and η = 8 so that extensive damage takes place before this critical displacement is reached. This will allow us to explore the transition towards catastrophic failure through progressive damage. Note that the type of mechanical response observed in Fig. 3 is comparable with the one observed for an array of brittle fibers loaded at the interface between two rigid blocks [START_REF] Delaplace | Damage cascade in a softening interface[END_REF]. This is consistent with our approach that considers that each elasto-damageable element represents a mesoscopic volume that contains an ensemble of smaller brittle elements.

For a homogeneous material, the snap-back position observed in Fig. 3 can be determined from the variations of the driving force with damage. If, at constant

∆/∆ 0 0 0.5 1 1.5 2 2.5 F/F 0 0 0.5 1 1.5 η = 0 η = 3 η = 6 η = 8
Figure 3. Effect of the hardening parameter η on the normalized mechanical response of the homogeneous material with a = -0.3. The normalization constants ∆ 0 and F 0 correspond, respectively, to the displacement and force at the end of the elastic regime. The snap-back instability along the force-displacement curve is indicated by a dot.

displacement, the difference Y -Y c between the damage driving force and the damage resistance increases with d, the failure process is unstable. For the criterion ( 12) to be satisfied, the macroscopic displacement must be decreased, leading to the characteristic snap-back behavior shown in Fig. 3. Hence, the snap-back condition writes as

∂(Y -Y c ) ∂d ∆ = 0 (14) which gives k (d sb ) k (d sb ) = η 1 + ηd sb (15)
This equality provides the critical level of damage at which snap-back occurs. For example, considering our chosen stiffness expression (10) with exponent γ = 3, it writes as

d sb = 1 η   1 -η 2 1 + a 3a -1   for a < 0 (16)
from which one can obtain the snap-back displacement using Eq. ( 13). As expected, when a > 0, the condition (15) required for snap-back instability is never reached. From the condition d sb ≤ 1, one can also discuss the existence of the snap-back instability for a ≤ 0. For -1 ≤ a ≤ -1 /4, the force-displacement response of the homogeneous material always display a snap-back, irrespective of the value of the hardening parameter η. On the contrary, when -1 /4 < a ≤ 0, this instability takes place only for sufficiently low hardening parameter η ≤ -6a 1 + 4a .

Interestingly, even though the mechanical response of the material depends on the value of a and in particular on its sign as illustrated in Fig. 2(b), the total mechanical energy required for breaking the specimen, which can be inferred from the area under the force-displacement curve, is constant. Indeed, the total energy stored in the material prior to the application of the loading and available for dissipation during damage growth is independent of a and follows E d tot = (1 + η/2)Y c0 as derived from the expression (11) of the damage energy. This is a remarkable property of the non-local damage model proposed here that allows to compare peak load and more generally macroscopic failure response of materials with different microscopic failure behavior for a constant total fracture energy. For example, we see in Fig. 2(b) that a non-hardening homogeneous material made of stable elements (a > 0) will resist to a larger force than if it was made of unstable elements (a < 0). However, the latter will sustain larger deformations.

4 Damage spreading and failure of heterogeneous materials: Numerical study

Interaction between elasto-damageable elements

We now explore the influence of heterogeneous properties of materials on their quasi-brittle failure behavior. For a homogeneous material, our approach reduces to a local damage model that results in a smooth evolution of the homogeneous damage field up to the snack-back instability taking place at the critical damage level given in Eq. ( 16). The introduction of material heterogeneities drastically alters this behavior as it allows interactions controlled by the function α to take place. Indeed, a heterogeneous distribution of damage is required to get d = d, and so Y = Y (see Eqs. ( 6) and ( 7)). We will see in Section 5.1 that the presence of a non-local damage driving force allows for the description of the load redistribution that takes place in the material after a damage event. This redistribution mechanism may result in cascades of failure events that build spatial and temporal correlations in the damage spatio-temporal structure, as the ones evidenced in Fig. 1(b)).

The interaction function α controls the load redistribution and, as shown in the following, controls the material failure behavior too. Its effect is explored by choosing the form

α( x) = α 0 exp - 1 2 x cos κ 2 x ( 17
)
where x represents the distance between elements. α 0 is a normalization constant fixed by the condition Σ α( x)d x = 1 that ensures energy balance during damage growth. The parameters and κ can be tuned to vary the spatial structure of the driving force redistribution after failure events. is an internal length that controls the range of the interactions. We consider situations where is larger than the heterogeneity size ξ. κ, defined as a positive constant, controls the shape of the redistribution function (see Fig. 4). Note the change of sign of the interaction functions with κ > 0. In the twodimensionnal version of the model, the interaction function solely depends on the distance between elements, irrespective of their orientation within the (x, y)-plane. 17)). For practical pruposes, one introduces the parameter κ c = 0.58.

Numerical implementation

Similarly to crack propagation problems, damage growth requires the equality Y = Y c between the rate of energy released and the rate of energy dissipation.

There are several ways to deal with this type of criterion in order to predict material evolution. One of them is to regularize this criterion by assuming an over-damped dynamics ∂d ∂t ( x, t) = µF( x) and use a very large mobility µ [START_REF] Gao | A first-order perturbation analysis of crack trapping by arrays of obstacles[END_REF][START_REF] Bonamy | Crackling dynamics in material failure as the signature of a self-organized dynamic phase transition[END_REF]). An alternative is to introduce a weakly rate-dependent critical energy release rate so that the damage evolution law takes the form of a kinetic relation Y ( x) = Y c x, ∂d ∂t that can be solved assuming the irreversibility condition ∂d ∂t > 0 [START_REF] Xia | Adhesion of heterogeneous thin films: Ii. adhesive heterogeneity[END_REF]. Here, we take advantage of the quasi-static loading conditions, and adopt a discrete dynamics inspired by [START_REF] Schmittbuhl | Interfacial crack pinning: effect of nonlocal interactions[END_REF] based on the following procedure: The imposed displacement is increased until the failure criterion ( 8) is reached by one element. The damage level of this element is then increased by δd 1 and the spatial distribution of driving forces in the material is recalculated using Eqs. ( 6) and ( 7). The redistribution of driving force results in additional damage and subsequently additional driving force redistribution until this cascade process stops when the driving force is below its critical value everywhere in the material. Those three approaches result in the same avalanche-like dynamics of the damage growth illustrated in Fig. 1(b): Damage increases by bursts that are localized both in space and time, characterized by power law statistics but that are characteristic of fracture processes in disordered materials [START_REF] Herrmann | Statistical Models for the Fracture of Disordered Media[END_REF][START_REF] Alava | Statistical models of fracture[END_REF][START_REF] Bonamy | Intermittency and roughening in the failure of brittle heterogeneous materials[END_REF]. Despite the presence of strong fluctuations during the growth of damage, we will see that main features of the macroscopic failure response of the material can be understood without taking into account this intermittency.

To simulate damage growth, we discretize a solid of size [0, L] D by considering discrete elements corresponding to individual heterogenities, so that the total number of elements is N = L/ξ. As shown later, a finer mesh, corresponding to more elements into a single heterogeneities, does not affect significantly the calculated material response, even in the softening section of the stressstrain curve. The following numerical algorithm is used: Starting from an initially intact material (i) the displacement ∆ is increased until one element satisfies to the condition Y ( x) = Y c ( x); (ii) the damage level at this point is increased by δd 0 ; (iii) the distribution of non-local damage driving force and fracture energy are recalculated; (iv) The step (ii) is repeated for all the elements verifying F( x, ∆) ≥ 0; (v) the steps (iii) and (iv) are repeated until no more elements verify the damage criterion; (v) We go back to step (i) and increase the external driving ∆ again until new damage events take place. The simulation stops when the damage variable is equal to one everywhere. This discrete dynamics allows the exploration of large specimens, typically N = 10 4 heterogeneities in 1D and N × N = 2.10 2 × 2.10 2 heterogeneities in 2D for a rather low computational time. As a result, several speciments with the same material properties, but different realizations of the disorder can be runned to achieve accurate average material response. We did not observe any significant effect of the value of N as long as the internal length is small with respect to the specimen size L. Note also that to remove edge effects, we use periodic boundary conditions.

To calculate the damage driving force defined by Eqs. ( 6) and ( 7), we use Fourier transforms since convolution product becomes simple products. First, the non-local damage parameter is obtained from

d( x) = F -1 α( q) × d( q) (18)
where F -1 denotes the inverse Fourier transform and α is the Fourier transform of the function α. The non-local damage parameter is then used to compute the stiffness and the driving force is obtained from the expression

Y ( x) = - 1 2 F -1 α( q) × k (d) . ( 19
)
The variations of the damage energy are also discretized by introducing a small but finite damage increment δd 0 that characterizes the amplitude of a single damage event taking place when the condition Y = Y c is reached. After each damage event, the fracture energy of the element is redrawn using the following expression

Y c ( x, n 0 ) = Y c0 [1 + g c ( x, n 0 ) + n 0 η] (20)
where n 0 = d/δd 0 is the total number of damage events experienced by the element since the beginning of the test, η = ηδd 0 is the modified hardening parameter and g c ( x, n 0 ) is a random number drawn in a uniform distribution within the interval -√ 3σ < g c < √ 3σ. In the following, we use δd 0 = 0.005 for the damage increment. This value is observed to have a minor effect on the simulation results as long as taken small with respect to unity.

Numerical results

We focus first on 1D specimens while 2D specimens, which display qualitatively similar behaviors, are left for Section 5.3. The normalized macroscopic response of weakly heterogeneous materials (σ = 0.001) is presented in Fig. 5(a) for different values of the shape parameter κ and compared with the response of the homogeneous material, reminiscent of a local damage evolution law. The force and displacement are normalized by their values at the end of the elastic regime. In order to smooth out the fluctuations observed from one material to another, the response of the heterogeneous material is shown after averaging over twenty different realizations of the disorder.

For small imposed displacements ∆ ∆ c , where ∆ c indicates the displacement at failure, the response of the heterogeneous material is identical to the one of the homogeneous material. As ∆ gets close to ∆ c , its behavior deviates slightly from the homogeneous response, as illustrated in the inset of Fig. 5(a) for κ = 2κ c . Finally, for ∆ = ∆ c indicated by the dots in Fig. 5(a), we observe a sudden drop of the force that reveals a sudden failure of all the material elements. For κ = 0, no deviation is observed until the catastrophic failure takes place at the snap-back instability. Hence, varying the value of κ shifts the failure instability that lies (i) at the snap-back point for κ = 0, (ii) between the peak force and the snap-back point for κ = 2κ c , and (iii) in the pre-peak region for κ = 4κ c . Therefore, at low disorder level, the shape of the interaction function controls the occurrence of the failure instability while maintaining a global behavior close to the one of a material with a local failure response. To determine the onset d h (resp. ∆ h ) of the deviation to the homogenous response, we introduce the intensity of the damage spatial fluctuations A first observation made from Fig. 6 is that the thresholds vary weakly with the material disorder: Low variations of the values of d h and d c are observed while the fluctuation amplitude σ of damage energy is varied from 0.1 % to 20 % of its average value. Another interesting observation is that the critical damage levels are rather independent of the internal length. As a result, the material response is essentially governed by the value of κ and two distinct regimes emerge:

∆/∆ 0 0 0.5 1 1.5 2 2.5 F/F 0 0 0.5 1 1.5 κ = 0 κ = 2κ c κ = 4κ c Homogeneous 1.9 2 1.2 1.3 (
σ d = ( d(x) 2 -d(x) 2 ) 2 (
• For κ ≤ κ c , catastrophic failure takes place at the snack-back instability (d c d sb ), and is not preceded by a phase where the response of the material deviates from the homogeneous material behavior. • For κ > κ c , the thresholds strongly depend on the shape of the redistribution function. A carefull comparison between Fig. 6(a) and (b) reveals that the deviation to the homogeneous response takes place close, but prior to complete failure (d h d c < d sb ). In this regime, failure occurs earlier when the value of κ is increased, taking place even before peak load for κ 3κ c

To emphasize the negligible influence of the spatial discretization on the calculated material response, we study the variations of the critical loading at failure as a function of the number N of heterogeneities in the material and the number of points N ξ used to discretize one single heterogeneity. In the latter case, the value of damage energy is drawn in a uniform distribution of variance σ at the center of each heterogeneity, and then interpolated in between so Y c ( x) is a smoothly varying function of the position. In both cases, we obtain rapid convergence as shown in Fig. 7 for κ = 0, = 2 ξ and σ = 0.001. Here, ∆ conv denotes the failure loading obtained for the finest mesh corresponding to N = 10000 and N ξ = 10. To further characterize the behavior of the heterogeneous materials, we explore the spatial structure of the damage field close to failure. It is convenient to introduce the distance to failure

δ = d c -d d c . ( 22
)
This parameter δ ranges from unity at the end of the elastic regime to zero at catastrophic failure. The power spectrum

P δ (q) = | dδ (q)| 2 N ( 23
)
is computed at different distance δ from failure. Here, dδ (q) is the Fourier transform of the damage field and q is the Fourier mode corresponding to the wavenumber λ = 2π/q. The power spectrum is used to characterize the spatial structure of the damage field and its evolution. If the damage fluctuations are random, all wavenumbers contribute equally and the power spectrum is flat.

On the contrary, if the damage spreads favoring a certain mode, the spectrum shows a peak. A typical power spectrum evolution is shown in Fig. 8(a) for κ = 2κ c and = 5 ξ. As highlighted in the inset, far from failure (δ ≥ 0.55), the power spectrum is fairly flat which indicates random fluctuations of the damage field. On the contrary, as the specimen gets closer to failure ( δ ≤ 0.3), a peak emerges in the Fourier spectrum, indicating that damage spatially organizes and develops over the mode q c . Its amplification when δ → 0 suggests that a characteristic length scale 2π/q c emerges from the damage spatial structure. This observation is consistent with the deviation from the homogeneous material response brought out in Fig. 6 close to failure.

Figure 8(b) shows the effect of κ on q c for two values of the internal length, = 5 ξ and = 10 ξ, and for a large range of disorder strength 0.001 ≤ σ ≤ 0.2. The two regimes identified from the variations of the thresholds d h and d c with κ are also revealed by the value of the characteristic mode as for κ ≤ κ c , the power spectrum does not reveal any characteristic length scale whereas for κ > κ c , a peak emerges and its position q c increases with κ. Furthermore, for κ κ c the characteristic mode q c increases linearly with κ, following the variations of the mode q α = κ/2 introduced through the interaction function α (see Eq. ( 17)). In summary, the numerical investigation of the damage evolution in heterogeneous materials with interacting elements shows that their response can be captured by the homogeneous interaction-free material behavior up to some critical loading where the two responses deviate from each other, quickly followed by an abrupt failure. Both onsets can take a broad range of values, from the snap-back position to the pre-peak regime, depending on the value of the shape parameter κ of the interaction function. Two distinct regimes can be evidenced: • For κ < κ c , the damage field does not display a characteristic mode until unstable failure takes place close to the snap-back instability without any deviation to the homogeneous material response. • For κ > κ c , we evidence a deviation to the homogeneous response, quickly followed by an abrupt failure. Evolution towards failure is characterized in this regime by the rapid growth of a characteristic mode in the damage field.

These observations call for the following questions: Can we explain deviation from the homogeneous response and failure and predict their onset? Where does the characteristic length scale of the damage field close to failure emerge from? What is the origin of the two regimes displaying distinct failure behaviors as the interaction function is varied?

5 Damage spreading and failure in heterogeneous materials: Theoretical study

Stability analysis of the homogeneous damage states

To address the questions raised by the numerical study, we carry out a stability analysis of the states reached by the homogeneous material during its evolution. Here, we do not restrict our analysis to 1D specimens, and explore the stability of the solutions obtained in Sec. 3 for any dimension. We assume small variations of the damage field around a homogeneous value and decompose it into two contributions: A uniform contribution d 0 solution of the homogeneous problem and a heterogeneous contribution δd( x)

d( x, ∆) = d 0 + δd( x) ( 24
)
where δd is assumed to be very small with respect to d 0 . To determine the evolution of the perturbation δd, the total damage driving force

F = Ȳ -Y c is also decomposed into two contributions F[d( x, ∆)] = F (0) (d 0 ) + F (1) [δd( x)] (25)
where the constant term F (0) is the contribution of the homogeneous part of the damage field while the heterogeneous part of the driving force F (1) results from the perturbation δd. The expression of each of these terms is determined by linearizing the expressions ( 6) and ( 11) of the nonlocal energy release rate and the fracture energy, respectively. The first step is to decompose the nonlocal damage parameter into two terms

d( x, ∆) = d 0 + (α * δd)( x) (26)
which leads to the following expansion of the stiffness derivative

k [d( x)] = k (d 0 ) + k (d 0 )δd( x). (27)
In this expression where only first order terms in δd have been kept, we have introduced the non-local perturbation parameter δd( x) = (α * δd)( x). Using Eq. ( 27) into Eq. ( 6), one obtains the first order expansion of the non-local energy release rate at a given applied displacement ∆

Y ( x) = Y (d 0 ) - ∆ 2 2 k (d 0 )(α * δd)( x). ( 28
)
Since the fracture energy also decomposes into two contributions

Y c ( x, d) = Y c (d 0 ) + Y c0 ηδd( x), (29)
one gets the zeroth and first order contributions of the total driving force

       F (0) (d 0 ) = Y (d 0 ) -Y c (d 0 ) F (1) [δd( x)] = - ∆ 2 2 k (d 0 )(α * δd)( x) -Y c0 ηδd( x). (30)
The zeroth order provides the equilibrium condition

F (0) = 0 ⇒ Y (d 0 ) = Y c (d 0 ). ( 31
)
It corresponds to the evolution equation of the homogeneous material obtained in Sec. 3 and provides the typical damage level d 0 (∆) in the material as a function of the external loading as given in Eq. ( 13).

We now focus on the perturbative term F (1) of the driving force. Its expression provides physical insights on the interaction function α introduced in the model definition. Consider a perturbation δd( x) ∼ δ( x -x 0 ) of the damage field resulting from a damage event located in x 0 , δ( x) being the Dirac function. The resulting perturbation on the field of damage driving force can be deduced from Eq. ( 30) and follows F (1) (r) ∼ (α * α)(r) where r = x-x 0 > 0 is the distance from the damaged element. For κ = 0, this term is positive, meaning that a damage event in the material is followed by an increase of the driving force everywhere else. On the contrary, for κ > 0, α and so F (1) change of sign with r, indicating that some elements will be unloaded while some others will be loaded further. The spatial structure of the load redistribution is actually set by the convolution of α with itself which behaves qualitatively as the function α represented in Fig. 4. Therefore, the interaction function α involved in the definition of the non-local damage parameter d plays a central role in our model by describing how the driving force is redistributed in the material after a failure event. Note that the damaging element itself is always unloaded, irrespective of the value of κ. Indeed, the sign of the load redistribution in x 0 is provided by F (1) (0) = ∆ 2 |k ( ∆ )|/2 -Y c0 η which remains negative before the snap-back instability, as shown by combining Eqs. ( 13) with ( 15). This illustrates the central role played by the hardening behavior of the material that prevents unstable failure of individual elements before the snap-back point.

We now use the expression of the driving force perturbation F (1) [δd( x] to determine the material stability. We consider harmonic perturbations

δd( x) = δd 0 (1 + cos( q 0 • x)) (32)
of mode q 0 and small positive amplitude δd 0 . This expression ensures the irreversibility condition δd ≥ 0 imposed to the damage evolution. Inserting this expression in Eq. ( 30) provides the driving force redistribution

F (1) [δd( x)] = δd 0 cos( q 0 • x) * - ∆ 2 2 k (d 0 )α( x) * α( x) -Y c0 ηδ( x) +δd 0 × - ∆ 2 2 k (d 0 ) -Y c0 η (33)
as the damage increases of δd( x). As the former expression comprises convolution products, it is more convenient to express it in Fourier's space

F(1) [ δd( q)] =δd 0 δ( q -q 0 ) × - ∆ 2 2 k (d 0 )α 2 ( q) -Y c0 η + δd 0 δ( q) × - ∆ 2 2 k (d 0 ) -Y c0 η . ( 34
)
Introducing the function

G( q, d 0 ) = - 1 2 ∆ 2 k (d 0 )α 2 ( q) -Y c0 η, (35)
the driving force perturbation takes the simple form

F(1) [ δd( q)] = δd 0 G( q, d 0 )δ( q -q 0 ) + G( 0, d 0 )δ( q) . ( 36
)
where the normalization condition α( 0) = 1 imposed to the interaction function have been used together with Eq. ( 13) relating ∆ to d 0 . Taking the inverse Fourier transform of the former expression provides the driving force perturbation

F (1) [δd( x)] = δd 0 G( q 0 , d 0 ) cos( q 0 • x) + G( 0, d 0 ) . ( 37
)
This expression provides the condition for a homogeneous damage growth: If F (1) ( x) is higher at the most damaged zones, i.e. for x• q 0 = 0, ±2π, ±4π..., the perturbation will grow and hence deviation from the homogeneous response will occur. On the contrary, if the driving force is higher on the elements with smaller damage levels, i.e. for x • q 0 = ±π, ±3π..., these elements will damage first and perturbation to the homogeneous damage field will be smoothed out.

Noticing that only the first term of F (1) in Eq. ( 37) varies with x and using the above considerations, we conclude that

G( q 0 , d 0 ) > 0 ⇒ Damage localization. ( 38
)
In the following, we refer to this process as damage localization, since it results into a concentration of the failure processes in some regions of the material. This condition is satisfied above a critical value of the imposed displacement, or equivalently, above some critical level of damage, for the most unstable mode q c that maximizes α. Indeed, according to Eq. ( 35), the position q c of the maximum of G( q 0 , d 0 ) is also the one of α, irrespective of the typical damage level d 0 . At some critical damage d h (resp. some critical loading ∆ h ), the function G that increases with d 0 (resp. ∆) hits zero. The perturbation then grows and leads to an heterogeneous damage field with a characteristic mode q c . On the force-displacement response, one expects this damage localization to produce a departure from the homogeneous material behavior of Eq. ( 13).

The condition (38) predicts the emergence of a heterogeneous mode of damage growth, but does not provide any clue on the abrupt failure observed in the simulations. Indeed, an unstable behavior requires an increase of the damage somewhere in the material to be accompanied by an increase of the driving force too at the same position. Considering again the positive damage increment δd = δd 0 (1 + cos( q 0 • x)), catastrophic failure takes place if F (1) ( x) is positive in at least one position. Using Eq. ( 37), this condition is first satisfied in x • q 0 = 0, ±2π..., leading to the criterion

G( q 0 , d 0 ) + G(0, d 0 ) > 0 ⇒ Unstable failure. ( 39
)
Predicting the onset d c of unstable failure requires at first the determination of the most unstable mode q c . As for localization, it is given by the position of the maximum of G( q 0 , d 0 ), so the localized mode q c is the failure mode too.

Note that the inequality G(0, d 0 ) < 0 derived from the expression (15) of the snap-back point ensures that localization always takes place prior to abrupt failure. Indeed, it implies that the localization threshold G( q c , d 0 ) > 0 has to be reached so that the failure condition G( q c , d 0 )+ G(0, d 0 ) > 0 can be fullfilled.

An important point is that, contrary to the localization prediction, the onset of failure determined from Eq. ( 39) predicts only approximately the actual failure threshold, even in the limit of weakly heterogeneous materials. Indeed, this criterion is derived from a perturbation of the homogeneous material response, hence assuming a homogeneous distribution of damage in the material. Since localization might take place prior to final failure, damage might be distributed heterogeneously when unstable failure takes place. We will see however that the predictions made following this approach describes rather well the numerical observations, especially if κ is not too large with respect to κ c , i.e. that localization does not occur much earlier than failure so that the assumption of a relatively homogeneous damage field at failure is not violated.

From criteria ( 38) and ( 39), one can now predict the specimen evolution. Two cases need to be distinguished:

• The most unstable mode corresponds to the homogeneous perturbation q c = 0. Damage localization is here simultaneously accompanied by unstable failure since both criteria simplify to G(0, d 0 ) > 0 that is the snap-back instability condition. As a result, the material response is the one of the homogeneous material where catastrophic failure takes place without any prior localization phase at the snap-back point

d sb = d h = d c .
• The mode q c that maximizes α( q) is different from the homogeneous mode. This implies α( q c ) > α( 0), or equivalently, G( q c , ∆) > G(0, ∆) using the expression (35) of G. This inequality implies that localization is strictly anterior to failure and that both occur before the snap-back instability of the homogeneous material, leading to d h < d c < d sb . Interestingly, this regime is characterized by the growth of localized damage modes that precede abrupt failure.

The prediction of the onsets of localization and failure are now compared with the simulation results.

Application to 1D specimens

The predictions of the linear stability analysis are first compared with the numerical results for 1D specimens for which all calculations can be performed analytically. The first step is the determination of the critical mode q c that maximizes the redistribution function G(q, d 0 ). According to Eq. ( 35), this amounts to determine the maximum of α(q). The Fourier transform α(q) = +∞ -∞ e -iqr α(x)dx of the interaction function α given in Eq. ( 17) reads

α(q) = (1 + κ 2 )(1 + κ 2 + (2 q) 2 ) (κ 2 + 1) 2 + 2(2 q) 2 (1 -κ 2 ) + (2 q) 4 . ( 40
)
from which G(q, d 0 ) can be deduced. Figure 9 represents its variations as a function of q/q c for different damage levels for the parameter values κ = 2κ c , = 5 ξ and η = 8. The position of the maximum, q c 0.55/ for these parameters, does not vary with the loading. On the contrary, the maximum value G(q c , d 0 ) of the redistribution function does increase and eventually hits zeros at the damage level d h 0.466.

The analysis of the function α defines the critical value κ c = 1/ √ 3 that distinguishes two different behaviors:

• For κ ≤ κ c , the Fourier transform of the interaction function is maximum for q c = 0. Failure is expected to take place at the snap-back point, without any prior damage localization. • For κ > κ c , the maximum is reached for q c > 0 that follows

q c = 1 2 2κ √ 1 + κ 2 -(1 + κ 2 ). ( 41
) q/q c 0 2 4 6 8 G(q, d 0 ) -10 -8 -6 -4 -2 0 2 d 0 =0.100 d 0 =0.200 d 0 =0.300 d 0 =0.400 d 0 =0.466 d 0 =0.517
Figure 9. Prediction of the localization for a 1D specimen with κ = 2κ c , = 5 ξ and η = 8. This graph shows the amplitude of the driving force perturbation for a sinusoidal perturbation of mode q of the damage field (Eq. ( 35)). Growing modes correspond to positive amplitude, so localization takes place at the damage level d 0 0.466 for the mode q c given in Eq. ( 41).

The criterion G(q c , d h ) = 0 of Eq. ( 38) provides the damage level at localization

d h = 1 η α2 (q c ) 2α 2 (q c ) -1   1 -η 2 2α 2 (q c ) -1 α4 (q c ) 1 + a 3a -1   (42)
with

       α(q c ) = (1 + κ 2 ) 3/2 4κ(1 + κ 2 -κ √ 1 + κ 2 ) for κ > κ c α(q c ) = α(0) = 1 for κ ≤ κ c . (43)
For κ ≤ κ c , the onset d h becomes independent of the redistribution function and is equal to the position of the snap back instability (Eq. ( 16)): We recover the behavior of the homogeneous material studied in Sec. 5.1. On the contrary, when κ > κ c , the most unstable mode q c > 0 is not homogeneous, so the shape of the redistribution function does affect the onset d h of localization through the parameter κ. The internal length introduced in the interaction function sets the wavelength of the unstable mode 1/q c ∼ , but does not affect the localization threshold. The introduction of an internal length in the redistribution function reflects on the characteristic wavevector of the instability and hence sets the width of the 'localization bands' emerging close to failure. The scaling 1/q c ∼ was expected as the internal length is the only relevant length scale -the system size L can be taken infinitely large while the heterogeneity size ξ does not affect the stability analysis that considers a homogeneous solid.

These predictions are now compared to the numerical results. Figure 6(a) shows the variations of the localization threshold d h with κ (solid line). The theory captures nicely its decay in the regime κ > κ c , irrespective of the disorder amplitude.3 . In Fig. 8(b), the comparison made on the wavenumber q c of the localization mode for two values of internal length and several disorder levels is also successful, confirming especially the scaling q c ∼ 1/ .

The thresholds are also compared to the mode q α = κ/(2 ) introduced in the interaction function expression. Interestingly, the localization mode q c is not equal to the characeristic mode: q c is indeed smaller than q α for the shape parameters κ κ c , and tends asymptotically towards q α for κ κ c .

We now tunr to the prediction of catastrophic failure. To determine d c , the failure criterion G(q c , d c ) + G(0, d c ) = 0 of Eq. ( 39) is applied using the critical mode q c determined previously. This gives

d c = 1 + α2 (q c ) 2η α2 (q c )    1 - 1 + a 3a 2η α(q c ) 1 + α2 (q c ) 2 -1    (44)
where the expression of α(q c ) is provided in Eq. ( 43). If κ ≤ κ c , we verify that d c = d h = d sb . In that case, perturbation growth and unstable failure occur simultaneously at the instability point of the interaction-free material. In the other regime, when κ > κ c , we obtain d h < d c < d sb , in agreement with the analysis performed for any dimension in Sec. 5.1.

The predicted failure thresholds are shown in solid line in Fig. 6(b). The equation (44) slightly overestimates the onset of failure found in the simulations, in particular for κ κ c . This was expected since the theory relies on a perturbation analysis of the homogeneous damage state, while localization already took place when the failure threshold is reached. Nevertheless, we note that this approach captures qualitatively well the behavior observed in the simulations, and in particular the decrease of d c with the interaction parameter κ. An improved prediction of the onset of failure based on a global energy minimization will be provided in Sec. 5.4.

Application to 2D specimens

We now apply the stability analysis of Sec. 5.1 in the context of 2D specimens. In order to apply the criteria derived in Eqs. ( 38) and (39) for localization and failure that involve α(q), we notice that the interaction function α( x) = α(r) is a function of the distance r = x only, so its Fourier transform α( q) = α(q) is a function of the norm q = q of the wavenumber. Hence, the stability analysis of a 2D specimen reduces to a 1D problem where all the perturbation modes δd( x) = δd 0 (1 + cos( q 0 . x)) of same norm q 0 share the same behavior, irrespective of their direction. The localization and failure thresholds are then determined following the following two steps procedure:

(i) The norm q c of the most unstable modes is determined from the maximum of G(q 0 , d 0 ). According to Eq. ( 35), this amounts to determine the maximum of the Fourier transform of the interaction function α( q) = +∞ -∞ α( x)e -i q. x d x.

Using the isotropy α( x) = α(r) of the interaction function, the Ankel transform gives α(q) = 2π +∞ 0 α(r)J 0 (qr)rdr where J 0 is the Bessel function of the first kind of order zero. This integral is computed numerically to determine its maximum α(q c ) and the norm q c of the unstable modes. (ii) From the values of q c and α(q c ) calculated numerically, the critical damage levels for localization and failure are determined from the criteria G(q c , d h ) = 0 and G(q c , d c ) + G(0, d c ) = 0, respectively. Since the function G introduced in Eq. ( 35) depends on the dimension of the specimen only through α(q), the expressions ( 42) and ( 44) of both thresholds remain unchanged.

This procedure is applied for different values of κ. For 2D geometries, the normalization constant α 0 of Eq. ( 6) is non-negative only for κ < 1. Above this value, the redistributions defined by the model are non-physical and hence we limit our analysis to the range 0 ≤ κ < 1. The obtained variations of the critical wavenumber q c with the shape parameter κ of the redistribution function is represented in Fig. 10(a). 4 As for a 1D geometry, one obtains two regimes: For κ ≤ κ 2D c , where κ 2D c 0.42, the most unstable mode is the homogeneous perturbation q c = 0 while for κ > κ 2D c , the norm q c of the most unstable wavevectors increases continuously with the parameter κ.

The critical damage values for localization (solid line) and failure (dashed line) deduced from Eqs.( 42) and ( 44) and the value of α(q c ) are represented as a function of κ in Fig. 10(b) for a toughening parameter η = 8 and a = -0.3. Similarly to 1D specimens, we observe a transition from a response similar to the one of a homogeneous interaction-free material for low values of κ, to another regime where the deviation to the homogeneous material response is observed prior to catastrophic failure. The later intervenes prior to the snapback instability observed in the homogeneous material response, resulting in d h < d c < d sb 0.78. These theoretical predictions are now compared with direct numerical simulations of the damage spreading in 2D specimens with η = 8 following the procedure described in Sec. 4.2. Irrespective of the value of the disorder level σ, the theory captures well the decrease of the localization threshold d h with κ, shown by dot symbols in Fig. 10. The decrease of the failure threshold d c (cross symbols) is also described, but only qualitatively for the same reasons invoked in the former section in the context of 1D geometries.

Even though failure behaviors are qualitatively similar in one and two dimensions, we would like however to emphasize some quantitative differences. First, the threshold value κ 2D c 0.43 in 2D is significantly smaller than its value in 1D for which κ c = 1 / √ 3 0.58. Second, the critical damage values for both localization and failure go to zero for κ → 1, whereas in 1D these quantities go to zero for κ → +∞. This latter case where d h = d c = 0 corresponds to a failure taking place right after the elastic regime while, on the contrary, κ = 0 ensures the existence of an extended softening regime until abrupt failure takes place at the snap-back point with d h = d c = d sb . This observation illustrates the centrale role played by the spatial structure of the stress redistribution on the emergent failure behavior of the material, as changing the shape parameter of the interaction function from zero to one shifts the response of the material from elasto-damageable to purely brittle. This also highlights the collective nature of the damage growth process in quasi-brittle solids: Although the properties of its constitutive elements are kept simple, the material displays a wide range of failure behaviors that emerge from the their interactions controlled through the load redistribution process.

Prediction of damage localization and failure from a global energy minimization

We now explore an alternative approach for predicting the failure behavior of elasto-damageable materials. In the previous section, we investigated the stability of the damage growth process through the analysis of the distribution of damage driving force at the local scale in the material. It allowed us to determine the mode q c over which the damage field builds up and subsequently, the critical load at localization for which this mode starts to grow. Despite the assumption of a slightly heterogeneous damage state, we applied this approach to predict the onset of failure that takes place after damage has localized, and showed that it provides nevertheless a reasonable estimate of the failure threshold measured numerically.

However, in an attempt to provide a better estimate of the failure threshold, we propose here a complementary approach based on the minimization of the global energy of the damaged material. The other motivation is to gain physical insights on the transitions revealed by the stability analysis, as an energetic approach might shed light on their nature. Here, we take advantage of the formulation of our damage model that has been derived from energy conservation, and determine, at each imposed loading conditions, the optimal distribution of damage in order to minimize the global specimen energy comprising the stored elastic energy and the damage energy.

For the sake of simplicity, this analysis is applied to 1D materials only. To resolve the onset of localization, we determine whether or not a small perturbation of the damage field is energetically favorable, i.e. results in a lower system energy than a homogeneous state. We assume at some prescribed displacement ∆ that the homogeneous damage field is perturbed by an harmonic mode q 0 : d(x) = d 0 + δd 0 cos(q 0 x) and determine the couple (d 0 , δd 0 ) that minimizes the total energy, calculated as

E ∆ (d 0 , δd 0 ) L = 1 2 ∆ 2 k(d 0 )+Y c0 d 0 + 1 2 ηd 2 0 +δd 2 0 1 8 ∆ 2 k (d 0 )α 2 (q 0 ) + 1 4 Y c0 η .
(45) This equation derives directly from the expression (4) of the total energy linearized to the second order in δd 0 where the work of the external force is not taken into consideration as we seek to determine the damage distribution at a fixed displacement. Noticing that the first two terms are positive and that the prefactor of δd 2 0 equals to -G(q 0 , d 0 )/4, we conclude that δd 0 > 0 results in an increase of energy as long as the redistribution function G(q 0 , d 0 ) is negative. In that regime, the energy is hence minimized for δd 0 = 0 and d 0 verifies ∂E ∆ /∂d 0 | (d 0 ,0) = 0 that corresponds to the equilibrium equation Y (d 0 ) = Y c (d 0 ) derived for a homogeneous material. We therefore retrieve the localization con-dition (38) obtained from the stability analysis that the damage field remains homogeneous as long as G(q 0 , d 0 ) < 0. This conditions provides the loading ∆ h at localization at which the energy is no longer minimized by a homogeneous damage field, but also the localization mode q c that maximizes G.

We now consider loadings ∆ > ∆ h for which a localized damage field develops. In that case, we must ensure the irreversibility of the damage perturbation that can be done by considering a perturbation d(x) = d 0 + δd 0 (1 + cos(q 0 x)). Since the localized mode has been previously determined, we consider q 0 = q c . The total energy follows

E ∆ (d 0 , δd 0 ) = A(d 0 ) -B(d 0 )δd 0 + C(d 0 )δd 2 0 -D(d 0 )δd 3 0 ( 46
)
where the coefficients A, B, C and D are provided in appendix B. The extrema of E ∆ can be determined from the following system of equations

               ∂E ∆ ∂d 0 δd 0 = 0 ∂E ∆ ∂δd 0 d 0 = 0. ( 47
)
For ∆ > ∆ h , the homogeneous solution (d 0 , 0) is no longer a minimum of E ∆ . Instead, a heterogeneous damaged state with a positive perturbation amplitude We first test the accuracy of our calculation. Figure 11(a) shows the energy landscape caculated from its expression (46) as a function of δd 0 at different loading increments, from δ∆ = 0 at localization to δ∆ c = 0.088 at failure, using the parameters η = 8 and κ = 2κ c . The level of damage d min 0 = d h + T δ∆ has been fixed for each loading to the value minimizing E ∆ (d 0 , δd 0 ), as determined from the procedure described previously. The position of the minimum predicted from Eq. ( 48) is indicated by dots. Due to the linearization performed around the localization threshold in our calculation, the failure loading is slightly overestimated, as it can be noticed that the energy has actually no local minimum for δ∆ > 0.08. However, the agreement with the numerical determination of the minimum of E ∆ (d 0 , δd 0 ) that corresponds to the equilibrium configuration of the damage field remains excellent.

δd 0 = C(d 0 ) 3D(d 0 )   1 -1 - 3D(d 0 )B(d 0 ) C 2 (d 0 )   . ( 48
Figure 6(a) compares the predicted values of the damage level at failure for different values of κ (dotted line) with the results of the simulations performed in Sec. 4 while Fig. 11(b) compares the failure loads. Failure prediction is significantly improved. This can be explained by the difference between both methods: The energy minimization procedure allows the determination of both d 0 (∆) and δd 0 (∆) in the post-localization regime while the homogeneous solution d 0 (∆) = d h 0 (∆) stricly valid for ∆ < ∆ h only is assumed to hold after localization in the stability analysis. But overall, both approaches provide reliable estimates of the failure properties of the elasto-damageable materials, demonstrating that a linear stability analysis and a procedure based on global energy minization are both appropriate to predict localization and failure from damage models derived from thermodynamic principles. Furthermore, the energy minimization procedure sheds light on the nature of the transitions. The case κ > κ c when a stable post-localization regime does exist is considered in Fig. 12(a) that schematically represents the position of the energy extrema in the (d 0 , δ 0 ) plane. At the localization threshold, the equilibrium position bifurcates from a homogeneous to a heterogeneous damage state. When d 0 increases further, the amplitude of the localized mode increases also, starting from δd 0 = 0 at ∆ = ∆ h , revealing that the localization is a continuous transition, as also observed in the context of gradient damage models [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in the uniaxial tests with gradient damage models[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]). An unstable equilibrium position, represented by the dotted branch Fig. 6(a), co-exists with the stable one for d 0 < d c . Both branches merge together at d = d c and subsequently vanish: The material abruptly fails as it does not admit any equilibrium position.

The case κ ≤ κ c of Fig 12(b) shows a different behavior as the stable and unstable branches merge before the bifuraction to a heterogeneous equilibrium stage could take place. The material therefore fails as the damage field is still homogeneous. However, the transition to failure is of the same nature in both cases as it results from the loss of equilibrium position.

Discussion

We proposed a damage model based on a few, well identified and tunable sets of physical ingredients (disorder, interaction, hardening) that was shown to capture the main features of quasi-brittle fracture, namely localization and unstable failure. Sufficiently simple to be tractable analytically, it provides several insights on these transitions.

At the element scale, the critical load at failure emerges from the competition between the degradation of the elastic properties that promotes catastrophic failure with the toughening process that delays it. In our model, the stability of a single element is set by the parameter a that describes the variation of the stiffness with damage (a negative value a < 0 ensures an unstable damage growth) while the toughening is controlled by the parameter η. Therefore, the critical failure load is a decreasing function of |a| and an increasing function of η. An interesting observation is that hardening is required in our model to achieve an extended regime of stable damage growth. First, it should be noticed that this property does not result from the specific conditions chosen in our work, as elements arranged in series or loaded under force imposed conditions would both lead to a similar behavior. This observation is compatible with [START_REF] Girard | Failure as a critical phenomenon in a progressive damage model[END_REF]'s findings who did not include hardening in their model and who observed that localization takes place after, but close to the first damage event -they do observe a small regime of stable damage growth, but we expect it to disapear in the homogeneous material limit, as in our model too, disorder tends to delay catastrophic failure. In the model proposed by [START_REF] Zapperi | Plasticity and avalanches behavior in microfracturing phenomena[END_REF], the situation is different, as they do observe an extended stable post-peak regime. This can be explained from a different choice of damage criterion. In their simulations, they consider that an element starts to damage when the local tensile stress exceeds some critical threshold reminiscent of the material resistance. In our model, this would amount to consider a damage criterion of the type F = k(d)∆ > F c instead of an energy-based failure criterion. Using their criterion in our model would drastically increase stability, as it would result in the stability condition ∂F/∂d| ∆ = k (d)∆ = 0 that is always satisfied.

At the sample scale, the competition between stiffness degradation and hardening is more complex as it involves the interaction between the constitutive elements. As a result, macroscopic failure emerges from the collective response of these constitutive elements. Note that this point of view is different the one adopted by [START_REF] Rudnicki | Conditions for the localization of deformation in pressure-sensitive dilatant materials[END_REF] who derived a localization criterion from the response of a single element. The mechanism at the origin of these interactions is the load redistribution process. It plays a central role by selecting the macroscopic failure behavior from the material characteristics at the small scale. It can be described by the function G that provides the damage driving force redistribution δF( x) = (G * δd) ( x) resulting from a damage increment δd( x). In our model, the function G was chosen to decay exponentially fast with the distance to the damage event, as it was set by the function α involved in the definition of the nonlocal damage variable.

But what is the mechanism at the origin of the load redistribution process? Recently, Démery et al. calculated explicitely the redistribution of damage driving force in a 3D elasto-damageable material governed by 3D elasticity and a damage criterion similar to the one employed in this study. It is shown that the redistribution kernel G resembles to the Eshelby solution for the mechanical fields around a soft inclusion embedded in an infinite elastic medium (Eshelby, 1957) and displays a quadrupole symmetry G( x , θ) ∼ cos(4θ x ). Note that based on numerical simulations [START_REF] Falk | Dynamics of viscoplastic deformation in amorphous solids[END_REF][START_REF] Maloney | Amorphous systems in athermal, quasistatic shear[END_REF], a similar function had been proposed to describe the stress redistribution following a plastic event in amorphous materials ( [START_REF] Vandembroucq | Mechanical noise dependent aging and shear banding behavior of a mesoscopic model of amorphous plasticity[END_REF][START_REF] Lin | Scaling description of the yielding transition in soft amorphous solids at zero temperature[END_REF][START_REF] Bouil | Emergence of cooperativity in plasticity of solt glassy materials[END_REF].

The differences and similarities between the redistribution function chosen in our model and the one calculated in the context of 3D elasto-damageable solids are the following:

• Both redistribution funtions are strongly heterogeneous in space. In particular, the smaller the distance is from the damage event, the stronger the reloading is. However, our model considers a finite range of interactions controlled by the internal length , while the kernel derived from 3D elasticity follows a power law decay. • Both redistribution funtions show variations in the sign of the reloading that depends on the position in the damaged solid with respect to the damage event. This property is a key feature as it allows the unloading of some regions of the material while other ones are reloaded. However, our model considers redistribution kernels the sign of which changes with the distance to the damage event, while the quadrupolar symmetry of the Eshelby-like redistribution function results in a variation of the sign with the direction with respect to the main loading axis.

Due to these differences, our work captures only qualitatively the mechanisms at play during the progressive damage of quasi-brittle solids. However, it provides rich insights on the effect of the shape of the redistribution function on the quasi-brittle behavior of solids as α( x), and so G( x), can be varied systematically. In particular, our work shows that the resistance to failure decreases with the shape parameter κ. This behavior might be counter-intuitive as large κ means an extended unloaded region after failure events (see Fig. 4). However, this unloading is also accompanied by a stronger reloading of some specific regions that favors the growth of a localized damage mode, and ultimately, the failure of the material along these weakened parts of the material. Indeed, the redistribution function of periodicity λ α = 2π/q α structures the damage field over some material regions where damage accumulates. The change of failure behavior observed at κ = κ c reveals a transition from a reloading mechanism dominated by the positive exponential part of G to another reloading mode dominated by the sign-changing sinusoidal part of G resulting in a structuration of the damage field. As a result, the most unstable mode for κ < κ c corresponds to the homogeneous perturbation q c = 0, so macroscopic failure takes place at the failure point of the individual elements without prior localization regime, as expected for an interaction free material. On the contrary, for κ > κ c , failure takes place earlier because of the collective effects generated by the interactions within the material. Indeed, the sinusoidal load redistribution favors a damage accumulation over some elements only, threatening the integrity of the whole material instead of sharing the applied load over all the elements. Once damage localizes, the load redistribution is more intense on the elements already damaged, resulting in an early overall failure. In our model, the interactions have a clear detrimental incidence on the resistance of the material. This observation reflects a general feature of the failure behavior of heterogeneous materials as the actual resistance of solids is generally lower than the one of their constitutive elements due to the emergence of cracks or localization bands.

Conclusions

In this study, we have explored the path towards failure in damageable heterogeneous solids. Our analysis is based on the description of damage spreading through a continuum medium with a spatially varying resistance to failure. The introduction of a non-local damage variable enables to capture the mechanism of load redistribution following individual damage events. In our description, the transfer of mechanical energy into fracture energy takes place through cascades of damage events resulting from this redistribution process and the presence of material heterogeneities. However, the level of disorder is shown to have a weak influence on the overall failure behavior of the material.

In addition, it can be well captured analytically by exploring the material response close to the homogeneous material solution.

In particular, we have shown that both localization and failure observed during damage growth emerges from the interaction between the constitutive elements of the material. As a result, changing the shape of the redistribution function changes the failure properties of the material. Two distinct behaviors are indeed observed: For mainly positive reloadings, the damage field grows rather homogeneously until abrupt failure takes place. This failure mode is similar to the one of a homogeneous medium without interaction and the critical load coincides with the snap-back instability of the individual elements.

On the contrary, when the redistribution mechanism results in a sufficiently strong unloading of some regions of the material, catastrophic failure is preceded by a stable localization regime where the damage field structures around a particular mode set by the shape of the redistribution function. In that case, the onsets of damage localization and failure strongly depend on the nature of the interaction, and may take place either in the pre-peak or the post-peak regime.

This study sheds light on the failure behavior of quasi-brittle solids, by showing how damage localization and catastrophic failure emerge from the interaction between the constitutive elements of the material. But it can also guide the design of stronger solids throug the control of their collective behavior. The fabrication of increasingly complex architectures opens new perspectives for the control of interactions in solids, and we hope that our approach that bridges material microscopic features to their macroscopic failure properties can be used as a guideline for the design of materials with improved failure properties.

A Derivation of the evolution law of the damage field This Appendix details the calculations of the damage evolution law of Eq. ( 8).

We first show that it can be obtained from the conservation of energy when the loading is varied, ∂E ∂∆ = 0. (A.1)

We then show that it can alternatively be obtained when considering energy conservation for damage variations at constant applied displacement

δE δd = 0 (A.2)
where a fonctional derivative is used for the computations. We start with the derivation with respect to displacement. In Eq. (A.1), we replace the total energy E by its expression given in Eq. ( 4) This condition is now derived through another approach using a functional derivative of the total energy with respect to damage. Let d 0 ( x) be the damage field at a given displacement ∆. We consider a variation δd 0 δ( x -x 0 ) located

∂ ∂∆   Σ 1 2 ∆ 2 k[d( x,

B Expression of the coefficients used for the global energy minimization

The coefficients involved in the expression (46) of the total energy are given by

A(d 0 ) = 1 2 ∆ 2 k(d 0 ) + Y c0 (d 0 + 1 2 ηd 2 0 ) (B.1a) B(d 0 ) = - 1 2 ∆ 2 k (d 0 ) -Y c0 (1 + ηd 0 ) (B.1b) C(d 0 ) = 1 4 ∆ 2 k (d 0 ) 1 + α2 (q c ) 2 + 3 4 Y c0 η (B.1c) D(d 0 ) = - 1 4 ∆ 2 k (d 0 ) 1 3 + α2 (q c ) 2 . (B.1d)
The coefficient T that relates the loading increment δ∆ and the damage increment δd since localization (see Eq. ( 50)) writes as

T = - ∆ 2 h k (d h )k (d h ) 1 + α2 (q c ) 2 + 2Y c0 ηk (d h )(α 2 (q c ) -1) Y c0 η∆ h k (d h ) 1 - 1 α2 (q c ) 1 + 3α 2 (q c ) 2 . (B.2)
Finally, the loading increment at failure follows

δ∆ c = -Y 2 c0 η 2 1 - 1 α2 (q c ) 2 / 2Y c0 η∆ h 1 - 1 α2 (q c ) 1 + α2 (q c ) 2 1 2 T ∆ h k (d h ) + k (d h ) -3∆ 2 h k (d h ) 1 3 + α2 (q c ) 2 T Y c0 η 1 - 1 α2 (q c ) + ∆ h k (d h ) . (B.3)

Figure 1 .

 1 Figure 1. (a) Schematic representation of the loading conditions and the interactions through the material in 1D. The individual elements are represented here as springs.The damage driving force that applies on the central element results from a spatial average of the driving force that applies on the other elements weighted by the interaction function α represented here. This interaction function also interprets as the spatial structure of the load redistribution that takes place in the material after the central heterogeneity is damaged. (b) Spatio-temporal evolution of the damage field for a 1D specimen of N = 2000 elements. Initially, damage is zero everywhere (d = 0 for a driving ∆ = 0, upper line), until catastrophic failure of the specimen for ∆ = ∆ c (d = 1, lower line).

Figure 4 .

 4 Figure 4. Interaction function α for three different values of the shape parameter κ (see Eq. (17)). For practical pruposes, one introduces the parameter κ c = 0.58.

  Figure 5. (a) Normalized Force-Displacement response of a weakly heterogeneous material (σ = 0.001) for = 5 ξ and η = 8 and three different values of κ after averaging over twenty realizations of the disorder, and comparison with the response of the material with local failure behavior (homogeneous material). Catastrophic failure is indicated by dots. The inset is a zoom of the response with κ = 2κ c that shows the deviation from the local behavior close to failure; (b) Methodology used to determine the onset d h of the deviation to the homogeneous response: The typical fluctuations of the damage field characterized by its standard deviation σ d shows a sudden increase indicated by the black dot for the damage level d h . The catastrophic failure is indicated here by the blue dot and correspond to d c .

  21) that we represent as a function of the average damage level d = 1 L Σ d(x)dx in Fig. 5(b) for κ = 3κ c . After a short transient regime, the level of the fluctuations in the damage field remains rather constant before showing a rapid increase from which the onset d h (resp. ∆ h ) is defined: We fit this increase by a linear behavior (green dashed line) and define d h (resp. ∆ h ) from the intersection with the average damage fluctuations (red thick line) far away from failure d d h . The onsets d c and ∆ c of failure are more obvious and they are defined as the average damage level in the material (resp. loading) prior to the final catastrophic failure event. Note that the full failure of one single element is systematically accompanied by the full failure of the whole material, so ∆ c is reached when the first element fully fails. The variations of the onsets d h and d c with κ are shown on Figs. 6(a) and (b) for different disorder levels σ and for an internal length = 5 ξ. The results obtained for = 10 ξ and for a disorder σ = 0.001 are represented with square symbols. The values plotted in Figs. 6 correspond to an average obtained over the different realizations of the disorder, leading to an error bar smaller than the symbol size. For comparison, the damage levels d p and d sb corresponding to the average damage at the force peak and the snap-back instability are represented by horizontal dotted lines. Note that similar graphs showing the variations of the critical loadings ∆ h and ∆ c with the parameters κ and σ can be drawn.

Figure 6 .

 6 Figure 6. Effect of the interaction function on localization and failure: (a) Variations of the localization threshold d h for deviation from the homogeneous material response with the interaction parameter κ; (b) Variations of the failure threshold d c with κ. Dots and squares correspond to internal lengths = 5 ξ and = 10 ξ, respectively, for different disorder levels σ. The position of the snap-back d sb and peak force d p for the interaction free material are indicated by horizontal dashed lines. The solid and dotted lines correspond to the analytical predictions made from the stability analysis of Sec. 5.1 and the energy minimization of Sec. 5.4, respectively.

Figure 7 .

 7 Figure7. Effect of the spatial discretization on the calculated failure response: Comparison of the converged value of the critical loading at failure with the one obtained when varying the number N of heterogeneities in the specimen while keeping N ξ = 1 (main panel) and the number of points N ξ used to discretize a single heterogeneity while keeping the heterogeneity number equal to N = 1000 (inset). The material parameters used for this convergence study are κ = 0, = 2 ξ and σ = 0.001.

  Figure 8. (a) Power spectrum of the damage field plotted as a function of the normalized wavelength λ/ξ = 2π/(qξ) at different distances δ to failure for an interaction function with κ = 2κ c and = 5 ξ.It shows the emergence of a peak for the mode q c . The inset shows a zoom highlighting the far from failure regime δ ≥ 0.55, for which all modes contribute about equally; (b) Variations of this mode q c with the interaction parameter κ for different disorder levels σ and comparison with the analytical predictions (see Sec. 5.1). The mode q α = κ/2 introduced through the interaction function is also shown for both values of .

Figure 10 .

 10 Figure 10. Failure behavior of a 2D specimen: (a) Evolution of the norm q c of the critical wavevector as a function of the shape parameter κ of the interaction function; (b) Variations of the thresholds for heterogeneous damage growth d h (circle) and failure d c (cross) as a function of κ for η = 8 and different disorder levels σ.

Figure 11 .

 11 Figure 11. Energetical interpretation of localization and failure in presence of a stable post-localization regime (κ > κ c ): (a) Energy landscape as a function of the perturbation amplitude δd 0 for different loading increments δ∆ = ∆ -∆ h from the localization threshold. For each loading level, the damage is fixed to the value d 0 = d min 0 determine theoretically to minimizes E ∆ . The position of the minima determined theoretically are indicated by dots; (b) Evolution of the critical displacements ∆ h at localization and ∆ c at failure for η = 8, κ = 2κ c and σ = 0.001.The theoretical values at localization are shown by the solid line, at failure in black dashed for the predictions from the linear stability analysis (referred to as Meth. #1) and red dashed for the predictions from the global energy minimization (referred to as Meth. #2). The horizontal dotted lines indicate the loadings at the peak force and at the snap-back instability.

Figure 12 .

 12 Figure12. Equilibrium states of the material during damage growth in the (d 0 , δd 0 ) plane (see Eq. (32)): (a) For a load redistribution with strong unloading (κ > κ c ), the homogeneous damage state bifurcates continuously to a localized damage stage for d = d h . This stable branch vanishes when it merges with an unstable branch at d = d c , leading to catastrophic failure. (b) For a load redistribution dominated by reloading (κ ≤ κ c ), the stable branch vanishes before any bifuraction to another stable branch took place: The material fails as the damage field is still homogeneous.

  ∆ Y ( x, ∆) -Y c ( x, d) d x = 0 (A.10)Therefore, the equilibrium condition writes asY ( x, ∆) = Y c ( x, d). (A.11)

Catastrophic failure for non-hardening (η = 0) materials takes place right after the elastic regime only for weakly heterogeneous materials. For strong heterogeneities, abrupt failure is delayed and takes place at a higher applied load(Shekhawat et al.,

2013).

This remains only true in absence of heterogeneities. For heterogeneous materials, different elements would break at different loadings, resulting in a progressive increase of the damage level in the material and softening without involving any elastic interactions.

For κ ≤ κ c , we could not identify the localization threshold in the simulaions using the definition given in Fig.5(b), as expected since catastrophic failure occurs without prior deviation to the homogeneous material response

Since the interaction function writes as α( r / ), its Fourier transform writes as α( q), so the normalized quantity q c is naturally used in Fig.10(a) as in Fig. 8(b).

Acknowledgements

The support of the city of Paris and the University Pierre et Marie Curie through the Emergence programs are gratefully acknowledged. We also wish to thank C. Dascalu and D. Kondo for unvaluable discussions.

in x 0 where δd 0 is a small parameter. The damage field hence writes as

which results in a non-local damage field of the form

Expanding the stiffness and fracture energy fields

at first order in δ 0 , we obtain the functional derivative of the elastic energy lim

We obtain similarly the rate of dissipated energy by damage as lim

We obtain from both previous equations the equilibrium condition .17)