A review of non local continuum damage: Modelling of failure?

Failure of quasi-brittle materials such as concrete needs a proper description of strain softening due to progressive micro-cracking and the introduction of an internal length in the constitutive model in order to achieve non zero energy dissipation. This paper reviews the main results obtained with the non local damage model, which has been among the precursors of such models. In most cases up to now, the internal length has been considered as a constant. There is today a consensus that it should not be the case as models possess severe shortcomings such as incorrect averaging near the boundaries of the solid considered and non local transmission across non convex boundaries. An interaction-based model in which the weight function is constructed from the analysis of interaction has been proposed. It avoids empirical descriptions of the evolution of the internal length. This model is also recalled and further documented. Additional results dealing with spalling failure are discussed. Finally, it is pointed out that this model provides an asymptotic description of complete failure, which is consistent with fracture mechanics.


1.
Introduction. Computational modelling and failure analysis of quasi-brittle materials such as concrete is facing the issue of strain softening due to distributed cracking. Upon mechanical loading, these materials undergo progressive microcracking because locally, emerging micro-cracks are being arrested e.g. by heterogeneities or voids. Micro-cracks may thus develop nearby and failure is the result of the propagation of a large quantity of micro-cracks, which, eventually, coalesce in order to form a macro-crack. Fig 1 is an example of such a failure process in a notched three point bending beam, taken from Ref. [18]. The notch is represented by the dark vertical line in the middle of the figure. Due to the tensile stresses at the notch tip, cracking occurs and a macro-crack propagates vertically from the notch towards the top of the beam. In this experiment, the acoustic events generated by micro-crack openings have been recorded and localized on the specimen. Each point on this figure corresponds to a single micro-cracking event. A large number of events has been recorded during loading, indicating that there is a so-called fracture process zone with many micro-cracking events ahead of the macro-crack tip.
Progressive micro-cracking yields a characteristic mechanical response which exhibits strain softening, e.g. a tangent operator in the incremental stress-strain relation which is no longer positive definite. Strain softening causes a loss of well posedness of incremental boundary value problems and, more importantly, it yields prediction of failure without energy dissipation [1]. The governing equations of Figure 1. Acoustic events recorded during the failure of a three point bending beam (reproduced from Ref. [18]). equilibrium loose ellipticity in statics. In dynamics, the incremental equations of motion loose hyperbolicity and waves may no longer propagate.
The remedy to this physically unrealistic feature was found in the early 90's, inspired by the extension of earlier theories of materials with microstructure to cases where the material behaviour is not reversible anymore. Among such models, called enriched models for failure analysis, are the non local damage model [38], gradient plasticity or gradient damage approaches [10,34] and Cosserat models [9]. Their common feature is the incorporation in the constitutive relations of an internal length, which controls the failure process and thus precludes any dissipative process to occur in a region of the solid of zero volume.
It is also interesting to observe that in computational fluid mechanics similar issues were faced with transonic flow and boundary layer problems [27]. The most popular techniques for circumventing those problems were the introduction of a characteristic time, for instance by considering that the fluid is viscous. A small amount of viscosity (damping) provides a regularisation of the governing equations. Similar solutions exist in solid mechanics. Needleman [32] and Sluys [48] demonstrated that material rate dependency is a proper regularisation method, i.e. that it restores well posedness at the inception of softening. Although the technique is very simple, some inconsistencies remain because the response of a quasi-brittle material is slightly more complex compared to simple fluids: rate dependency provides for instance an evolution of the material strength for homogeneously strained specimens at different loading rates. At the same time, it provides also a width of the fracture process zone and eventually controls the fracture energy of the material. The fracture energy and the evolution of the material strength with the loading rate are two characteristics, which should be fitted with a single parameter (the viscosity), which is not always in accordance with experimental data.
Continuum failure models with an internal length have been applied on a wide range of problems involving the safety of concrete and reinforced concrete structures, fracture of sea ice, or the structural strength of composite elements. Recently, it became apparent, however, that such regularized models, e.g. the non local damage model, exhibited some inconsistencies such as (i ) incorrect crack initiation, ahead of the crack tip; (ii ) propagating damage fronts after failure due to non local averaging, (iii ) incorrect shielding effect with non local interactions transmitted across a crack; (iv ) deficiencies at capturing spalling properly in dynamics, with spalls of zero thickness when the expected spall size is below the internal length of the model (see e.g. [47,29,19,26]). Usually, the internal length is considered to be a constant in the formulation and there is a consensus today that some of the shortcomings of the model are due to this assumption. Geers and co-workers were among the first to propose a gradient damage model with evolving internal length [13]. They argued that it is necessary that the internal length decreases in order to model complete failure, i.e. a complete separation of the crack surfaces at failure. This argument was subsequently recovered by Dufour and co-workers among others [19,14,42].
In the meantime, boundary effects are expected to occur in non local formulations, typically, when the neighborhood which influences the response of a material point protrudes the boundary of the structure. This occurs when damage initiates from, or arrive close to a material boundary. Usually, it is treated rather arbitrarily: in integral formulations (as we will see next) the weight function is normalized so that the non local averaging is invariant with respect to a homogeneous distribution of the local variable; in gradient formulations, additional boundary conditions are required for the non local variables (typically, zero gradient conditions normal to the boundary). Considerations from micro-mechanics and size effect analyses indicated, however, that non locality requires a special treatment near boundaries, essentially that the response of the material should become local [29,5,42].
The purpose of this paper is first to review the basic equations of non local continuum damage and the properties of the formulation upon localisation of strain and damage. Then, we examine the properties of a new interaction-based non local model in which the internal length is not constant. The weight function involved in the non local formulation is constructed explicitly from the computation of the elastic interactions between material points. The basic properties of the formulation are reviewed, using simple one-dimensional computations.
This contribution is a review of existing results on non local damage models [38,40,12,39] and on the interaction-based non local model [46]. Subsequently, new results illustrate the shortcomings of the initial non local damage formulation, namely the spurious propagation of damage fronts, and boundary effects, as treated with the interaction-based model, are investigated.
Although the focus of this paper is clearly on the modeling of the mechanical response of concrete, with arguments dealing with the physical aspect of fracture, it is worthwhile to point out that the difficulties encountered in modeling the response of concrete upon strain softening could not have been elucidated without some basic considerations about well-posedness of boundary (or initial) value problems, going back to some fundamental results on partial differential equations. The results presented in section 3 are the outcome of very rich exchanges between mechanics and mathematics which took place in the 90 s, inspired from contributions by Hill [23,24] about uniqueness and stability in elastic-plastic solids and enriched by e.g. the works of Geymonat and co-workers [6].
2. Non local damage model. We shall recall here the simplest version of continuum damage based models, involving a single additional variable, which measures damage as a function of the degradation of the elastic modulus of the material [30]. Although the constitutive relations are discussed in a full 3D format, we will restrict applications to the one-dimensional case, which is more easily tractable and quite simple to implement.
2.1. Isotropic (scalar) local damage model. The classical stress-strain relation for this type of model reads: In Eq 1, σ ij and ε kl are the components of the stress and strain tensors, respectively (i, j, k, l ∈ [1, 3]) and C ijkl are the components of the fourth-order elastic stiffness tensor. The damage variable D represents a measure of material degradation which grows from zero (undamaged material with the virgin stiffness) to one (at complete loss of integrity). The material is isotropic, with E and ν the initial Young's modulus and Poisson's ratio respectively. For the purpose of defining damage growth, a scalar equivalent strain ε eq is introduced, which quantifies the amount of tensile deformation in the material. In this contribution, Mazars' definition of the equivalent strain is used [31]: In Eq 2, ε i + are the positive principal strains. Damage growth is governed by the loading function: g(ε, k) = ε eq (ε) − k .
In Eq 3, k equals the damage threshold ε D0 initially, and during the damage process it is the largest ever reached value of ε eq . The evolution of damage is governed by the Kuhn-Tucker loading-unloading condition: The damage variable D is determined as a linear combination of two damage variables D t and D c , that represent tensile damage and compressive damage respectively, with the help of two coefficients α t and α c which depend on the type of stress state [31]: .
In this paper, we will have only cases where α t = 1 and α c = 0. Standard values of the different model parameters involved in the isotropic local damage model have been given in Ref. [31].

Non local integral formulation.
In the integral-type non local damage models, the local equivalent strain is replaced by its weighted average: In Eq 7, Ω is the volume of the structure and ψ (x, ξ) is the weight function. It is required that the non local operator does not alter the uniform field, which means that the weight function must satisfy the condition: For this reason, the weight function is recast in the following form [38]: In Eq 9, Ω r (x) is a normalization factor and ψ 0 (x, ξ) is the non-normalized weight function, which is often taken as a polynomial bell-shaped function [3], as a Green distribution function [15] or here as a Gauss distribution function: In Eq 10, l c is the internal length of the non local continuum. Preserving the uniform field in the vicinity of the boundary makes the averaging in Eq. (9) not symmetric with respect to its arguments x and ξ. This lack of symmetry leads to the non-symmetry of the material tangent operator [4,41,25].

Associated gradient formulation.
A simple method to transform the above non local model into a gradient model is to expand the equivalent strain into Taylor series truncated for instance after the second order. If we restrict ourselves to an infinite body, we have: ε eq (x, s) = ε eq (x) + ∂ε eq ∂x ds + 1 2 ∂ 2 ε eq ∂x 2 ds 2 + ... .
Substitution in Eq.(7) and integration with respect to variable s yields: In Eq 12, c is a parameter which depends on the type of weight function in Eq. (7). Its dimension is m and it can be regarded as an internal length. Substitution of the new expression of the non local effective strain in the non local damage model presented above yields a gradient damage model. Computationally, this model is still delicate to implement because it requires higher continuity in the interpolation of the displacement field. This difficulty can be solved with a simple trick: let us take the Laplacian of the right and left hand-sides of Eq. (12). Because the Taylor expansion in Eq. (11) was truncated to the second order, higher derivatives can be assumed to be negligible. It follows that: Therefore, Eq.(12) becomes: This implicit equation which defines the non local effective strain as a function of the local effective strain is easy to discretise in a finite element scheme, adding boundary conditions in order to solve it. The implementation of the gradient damage model becomes in fact similar to the implementation of a thermo-mechanical (local) model in which the non local effective strain replaces the nodal temperatures. In the case of an infinite solid, it was demonstrated in Ref. [35] that the gradient equation Eq. (14) is strictly equivalent to the integral formulation, with a specific weight function, i.e. the Green-type weight function which will be used later on in this paper.
3. Failure due to strain localisation.
3.1. Analytical results. In the original paper by Rice [45], strain localization was defined by the bifurcation of the incremental (rate) problem for the linear comparison solid as defined by Hill [24]. In a small strain context, bifurcation occurs whenever strain softening is encountered. Later on, these considerations were given a rational background with the help of basic results on partial differential equations. Upon strain softening, the governing equations have, in statics, an infinite number of solutions and therefore the corresponding boundary value problem is ill-posed. This ill-posedness derives directly from the fact that the tangent (incremental) operator in the constitutive relations is no longer positive definite, thus recovering a result which had been established at the beginning of the 20th century by Hadamard [21,22].
In the case of a non local damage model, and because the non local average introduces a regularization, strain localisation cannot be defined at the onset of a discontinuity of the incremental strain distribution. It was demonstrated in Ref. [40] that such a solution cannot satisfy the governing equations of the problem (constitutive relations and equilibrium). It remains however that strain softening may produce a loss of uniqueness of the solution to a boundary conditions problem. The detection of possible bifurcation points cannot be carried out analytically in the general context. This analysis can be performed in the case of an infinite body only using a linear comparison solid [24] and assuming that the boundary conditions are such that the deformation and damage variables are initially homogeneous over the solid. Hence, we study the conditions of uniqueness and admissibility of small perturbations which satisfy the incremental constitutive relations: In Eq 15, D o and ε o kl are the initial states of damage and strain respectively from which the incremental perturbation occurs, and the momentum equations: In Eq 16, ρ is the mass density and du i is the unknown perturbation (incremental displacement components). Substitution of the equations governing the growth of damage in Eq.(15) provides the following integral form to be inserted in the equations of motion: In Eq 17, h(ε o eq ) is the derivative of the function, which relates damage to the non local equivalent strain computed at the initial state of damage D o (derivative of Eqs.(6) under the assumption of monotonically increasing strain). Substitution of the incremental constitutive relations in Eq.(16) yields an integro-differential system of equations.
Solutions of this system are harmonic waves, propagating in direction n, of amplitude A and phase velocity ζ (for the sake of simplicity of the notation, we leave out the indicial notations in the following): Substitution of these solutions into the equation of motion yields the linear algebraic system: In Eq 19, I is the 2x2 identity matrix and C * o is the tangent operator governing the incremental constitutive equations: In Eq 20, × denotes the tensorial product. n.C * o (ζ). n can be regarded as a pseudo acoustic tensor, function of ζ entering inψ(ζ), which is the Fourier transform of the weight function:ψ The condition of bifurcation, i.e. the condition of existence of harmonic waves is (a homogeneous deformation with constant velocity is already a trivial solution): We may now restrict the analysis to statics (c = 0). Whenever Eq. (22) is satisfied, there exists a non homogeneous incremental solution to the governing equations, in addition to the homogeneous one.
Let us now consider the case of a local continuum, which corresponds to a vanishing internal length. When Eq.(22) is verified, there is an infinite number of solutions, with arbitrary values of ζ since it no longer appears in the equations. This is exactly the criterion of bifurcation introduced by Rice [45]. The incremental problem has an infinite number of solutions and it is ill-posed. When the internal length is non zero, Eq. (22) is verified for specific values of the phase velocity ζ. At the inception of failure, there is in fact a single possible wave number ζ, in addition to the trivial solution [39]. Therefore, the role of the internal length is to control the possible bifurcation modes and also to avoid discontinuous modes (which yield to zero dissipation of energy).
As a remark concluding this paragraph, it is interesting to note that non locality of damage does not completely solve the issue of ill-posedness of the boundary value problem. Strictly speaking, well-posedness can be defined by the existence and uniqueness of the solution, which is not sensitive to a small perturbation of the boundary conditions. Same as in the case of buckling, uniqueness of the solution is lost at the onset of the localization of strain and damage. Non locality prevents zero energy dissipation only. It means that special numerical schemes are necessary if one wants to investigate the possible equilibrium paths upon bifurcation [44].

3.2.
Dynamic failure of a rod. Let us now illustrate the response of the non local damage model upon localisation on a simple one-dimensional case [38]. A bar of length L, in which two constant strain waves converge toward its center is considered (see Fig. 2). The amplitude of the wave is 0.7 times the deformation at the peak load in tension. When the two waves meet at the center, the strain amplitude is doubled, the material enters the softening regime suddenly and failure occurs.
To study damage front propagation and to avoid wave reflections on the boundaries, a long bar is considered (L = 300 cm). The parameters used in this example are: mass density ρ = 1 kg/m 3 , Young's modulus E = 1 MPa and a velocity boundary condition of v = 0.7 cm/s applied at the bar ends. The other model parameters are A t = 1, B t = 2, ε D0 = 1 and the internal length l c is 4 cm (there is no damage in Figure 2. Dynamic failure of a rod: test description and time evolution of the strain amplitude repartition along the rod (reproduced from [19]). compression). A fixed mesh of 999 constant size elements is used. Time integration is performed according to an explicit, central difference scheme. The time step is chosen to be equal to the critical time step of the elements (∆t = 1 s). The two waves meet at the center of the bar (x = 150 cm) at time t = 150 s. Fig. 3 shows the evolution of damage over the bar at different time steps according to the original non local damage model with both Gauss-type or Green-type weight functions (Eq 23).
In Eq 23, R c is chosen in such a way that the integral of both weight functions on a infinite domain are equal (R c = l c √ π 4 ). For both cases, the damage band has a finite width. A singular feature observed in Fig. 3 is that when damage is complete in the middle of the bar (damage= 1), the width of the damage band increases. In other words, there is a propagation of the damage front, even though the bar does not experience any stress in the vicinity of this front. This is due to the non local averaging. At material points located just ahead of the damage front, in the region which is not damaged, the non local equivalent strain grows because averaging encompasses the strain at neighbors located inside the damage band. At these neighbor points, the strain is very high and therefore, the non local equivalent strain just ahead of the front reaches the damage threshold eventually. The propagation of the damage front after failure is obviously a spurious consequence of averaging which ought to be circumvented. To this end, let us look at the physical motivations for introducing non local damage.
It is commonly accepted that non locality is the result of the interactions between material points undergoing damage in the course of failure [2,37]. Several interaction mechanisms may be distinguished, which provide some requirements for the consistency of the constitutive relations: (i ) an interaction exists if there is damage and it grows with damage [43]; (ii ) shielding occurs when two points are located apart from a crack (or a non convex boundary), they do not interact whatever their distance; (iii ) on free existing or evolving boundaries, and along the normal to these boundaries, non local interactions should vanish as demonstrated in [42].
The weight function which is implemented in the non local averaging should yield results which capture at the continuum level these mechanisms because it is this function which carries the non local interactions from one point to its neighbor (see for instance the comparisons in [19]). In most cases, this weight function is defined a priori, possibly with some variations of the internal length [14]. In the following Damage: D(ε(X, t)) in Eq 6,ε from Eq 7. section of this paper, we discuss a model in which the weight function is not defined explicitly. It is upscaled from the description of elastic interactions and therefore evolves during the failure process.

4.
Interaction-based non local weight function. The aim is here to estimate the non local weight function ψ presented in Eq. (7) directly from elastic interactions between material points. Let us consider two material points (x, ξ) of a quasi-brittle structure and an infinitesimal strain perturbation locally produced at point ξ. Considering the interaction between x and ξ consists in estimating the non local contribution of the Figure 4. Non local contribution seen by a point x when a perturbation is produced in ξ (reproduced from [46]).
strain perturbation of ξ on x. This problem can be written as: In Eq. (24), Ω is the volume of the structure, δ is the Dirac distribution, ε * is a given strain field such as ε * ξ = δε * represents the local strain perturbation centered in ξ, ψ int x is the interaction-based non local weight function to be constructed centered in x, and ε ξ (x) represents the non local contribution of the strain perturbation of ξ on x.
Following Eq. (25), Eq. (24) can be rewritten as: Therefore, the interaction-based non local weight function ψ int x , centered on x, can be simply reconstructed by considering on each point ξ the measure of the non local contribution seen by a point x when a perturbation is produced in ξ.
Practically, the interaction-based non local weight function ψ int x is reconstructed numerically. The strain perturbation is produced by a thermal expansion: ε * = α∆T I. This thermal expansion is imposed on each integration points contained in a sphere centered in ξ (see Figure 4). The radius a of the sphere is a model parameter and is related to the size of the fracture process zone, i.e. to the maximum aggregate size. The measure of the non local contribution is chosen (arbitrarily) equal to the Euclidean norm of the eigenstrains generated by the thermal expansion. The interaction-based weight function writes: In Eq. (27), ξ i (x) represents the value at point x of the ith-eigenstrain when a thermal expansion is imposed at point ξ.
Finally, the classical non local integral formulation is modified by taking into account the new integration-based weight function presented in Eq. (27). Similarly to Eq. (7), a normalisation is introduced in such a way that the non local operator is invariant with respect to a uniform field of equivalent strain at the onset of damage in the structure. The non local formulation writes: In Eq. 28, ψ int x is the new interaction-based weight function defined in Eq. 27 and ψ 0 x represents the same function reconstructed when no damage occurs in the structure (typically at the beginning of the computation).
Since the reconstruction process of the interaction-based weight function is kinematically driven by the successive thermal expansion imposed at each point of the structure, we choose that all boundaries are clamped in order to avoid any perturbation of the kinematics on the boundary.

Clamped bar in tension.
We consider here the simple problem of a tensile bar (see Fig. 5) in order to look at the properties of the interaction-based weight function. The response of the material is assumed to be purely elastic excepted in a damaged area located in the center of the bar. The structure is modeled as a succession of circular inclusions (radius a) which will be dilated successively in order to reconstruct the interaction-based non local weight function. Note that although this is a 1D problem, the calculation of the interaction-based weight function follows a 2D description, with a line of circular inclusions of radius a located along the bar, centered in the middle of its cross section. Figure 5. Simple problem of a clamped bar in tension (L = 10cm, σ Y = 3M P a, E = 30GP a, ε D0 = 10 −4 , σ 0 = 1.8M P a, D = 70%)(reproduced from [46]).
At the left bar end, we have placed an inclusion in order to check the evolution of the weight function close to the boundary. Fig. 6 shows the various weight functions obtained for a decreasing radius of the inclusion. The weight function tends to a Dirac distribution as the size of the inclusion decreases. It spans over less and less elements and the amplitude is increasing. Because the construction of the weight function results from a succession of Eshelby-like problems (inclusions inserted in  an elastic matrix), this conclusion could have been inferred directly. The difficulty here is that we deal with a finite size solid and analytical results are not available. We may conclude that if the non local interactions should vanish near the end of the bar as derived in [42], the size of the inclusions considered for the calculation of the weight function should decrease when material points are getting closer to the boundary and should vanish at the boundary. As we will see further, a line where damage is equal to one is also a boundary and therefore, the size of the inclusion, which is the counterpart of the internal length entered in the original non local formulation, should decrease with increasing damage.
When the inclusion stays inside the undamaged part but gets close to a damaged area, Fig. 7 shows that a shielding effect is observed whereas the non local contribution are higher inside the damaged area with the original formulation. The only restrictive condition is that the inclusion size should stay smaller than the damaged zone, otherwise some information may be transferred through the inclusion itself.
These first observations on a very simple example of a clamped bar in tension give some insight about the evolution of the size of the sphere containing the integration points where the thermal expansion is imposed in order to construct the interactionbased weigh function. This sphere has to decrease in size when it gets close to a boundary in order to recover a local behavior. Moreover, this sphere has to decrease in size when damage increases to ensure an efficient shielding effect.
To meet these requirements, we have chosen the following empirical rules: In Eq. 30, a(x) represents the radius of the sphere containing the integration points where the thermal expansion is imposed to construct the interaction-based weight   function, a 0 is a model parameter related to the maximum aggregate size, D is the local damage, d is the minimal distance from any boundary of the structure to point x.

4.2.
Dynamic failure of a rod. The simple one-dimensional case presented in section 3.2 is now considered with the interaction-based model. The main case characteristics are collected in Table 1. Fig. 8 presents the evolution of the damage profile with time over the bar according to the interaction-based model. The damage band has still a finite width but there is no longer propagation of the damage front, which is provided by the shielding effect appearing when the middle element is fully damaged. Damage: D(ε(X, t)) in Eq 6,ε from Eq 28. This shielding effect is highlighted by Fig. 9 representing the normalized weight function estimated close to the middle fully damaged element (at x = 14 cm for a middle element located ar x = 15 cm). Whereas the original non local weight functions (Gauss-type and Green-type) are not influenced by the evolving damage in the structure, the new interaction-based non local weight function encompasses intrinsically a shielding effect close to a damage area. An animated version of Fig. 9 is proposed as supplementary material 1 showing the evolution upon damage of the new interaction-based weight function. Fig. 10 represents the local strain, the non local strain, the damage and the non local contribution estimated at the middle of the bar at the end of the computation when the middle element is fully damaged. After failure, the middle of the bar should act as a new boundary and non local interactions should vanish and a local behavior should be recovered. Is is observed only with the new interaction-based formulation, where the non local strain tends to zero at total failure. An animated version of Fig. 10 is proposed as supplementary material 2 showing the evolution upon damage of the local strain, the non local strain, the damage and the non local contribution estimated at the middle of the bar.

4.3.
Spalling test. A second one dimensional example is used to test the response of the model close to a boundary. This is a simplified description of a spalling test based on a split Hopkinson pressure bar test primarily developed by [28] for material dynamic behavior characterization, but often adapted for dynamic fracture testing [17,16]. A striker bar generates a square compressive wave that propagates along the bar in the linear elastic regime. When this compressive wave reaches the free   extremity of the bar, it is converted into a tensile wave and added to the incoming compressive wave (see Fig. 11 and Table 2). The resulting wave stays equal to zero until the tensile one reaches a distance from the boundary equal to half the initial signal length. Failure is initiated at this point if the absolute value of the amplitude of the wave is greater than the damage threshold, generating a spall at a controlled distance from the boundary that depends on the initial compressive signal duration. For all numerical studies, the time step is chosen to be equal to the critical time step of the corresponding element size. Fig. 12 shows that the spalling failure is better described with the interactionbased model. When the spall initiates far from the boundary (Fig. 12.a), all weight functions (Gauss, Green and interaction-based) provide relevant results. When the spall initiation goes closer to the right boundary, the final spall location is correctly predicted inside the bar with the interaction-based nonlocal model, whereas damage is maximum on the boundary with the original model (Fig. 12.c and Fig. 12.d). The original non local damage model cannot predict a spall size less than the order of the   Figure 11. Spalling test: test description and time evolution of the strain amplitude repartition along the rod (reproduced from [19]).  5. Continuum to discontinuum description of failure. In the course of damage, failure is captured by continuous functions. At the location of a macro-crack, the displacements are continuous although they should be discontinuous if the description of complete failure is expected. This is clearly a shortcoming of the non local damage approach, which sort of embeds continuity in a process which is expected to become discontinuous, asymptotically at least. In applications, it is impossible to obtain directly a crack opening displacement (COD) which is defined as the displacement jump across a crack. This quantity is important in order to evaluate the durability of concrete structures. Fluid flow inside cracks, or diffusion controlled processes such as carbonation are strongly related to the amount of crack opening for instance and it is important that the model can connect a distribution of damage to a crack opening.
To be consistent with fracture mechanics, at least in an asymptotic sense when damage has reached 1, the strain distribution at complete failure provided by the damage model should tend towards a Dirac distribution, and numerically it should be as close as possible of it. If it is not, it means that the damage model cannot converge toward a full fracture formation process. 5.1. Numerical considerations. The crack opening displacement can be estimated using the method proposed by [12,11]. The idea is that instead of comparing the strain distribution to a Dirac distribution, it is easier to compare the distribution of the non local equivalent strain to that resulting from a strong discontinuity (discontinuity of the displacement field) (see Eq. 31). The comparison provides two results. First, the minimisation of the distance between the two distributions yields a crack opening displacement, that is the best value of the displacement jump according to the non local damage model. Second, the difference between the two distributions indicates how close (or far) the two non local equivalent strain distributions are. This is summarized in the following expressions: [U ](x) is the COD computed at the point of coordinate x that corresponds to the location of the crack, ∆(x) represents the distance between the computed distribution of non local equivalent strain and the ideal one obtained for a strong discontinuity (single crack at point of coordinate x), ε eq is the non local strain and ε sd is a non local strong discontinuity-based strain profile resulting from a crack located at coordinate x, e.g., the non local equivalent strain corresponding to a local strain described by a Dirac distribution. Practically, the Dirac distribution is regularised by using the classical Gauss-type distribution function (Eq. 10) but the results do not depend on the choice of a particular weight function. Details may be found in [12]. Fig. 13 shows the calculation of the distance between the obtained distribution of the non local equivalent strain and the distribution which results from a strong discontinuity computed during the dynamic failure of a bar investigated in section 4.2. The original non local damage model and the interaction-based non local damage model have been considered. The failure process is well described with the interaction-based model since the distance ∆(x) tends rapidly to zero. This is not the case with the original model with a gaussian weight function and a constant internal length. This result is also a consequence of the observed propagation of damage fronts. The original non local model does not converge towards an asymptotic description of cracking, at least according to the above definition provided by the calculation of the COD. At complete failure, the crack opening computed according to the above technique (see Eq. 31) should be independent of the element size. In a simple 1D setting, and assuming that the crack opening is smeared over the finite element that contains the discontinuous displacement at complete failure, the crack opening is also equal to the strain distribution times the element size. Therefore, after complete failure, the strain ε in the cracked element should evolve in inverse proportion of the element size h (for constant strain elements):  Figure 13. Dynamic failure of a rod: distance between the computed COD and an ideal opening profile obtained for a strong discontinuity versus time.
Fig. 14 shows that the complete failure is better described with the interactionbased model since the strain versus normalized element size curve follows a linear trend in a logarithmic plot. Moreover the slope is consistent with the numerical value of the COD estimated at complete failure. Note that in the original model, the element size is dimensioned by the internal length l c whereas in the new model, it is dimensioned by the characteristic length a 0 . For the integration-based model, a discontinuity is observed when the element size is approximately equal to the characteristics length a 0 in Fig. 14. It means that several elements are needed inside the inclusion where the perturbation is applied in order to construct the interaction-based weight function accurately.
The above considerations are based essentially on numerical analyses and a heuristic approach rooted in simple physical principles. From a theoretical point of view, however, the reasons why the original non local damage model fails at describing complete failure and why the interaction-based model seems to perform properly are not elucidated yet. Convergence towards a discontinuity upon complete damage, and the main ingredients at stake for a consistent result are yet to be clarified and identified. In a different setting, the numerical implementation of the variational approach of fracture (e.g. proposed by Bourdin, Marigo and Francfort and co-workers [8,7]), can be based on a regularization of the problem (which can be understood as the introduction of an internal length). Γ-convergence considerations ensure that the variational model converges toward a physically consistent description of fracture. Several other results (e.g. by Negri [33]) relies on the same sound framework for the non local approximation of a free discontinuity. Recently, the variational approach implemented for fracture has been extended to damage and gradient-damage models [36]. The step towards a complete demonstration of convergence towards a complete displacement discontinuity upon failure is still missing. In particular, it is important to determine which properties the weight functions ought to exhibit in order to achieve this convergence and whether there is an effect of the damage evolution law. Most probably, the mathematical setting provided by the variational approach lends itself better for this demonstration compared to the integral damage model discussed in this paper. 6. Conclusions. Because of strain softening due to progressive cracking, the mechanical response of quasi-brittle materials such as concrete needs to incorporate an internal length. The constitutive relations are non local and the internal length that enters in the model controls the process of damage and prevents bifurcation to (discontinuous) zero energy dissipation modes. The internal length controls the size of the damage band upon inception of bifurcation. In fact it controls also the size of the fracture process zone ahead of a propagating macro-crack in a more general context. Integral and gradient-type damage models have been proposed in the literature for this purpose. The early versions of the non local damage model relied on a constant internal length. It turns out, however, that this simple description has some serious shortcomings. Propagation of damage fronts upon failure is observed, in addition to other features such has incorrect initiation of damage from a notch as mentioned in the literature. This is the reason why some authors started to devise non local models with a varying internal length.
Non locality is motivated by the interactions between the distributed microcracks which develop in the process of failure, and between these micro-cracks and the boundary of the solid. We have reviewed a technique used to construct the weight functions which carry the non local interactions in the course of damage. This yields a new interaction-based non local formulation with the following features: • The scheme provides a shielding effect in the calculation of the non local average which avoids the propagation of a damage front. The computed weight functions evolve so that the effect of material points lying in a damage zone on points close to this zone vanish with growing damage. • Near boundaries, interactions are also decreasing, which is consistent with considerations from micro mechanics. • In the course of damage, the failure process is better described. The crack opening profile is very close to an ideal opening profile obtained for a strong discontinuity. Unfortunately, there is no theoretical proof at this time that a convergence towards a discontinuity upon complete damage is indeed reached.
Understanding why the interaction model provides a better result remains an open problem. • At complete failure, the crack opening is independent of the element size but the strain in the fully damaged element evolves in inverse proportion of the element size. • Spalling failure is better described with the interaction-based model. An arbitrary small spall thickness can be obtained whereas it is not possible to reach a spall thickness less than the internal length with the original non local damage model.
Finally, let us stress that the present considerations have been restricted to onedimensional cases even though interactions are computed according to a 2D description. The extension to 3D cases is possible but remains to be performed. Then it will be possible to tackle more challenging issues of quasi-brittle materials failure such as size effect on structural strength and the description of load-displacement responses of notched and unnotched geometrically similar beams [20].