Duc Trung Le

Jean-Jacques Marigo

Corrado Maurini

Stefano Vidoli

Strain-gradient vs damage-gradient regularizations of softening damage models

Local damage models with softening needs localization limiters to preserve the mathematical and physical consistency. In this paper we compare the properties of strain-gradient and damage-gradient regularizations. Gradient-damage models introduce a quadratic dependency of the dissipated energy on the gradient of the damage field and are nowadays extensively used as phase-field approximation of brittle fracture. Their key feature is to provide a smeared approximation of a crack as a band of localised damage with a finite energy dissipation per unit of surface, that can be identified with the fracture toughness of the Griffith model. Strain gradient models introduce a quadratic dependence of the elastic energy on the gradient of the strain field. A similar term can be physically interpreted as the presence in the material of linear, but nonlocal, stiffnesses, that can be eventually be affected by damage. Despite this attractive interpretation, we have found that strain-gradient regularized models can hardly be used to approximate brittle fracture, because smeared cracks with non-vanishing and finite dissipated energies are hardly obtained. Our analysis is based on variational models and focuses on the one-dimensional traction problem.

Introduction

Damage models describe material failure by an additional internal variable modulating the elastic stiffness and inducing an internal energy dissipation in the material. In order to model material failure observed in quasi-brittle materials, damage models must include stresssoftening, i.e. the reduction of the admissible stress domain for increasing damage. It can be shown that this constitutive property is associated to the loss of the uniqueness of the solution and the appearance of localized states, a key aspect of the behavior of real structures. It is widely recognized that in order to preserve the mathematical consistency and the capability of formulating mesh-independent numerical approximation, strain-softening damage models must be regularized by introducing some sort of non-local effects, [START_REF] Bazant | Continuum theory for strain-softening[END_REF]. Among the several regularization strategies proposed in the literature, we may first distinguish between those based on the introduction of smoothed variables in the form of convolution integrals of the local variables (nonlocal damage models, see e.g. [START_REF] Comi | A non-local model with tension and compression damage mechanisms[END_REF]) and those penalizing extreme localisations through gradient terms (gradient damage models). One can further classify the possible approches between those introducing the regularization through nonlocal terms in the damage variable (damage-regularized damage models, see e.g. [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF]) or in the elastic strain variable (strain-regularized damage models, see e.g. [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF][START_REF] Peerlings | Localisation issues in local and nonlocal continuum approaches to fracture[END_REF][START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF][START_REF] Peerlings | A thermodynamically motivated implicit gradient damage framework and its application to brick masonry cracking[END_REF]). In [START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF] the authors perform a throughout analysis of the possible nonlocal regularization strategies. They conclude that only a very small subset of the models proposed in the literature are really effective in providing mathematically well-posed and thermodynamically consistent problems. They suggest the use of gradient approches. In the present paper, we will focus on gradient-type regularizations and compare the properties of strain-gradient (SG) and damage-gradient (DG) models, discussing their appropriateness to model fracture of brittle solids.

Failure often manifests in the form of cracks, i.e. surfaces of discontinuity of the displacement field, where the material "breaks". Damage models are often regarded as smeared crack models smoothing out this discontinuity in bands of finite length. They are opposed to discrete approaches, which models the sharp discontinuities and the possible jumps of the displacement field explicitly. The most widely accepted discrete crack model is the Griffith model, which assumes that the creation of the crack is done at expenses of an energy dissipation proportional to the crack surface, and that the crack propagation is not possible if this energy is larger than the elastic energy rewarded during a virtual crack propagation. In the end of the nineties, this energetic theory has been put in a precise mathematical setting by [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], allowing for the generalization to the case of cracks of arbitrary shapes and with complex evolution in time.

The key advantage of the Griffith model with respect to the damage model is its simplicity when the crack path is postulated in advance. Vice-versa, damage models are much more practical when considering cracks of unknown and possibly complex shapes, because they do not require the explicit description of the crack geometry: the damage field is treated as a standard field modulating the stiffness and the cracks are identified a posteriori as the regions where the elastic stiffness vanishes. Another fundamental advantage of damage models is to be able to retrieve crack nucleation, i.e. the creation of a crack from a intact material with smooth boundaries.

The so-called phase-field models of fracture are gaining an increasing popularity because they combine the advantages of the smeared and the discrete approaches and give a precise meaning to the idea of using damage models to approximate sharp cracks, or vice-versa1 . These models have been independently developed in different contexts. In applied mathematics and theoretical mechanics, they arise through the transposition [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] of the regularized models for image segmentation [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF] to the variational formulation of fracture mechanics [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] . In physics they derive from the application to fracture [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF][START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF] of the Ginzburg-Landau theories of phase transformations. Phase-field models of fracture are also a special type of the DG models presented in [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF].

Previous studies [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] analyzed in depth the properties of DG models, showing that they can be regarded as a regularized version of the variational theory of brittle fracture, where the regularization parameter, the internal length, can be set to recover crack nucleation in agreement with experimental observations [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. In this paper we analyse whether a strain-gradient (SG) model can be similarly used to approximate brittle fracture à la Griffith. In particular, is it possible to recover with SG models the energetic equivalence with brittle fracture?

Our analysis relies on a variational approach and focuses only on a one-dimensional traction problem. Specifically we consider strain-gradient damage model where the elastic energy density is quadratic in both the strain and the strain-gradient. The corresponding local E(α) and nonlocal G(α) stiffnesses are assumed to depend on the damage variable α with the local stiffness vanishing for α = 1. In a one dimensional traction problem, we identify cracks as solutions with vanishing stress and a smeared displacement jump. We show that two fundamentally different qualitative behaviors are possible, depending on whether the limit value of the nonlocal stiffness G(α = 1) is vanishing or not:

1. when G(1) = 0 the regularizing term is not effective and the SG-model shows pathologies similar to the non regularized case, e.g. the creation of cracks with vanishing energy dissipation, resulting in mesh-dependence of numerical approximations;

2. when G(1) > 0 cracks, the creation of cracks requires to fully damage the whole bar and, hence a dissipation proportional to the bar length.

In both cases, and differently from DG models, it is not possible to recover an energetic equivalence with the Griffith model where the energy dissipation to produce a crack is a finite material constant G c , the fracture toughness. Our theoretical results are supported by numerical illustrations. The reader is referred to Section 7 for a more detailed summary of our findings.

The paper is organized as follows. In Section 2 we introduce the variational setting of the local damage model and discuss why regularizing terms are needed. Section 3 recalls the features of the DG models, whilst the SG regularization are introduced in Section 4. The core of the paper is contained in Sections 5 and 6 where both families of SG models are studied in details. Section 7 is devoted to draw our conclusions.

For the sake of reproducibility, open-source implementations of the strain-gradient damage models used to obtain the numerical results of this paper are provided as supplementary material.2

Local damage models with strain softening

We introduce here the basic ingredients of the local model at the basis of the regularizations discussed in the rest of the paper. We illustrate the response of the material element of a bar under uniaxial traction in a one-dimensional setting, denoting by σ, ε, and α the scalar stress, strain, and damage field, respectively. We consider the strain as a loading parameter, and assume that its parametrisation with the time, t → ε t , to be monotonically increasing. When there is no risk of confusion, the explicit reference to the the variable time t is omitted. Superimposed dots will denote derivatives with respect to t, α = dα t /dt, whilst a prime will denote a derivative of a function with respect to its argument, e.g. E (α) = dE(α)/dα.

Formulation

We follow the presentation proposed by [START_REF] Marigo | Formulation d'une loi d'endommagement d'un matériau élastique[END_REF] in the spirit of Generalized Standard Materials [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF]. Specifically, we consider a damage model characterized by the total energy density

W 0 (ε, α) = 1 2 E(α)ε 2 + w(α), (1)
sum of the elastic energy and the dissipated energy. Here ε denotes the strain and α denotes the scalar damage variable which is chosen in such a manner that it grows from 0 to 1, α = 0 corresponding to the undamaged state and α = 1 to the completely damaged state. The smooth monotonically decreasing state function α → E(α) gives the evolution of the Young modulus of the material with its damage state. We assume that

E(0) = E 0 > 0, E (α) < 0, ∀α ∈ [0, 1), E(1) = 0. (2)
Accordingly, we adopt the following hypotheses in terms of the compliance state function α → S(α) = 1/E(α):

S(0) = 1/E 0 > 0, S (α) > 0, ∀α ∈ [0, 1), S(1) = +∞. ( 3
)
As we will see below, the smooth monotonically increasing state function α → w(α) gives the evolution of the dissipated energy density by the material with its damage state. We assume that w(0

) = 0, w (α) > 0, ∀α ∈ [0, 1), w(1) = w 1 < +∞. ( 4
)
The damage evolution problem is characterized by the following items Stress-strain relation :

σ = ∂W 0 /∂ε = E(α)ε, (5)
Irreversibility condition : 0 ≤ α ≤ 1, α ≥ 0, (6)
Damage yield criterion :

∂W 0 /∂α = 1 2 E (α)ε 2 + w (α) ≥ 0, (7)
Consistency equation :

(∂W 0 /∂α) α = 0. ( 8
)
The damage law is standard in the sense that the damage yield criterion is stated in terms of the elastic energy release rate - 1 2 E (α)ε 2 . This type of criterion can be justified in a full three-dimensional setting by invoking Drucker-Ilyushin postulate, see [START_REF] Marigo | Constitutive relations in plasticity, damage and fracture mechanics based on a work property[END_REF]. In our uniaxial context, [START_REF] Comi | A non-local model with tension and compression damage mechanisms[END_REF] can also read as

|ε| ≤ ε D (α) := 2w (α) |E (α)| . ( 9
)
Accordingly, the set of admissible strain states is an interval which depends on the damage state. We will assume that this interval grows when the damage grows. This leads to the following Strain hardening condition :

α → w (α) |E (α)| is monotonically increasing. ( 10
)
This condition will allow us to obtain a unique response in uniaxial test under controlled strain path.

The damage yield criterion [START_REF] Comi | A non-local model with tension and compression damage mechanisms[END_REF] can be expressed in terms of the stress σ and read as

|σ| ≤ σ D (α) := 2w (α) S (α) . ( 11
)
Since we are only interested to softening behaviors, i.e. the case when σ D (α) is a monotonically decreasing function of α, we adopt the following Stress softening condition : α → w (α) S (α) is monotonically decreasing. [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF] Note that, by virtue of (3)-( 4), σ D (1) = 0 and hence the material cannot sustain any stress when it is completely damaged.

Response of the volume element

We examine the response of a volume element submitted to a monotonically increasing strain ε = t, using Eqns. ( 5)- [START_REF] Forest | Nonlinear regularization operators as derived from the micromorphic approach to gradient elasticity, viscoplasticity and damage[END_REF]. This solution maybe also be interpreted as the homogenous response of a bar of length L submitted to end displacements u(±L/2) = ±ε L/2.

The stress-strain relation reads as

σ = E(α) ε, ( 13
)
and the damage criterion gives:

α = α H (ε) :=          0, |ε| ≤ ε C w -E -1 ε 2 2 , ε C < |ε| < ε L 1, |ε| ≥ ε L , ( 14
)
where the function (-w /(E )) -1 denotes the inverse function of α → -w (α)/E (α), which is well-defined by virtue of the strain hardening condition [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. The material constants ε C and ε L are given by

ε C = ε D (0), ε L = lim α→1 ε D (α), (15)
with ε D defined in [START_REF] Francfort | The variational approach to fracture[END_REF]. These correspond, respectively, to the elastic limit strain and to the strain required to reach the fully damaged state; ε L may be finite or not depending on the model. Substituting ( 14) in [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], one gets the stress in the bar and hence the stress-strain relation for the homogenous response:

σ = σ H (ε) = E(α H (ε)) ε. ( 16
)
The total energy density for the homogeneous solution is given by

W H (ε) = E(α H (ε)) ε 2 2 + w(α H (ε)). ( 17
)
Because of the damage criterion, the derivative of the energy density with respect to the strain ε gives the stress-strain response:

W H (ε) = E (α H (ε)) ε 2 2 + w (α H (ε)) α H (ε) + E(α H (ε)) ε = σ H (ε). ( 18
)
Hence we can identify the strain energy with the area of the region bounded by the graph of the stress-strain response σ H (ε). The stress is limited by an elastic limit stress

σ C σ H (ε) ≤ σ C := E 0 ε C = 2w (0) S (0) , ( 19
)
the response being linear elastic with the stress increasing from 0 to σ C for ε ≤ ε C and with softening for ε > ε C where the stress decreases monotonically from σ C . Indeed:

σ H (ε) = W H (ε) = σ D (α)α H (ε) < 0, (20)
because of the stress softening [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF] and strain hardening (10) assumptions. We require as a constitutive assumption that the elastic energy in the fully damaged state is vanishing:

lim ε→∞ E(α H (ε)) ε 2 2 = lim α→1 E(α)w (α) -E (α) = 0, (21)
which is verified in particular if

ε L = ε D (1) < ∞ or if E (1) = 0 and ε D (1) < ∞.
With this assumption, a fortiori σ ∞ = lim ε→∞ σ H (ε) = 0, and we get that the limit value of the total energy for large strains is:

W ∞ H = lim ε→∞ W H (ε) = lim α→1 W H (ε D (α)) = w 1 , ( 22
)
which justifies the definition of w 1 as the specific fracture energy (see also [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]).

For the following developments, it is useful to explicitly distinguish between the two qualitatively different behaviours obtained depending if the deformation ε L = ε D (1) required to reach the fully damaged state is finite or not. Examples of the two different behaviours are given in Figure 1, which reports the homogenous stress-strain responses for the two models

(LS) : E(α) = 1 -α 1 + (k -1)α , w(α) = w 1 α, ( 23
) (NS) : E(α) = (1 -α) 2 , w(α) = w 1 α. ( 24
)
These are characterised by a linear elastic phase followed by a Linear Softening (LS) or a Nonlinear Softening (NS) phase, respectively. For w(α) = w 1 α, the strain-hardening and stress-softening condition simplifies into:

E (α) > 0, E (α) 2 -E (α)E(α) > 0, ∀α ∈ [0, 1), (25)
and we get ε L → ∞ if and only if E (1) = 0. In this case, the condition of vanishing elastic energy ( 21) is verified if E (1) < 0 or if E (1) = 0 and E (1) > 0.

Need for a regularization

Let us consider a traction problem on a one-dimensional bar of length L submitted to monotonically increasing end displacements. The quasi-static evolution of damage can be formulated as a minimality condition for the energy functional

E(u, α) = L/2 -L/2 W 0 (u (x), α(x)) dx. ( 26
)
where u denotes the displacement and ε(x) = u (x).

In a time-discrete setting3 , let ±t i L/2 be the displacements imposed on the ends ±L/2 of the bar at the time step i ≥ 0 and let (u i , α i ) be the current state of the bar. We assume that, at time t 0 = 0, the bar is undamaged and stress free, and hence that u 0 (x) = 0 and α 0 (x) = 0 everywhere in the bar. Let us denote respectively by C i and D i the sets of admissible displacements, u(±L/2) = ±t i L/2, and damage fields 1 ≥ α ≥ α i-1 ≥ 0. The current displacement and damage field (u i , α i ) are the elements of C i × D i which are local minimizers of the energy E in the following sense 4 :

∃h i > 0, ∀(u, α) ∈ C i × D i such that (u, α) -(u i , α i ) ≤ h i , E(u i , α i ) ≤ E(u, α). ( 27
)
If (u i , α i ) is a local minimum in the sense of [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] such that α i (x) < 1 everywhere in the bar, then (u i , α i ) must verify the following first order optimality conditions

E (u i , α i )(u -u i , 0) = 0, ∀u ∈ C i , E (u i , α i )(0, α -α i ) ≥ 0, ∀α ∈ D i , ( 28
)
where E (u i , α i )(v, β) denotes the directional derivative of the energy functional E at (u i , α i ) in the direction (v, β). In particular,

E (u i , α i )(v, 0) = L/2 -L/2 ∂W 0 ∂ε v (x) dx, E (u i , α i )(0, β) = L/2 -L/2 ∂W 0 ∂α β(x) dx. ( 29
)
By standard arguments of calculus of variations, one can show that (28) 1 implies the equilibrium equation σ = 0 with σ defined in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. The variational inequality (28) 2 yields as Kuhn-Tucker conditions the damage yield criterion [START_REF] Comi | A non-local model with tension and compression damage mechanisms[END_REF] and the time-discrete versions of the consistency equation ( 8) and irreversibility condition [START_REF] Braides | Approximation of free-discontinuity problems[END_REF], where α (α i -α i-1 )/(t i -t i-1 ).

The following proposition, proven in AppendixA, demonstrates that such a local damage model is mathematically ill-posed and it is not satisfactory from a physical point of view. Proposition 1. For the local model [START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF], every state (u i , α i ) satisfying (28) 2 as equality on a set with non-vanishing measure is unstable.

This fact implies a pathological mesh-dependence when trying to numerically minimize the functional [START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF] and call for its regularization, see e.g. [START_REF] Bazant | Continuum theory for strain-softening[END_REF][START_REF] Bazant | Nonlocal Continuum Damage, Localization Instability and Convergence[END_REF]. Indeed, when using finite-element discretizations, the deformations and damage completely localize in bands whose width depends on the mesh. Hence, the dissipated energy goes to zero when refining the mesh size. In 2D and 3D problems mesh dependency implies also the inability to predict the crack shape.

We illustrate this pathology, reporting the finite-element results of a traction test on an bar based on the local model [START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF]. The bar is clamped at the left-end, it has a varying elastic stiffness, see Fig. 2d, and it is loaded by an imposed end-displacement u(L) = t L. Due to the non-homogeneous stiffness, we expect the damage criterion to be satisfied at first at the middle of the bar. We solve the evolution problem obtained when increasing from 0 the end displacement by minimising at each time step the corresponding energy functional [START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF] under the irreversibility constraint on the damage. The minimisation is performed using the alternate minimisation algorithm [START_REF] Bourdin | Numerical implementation of the variational formulation for quasi-static brittle fracture[END_REF] and using a finite element discretization. For the energy functional to be defined, we use a linear interpolation for the displacement and piecewise constant approximation for the damage field. This choice for the discretisation appears as the natural one because the damage field can be discontinuous in absence of any control on its derivative in the energy.

Figure 2a shows that the damage field localizes in a band of one or two elements width. Figure 2a illustrates the evolutions of the elastic and dissipated energy with the loading. Figure 2c reports the dissipated energy after failure as a function of the mesh size for both (NS) and (LS) models. The results clearly prove the pathological, and well-known, mesh-dependence: in local damage model with softening the dissipated energy depends on the mesh size and vanishes when decreasing the mesh size. This issue persists even when trying different discretisation strategies.

Damage-Gradient (DG) regularization

We briefly review the regularization of the local damage model with softening [START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF] based on the penalization of the gradient of the damage field [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF]. This kind of models is extensively used as a regularized phase-field approach to brittle fracture [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. As before, we present this model in a simplified one-dimensional setting considering the same traction test as in Section 2.

Formulation

Rate-independent damage-gradient models can be formulated in a variational framework by introducing the total energy functional 2a shows the damage and displacement fields after failure, emphasising the discretisation (markers correspond to nodes in the finite element mesh). Figure 2b reports the evolution of the energy with the loading t. Before failure (t 1.0), the elastic energy is quadratic in t (elastic response) and the dissipated energy is zero. After failure, the elastic energy vanishes and the dissipated energy is strongly dependent on the mesh size. The damage localizes in a zone of a one or two elements and the dissipated energy scales linearly with the mesh size, as shown in Figure 2c for both the (LS) and (NS) models in [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. Here the dashed lines are the energy values for damage localized in one (w 1 h) or two elements (2w 1 h). Figure 2d illustrates the initial stiffness used in the simulations:

E(u, α) = L/2 -L/2 W 0 (u (x), α(x)) + 1 2 w 1 2 α (x) 2 dx, ( 30
(α) = (1 -α) 2 E 0 and w(α) = w 1 α. Figure
E 0 (x) = {1 for |x| ≥ d/2; 1 -0.05 (2 x/d) 2 for |x| < d/2}.
where W 0 is the energy density associated with the local mode. The second term in the integral is the regularizing gradient damage term, which forbids any discontinuity of the damage field, penalises sharp damage variations, and introduces an internal length , see e.g. [START_REF] Pham | The variational approach to damage: I. The foundations[END_REF][START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF].

At each time step t i , the displacement and damage fields (u i , α i ) are defined as local minimizers of the energy [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] in the sense [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], where now the admissible spaces C i and D i are

   C i = v ∈ H 1 (-L/2, L/2) : v(-L/2) = -t i L/2, v(L/2) = t i L/2 , D i = α ∈ H 1 (-L/2, L/2) : α i-1 (x) ≤ α(x) < 1, ∀x ∈ [-L/2, L/2] .
(31) and • denotes the natural H 1 norm of C i ×D i . Note that D i contains the unilateral constraint requiring that the damage field can only increase with time, which constitutes the condition of damage irreversibility.

In a series of recent works [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Pham | From the onset of damage to rupture: Construction of responses with damage localization for a general class of gradient damage models[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] the key properties of quasi-static evolutions verifying [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], are studied performing bifurcations and stability analysis for the 1D traction problem. We briefly recall the main results in the next subsection.

Remark 1.

As far as w is monotonically increasing, by a suitable change of variables, one can always set w(α) = β to rewrite [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] in the following equivalent form, as done in [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]:

Ẽ(u, β) = L/2 -L/2 E 0 Ẽ(β) 2 u 2 + w 1 β + w 1 2 F(β) 2 β 2 dx, ( 32
)
where the new constitutive functions Ẽ and F can be univocally computed from the definition of E and w in (30).

Main properties

Throughout this subsection, we consider a given time step i ≥ 1 and we drop the index i in (almost) all the current quantities at that time. If (u, α) is a local minimum in the sense of ( 27) such that α(x) < 1 everywhere in the bar (which means that there exists no fully damaged point), then (u, α) must verify the first order optimality conditions [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF].

After integration by part and use of the classical argument of the calculus of variations one can show that (28) 1 implies the mechanical equilibrium and that the stress field is constant throughout the bar,

σ = E(α(x)) u (x), (33)
whilst ( 28) 2 gives the following Kuhn-Tucker complementary conditions at every point x:

     f (α(x), u (x), α (x)) ≥ 0, α(x) -α i-1 (x) ≥ 0, f (α(x), u (x), α (x))(α(x) -α i-1 (x)) = 0 (34) with f (α, u , α ) = 1 2 E (α) u 2 + w (α) -w 1 2 α . ( 35
)
Equation ( 34) represents the damage evolution criterion. Writing the deformation u in terms of the stress σ, (34) 1 reads in the usual form in terms of the compliance:

1 2 S (α)σ 2 ≤ w (α) -w 1 2 α . ( 36
)
The analysis of the first order optimality conditions for the uniaxial traction problem gives the fundamental results for understanding the properties of the model [24, 27, 28]:

1. Elastic limit: The bar can support loading without damage (α = 0) for stress values up to the elastic limit σ C given by (19). 2. Fully localised solutions (cracks): In the case of materials with stress softening, see [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF], there exist solutions with vanishing stress and with damage localised in a band of width ∆ 0 . Within this band α varies from 0 (undamaged material) to 1 (broken material) with an energy dissipation G c . Both the width of the localization bands and the energy dissipation are explicitly given in term of the model constitutive parameters as:

∆ 0 = 1 0 1 2w(β) dβ, G c = 2 w 1 1 0 2w(β) dβ. ( 37
)
One can assimilate these solutions with cracks, energetically equivalent to those encountered within the Griffith fracture model, where G c stands for the fracture toughness.

The results above are derived in a one-dimensional setting. The variational approach to fracture [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Francfort | The variational approach to fracture[END_REF][START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] introduces the functional of the class (30) as a regularized approximation of the brittle fracture in arbitrary dimensions. In particular, exploiting the methods developed in the variational theory of free-discontinuity problems [START_REF] Braides | Approximation of free-discontinuity problems[END_REF][START_REF] Ambrosio | On the approximation of Free Discontinuity Problems[END_REF], it can be shown that, in the limit → 0, the global minimisers of (30) converges towards the global minimisers of the energy of a Griffith fracture model. We remark that the energy [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] together with the assumptions (24) constitute a modification of the classical Ambrosio-Tortorelli assumption, [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], used in most of the mathematical literature on phase-field models of fracture:

(AT) : E(α) = (1 -α) 2 , w(α) = α 2 . ( 38
)
The main motivation to use (LS)/(NS) instead of (AT) is to have a non vanishing elastic limit, as can be immediately deduced by replacing (38) in [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF].

To illustrate the behaviour of the DG model, we report in Figure 3 the result of the same traction test described in Section 2. These are obtained minimizing [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] after a finite element discretization with piecewise linear basis functions both for displacement and damage fields, which is the simplest one for respecting the H 1 regularity of the displacement and the damage fields that guarantees a finite energy value. The key conclusions are the following:

• After failure the damage localises in a band whose width is independent of the mesh size, for sufficiently small mesh sizes.

• After failure, the dissipated energy (t 1) is independent of the loading t, and the elastic energy is vanishing.

• The dissipated energy converges to a value independent of the mesh size for the mesh size going to zero. This value corresponds to the one calculated analytically in Eq. ( 37) and defines a well-definite fracture energy for creating a crack, which is the fracture toughness G c equal to π √ 2 w 1 for the (NS) model and 4

√ 2 3
w 1 for the (LS) model.

• Numerically we observe good and smooth convergence of the alternate minimisation algorithm.

Strain-Gradient (SG) regularization

We discuss in the present section and the next one a regularized variational damage model where damage gradients are replaced by strain gradients (SG), as in [START_REF] Placidi | A variational approach for a nonlinear 1-dimensional second gradient continuum damage model[END_REF]. To this aim we study the evolution problem for a one-dimensional traction of a bar. We remark that our analysis does not include the damage models proposed in [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Peerlings | A thermodynamically motivated implicit gradient damage framework and its application to brick masonry cracking[END_REF], where the strain-gradient regularization is introduced through an additional non-local equivalent strain.

Variational formulation

In the one-dimensional traction problem, the linearized strain is ε(x) = u (x) and its gradient gives ε (x) = u (x). Replacing in (32) the damage gradient regularization term by a strain gradient term, one obtains an energy functional in the following form:

E(u, α) = L/2 -L/2 E 0 2 E(α) u 2 dx Elastic energy E1 + L/2 -L/2 E 0 2 η 2 G(α) u 2 dx Elastic energy E2 + L/2 -L/2 w 1 α dx Dissipated energy G , ( 39
)
where η is a characteristic length and G(α) a constitutive function. The energy above allows for an attractive physical interpretation of the regularization term as a nonlocal stiffness. Indeed, (39) can be regarded as the energy of second gradient continuum with a local damage variable affecting the local and nonlocal stiffness through the functions E and G, respectively. As for DG models, we require E to be a monotonically decreasing smooth function of α, with E(0) = 1 and E(1) = 0, and to satisfy the stress-softening condition and the strain-hardening condition E (α) > 0 for each α ∈ [0, 1), see [START_REF] Pham | The variational approach to damage: I. The foundations[END_REF]. Setting a damage dissipation function linear in α is not a restrictive hypothesis (Remark 1 explains how this form can be always obtained through a suitable redefinition of the damage variable).

The variational formulation of the corresponding rate-independent time-discrete evolution problem is still ruled by a local minimality condition for the energy functional. The results are for = 0.1 L and a small imperfection on the Young modulus in x = 0.6L. The width of the localisation zone and the dissipated energies converge to well-defined finite values, see (37), when decreasing the mesh size. The dissipated energies after failure do not depend on the loading. Hence it is possible to associate a well-definite constant dissipated energy to a crack and the model is energetic equivalent with the Griffith fracture model.

The two fundamental families of SG models

For α → G(α), we can distinguish two families of functions that give rise to fundamentally different behaviors:

1. The first family (SG1) consists in functions G that are positive and bounded from below by a strictly positive constant for any value of α ∈ [0, 1]. For such functions the straingradient stiffness never disappears. To simplify the presentation we will only consider, in the present section, the particular case where G(α) is a constant. 2. The second family (SG0) contains the functions G that vanish at α = 1 and hence that lead to a loss of any stiffness at α = 1, the material being unable to sustain any stress. This family will be studied in the next section.

The first family of SG models forces an additional regularity of the strain field but a priori allows for less regular damage fields. Indeed, having a finite value of the functional (43) requires the strain to have square integrable first derivatives. On the other hand, the condition of finite energy allows the damage field to be discontinuous, in contrast with the case of DG models. Hence, we consider the following sets of admissible displacement and damage fields at time i:

   C i = v ∈ H 2 (-L/2, L/2) : v(-L/2) = -t i L/2, v(L/2) = t i L/2 , D i = β ∈ L ∞ (-L/2, L/2) : α i-1 ≤ β ≤ 1 in (-L/2, L/2) , ( 40
)
which differ from their homologue (31) of the DG models by the new regularity assumptions dictated by the finite energy condition. The linear space associated with the affine space C i is

C 0 = {v ∈ H 2 (-L/2, L/2) : v(±L/2) = 0}. ( 41
)
The local minimality condition for the displacement and damage fields (u i , α i ) can still read as [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], but the natural norm entering in the statement is now the H 2 norm for the displacements fields and the L ∞ norm for the damage fields.

For the second family of models, where G(1) = 0, the discussion of the regularity is subtler. In all the points where α = 1 the regularizing effect of the strain-gradient term vanishes and displacement jumps become possible.

The following sections 5 and 6 are devoted to the analysis of the two families of models.

SG models with constant strain-gradient stiffness

We consider here the case in which the nonlocal stiffness G(α) is a constant, independent of α. Without loss of generality, we take

G(α) = 1, (42)
and let η determine the magnitude of the nonlocal stiffness. Let us note that η has the dimension of a length and thus can be seen as the characteristic length of the SG model. Accordingly, throughout all the present subsection, the energy functional reads as

E(u, α) = L/2 -L/2 E 0 2 E(α) u 2 + E 0 2 η 2 u 2 + w 1 α dx. ( 43
)
To illustrate the behaviour of this class of models, we extend the analysis of the previous sections by deducing from the local minimal principle the first order optimality conditions that any solution must satisfy and determining the analytical solutions at fixed loading that are either homogeneous or localised in space.

First order optimality conditions. 5.1.1. Case of solutions with no fully damaged points.

Let us first consider a local minimum in the sense of [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], denoted (u i , α i ), such that α i < 1 everywhere in the bar. Then (u i , α i ) must still verify the first order optimality conditions [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], but because of the change of the energy functional, these conditions lead now to

L/2 -L/2 E 0 E(α i ) u i v + η 2 u i v dx = 0, ∀v ∈ C 0 (44a)
and to

L/2 -L/2 E 0 2 E (α i ) u 2 i + w 1 (β -α i ) dx ≥ 0, ∀β ∈ D i . ( 44b
)
For smooth solutions, applying classical arguments of the calculus of variations, condition (44a) implies the equilibrium equation and the natural boundary conditions

σ i = 0, in (-L/2, L/2), u i = 0 at ± L/2, ( 45
)
where the stress σ i is now given by

σ i = E 0 E(α i ) u i -η 2 u i , ( 46
)
which includes the term -E 0 η 2 u i coming from the by parts integration of the strain-gradient term in (44a).

In the case of non smooth solutions, (45) and ( 46) remain valid provided that the derivatives are interpreted in a weak sense. Accordingly, recalling that one assumes that α i < 1 everywhere, at a generic non-regular point in (-L/2, +L/2) where one investigates possible jump discontinuities, one deduces that the following continuity conditions must hold:

u i = 0 u i = 0, u i = 0, σ i = 0. (47)
In (47) • denotes the jump of the function • at the considered point. Therefore the stress is necessarily uniform all along the bar:

σ i (x) = σ i , ∀x ∈ (-L/2, L/2).
Let us note that the continuity of the strain field and the uniformity of the stress field will induce (in general) the continuity of the damage field. In other words, the optimality conditions force the regularity of the real damage field.

Still classical arguments of the calculus of variations show that the first order optimality condition (44b) is equivalent to the following set of local Kuhn-Tucker conditions at each point x where α i (x) < 1 (hence, at that point α i-1 (x) < 1):

       0 ≤ E 0 E (α i (x)) u i (x) 2 + 2w 1 0 ≤ α i (x) -α i-1 (x) 0 = E 0 E (α i (x)) u i (x) 2 + 2w 1 α i (x) -α i-1 (x) (48)
which are nothing but the damage evolution law at a not fully damaged point.

Case of solutions with fully damaged points.

Even if the solution contains fully damaged points, the Kuhn-Tucker conditions (48) remain valid at each point which is not fully damaged, i.e. at each x such that α i (x) < 1.

Let us now consider a point x 0 which is already fully damaged at step i -1, i.e. such that α i-1 (x 0 ) = 1. The irreversibility condition imposes that α i (x 0 ) = 1 and then, since every test field β ∈ D i must be such β(x 0 ) = 1, the optimality condition (44b) is automatically satisfied.

Finally, let x 1 be a point fully damaged at step i but which was not fully damaged at step i -1:

α i-1 (x 1 ) < α i (x 1 ) = 1. (49)
Then, by considering admissible test fields β such that β(x 1 ) < 1, the optimality condition (44b) gives

u i (x 1 ) 2 ≥ 2w 1 E 0 |E (1)| (50)
which can be seen as the necessary condition for a point to become fully damaged.

Homogenous solutions and their bifurcation analysis

The homogenous response presented above satisfies the first order minimality conditions (44) for the energy functional. During the elastic phase, the damage criterion is not attained and the condition (44b) is a strict inequality for any non vanishing variation β of damage. This implies that for ε < ε C the homogenous solution is a local minimum of the energy. For ε ≥ ε C , the damage criterion is satisfied as an equality everywhere in the domain. In this case, assessing the stability of the homogenous response requires the study of the sign of the second derivative of the energy functional near u(x) = ε x and α(x) = α H (ε). After relative cumbersome, but standard, calculations which are presented in AppendixC, one can show that the second derivative of the energy is positive definite for any admissible variation of displacement and damage if and only if the bar is sufficient short with respect to the internal length of the material, namely

L η < π E (α H (ε)) 2E (α H (ε)) 2 -E(α H (ε))E (α H (ε)) =    (LS) model : π √ k -1, (NS) model : π 4 √ 3 ε ε C (51)
The condition above on the sign of the second derivative of the energy gives a criterion of non-bifurcation (or local uniqueness) for the homogeneous solution as a function of the ratio L/η and the strain level ε. It also provides a sufficient condition for stability; to establish a necessary condition for stability would require the same analysis to be performed under the constraint of positive variations of the damage variable, β ≥ 0 in AppendixC. We refer the interested reader to [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] for further details on such an issue.

For the constitutive choices of E(α) given in [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF], Figure 4 reports the regions in the plane (L/η, ε) where the homogeneous response is (locally) unique. For the (LS) model, the homogenous response is always unique for short bars (L/η < π √ k -1) but it could bifurcate for ε > ε C for long bars (L/η > π √ k -1). Instead for the (NS) model the homogenous solution turns out to be (locally) unique at large strains, even for long bars. Even if not explicitly reported here, the stability analysis would give qualitatively similar results.

Non-homogenous solutions at fixed loading. 5.3.1. Assumptions

For sufficiently long bars, the homogenous solutions can become unstable and solutions with damage localization could emerge. We study here the properties of non-homogenous solutions of the traction problem in which the damage is non-null in a single non-vanishing interval. Specifically we make the following set of assumptions, see also Figure 5:

(1) To simplify the presentation we will only consider infinitely long bars, L = +∞. Accordingly, the boundary conditions are replaced by the conditions that the strain field ε i must tend at infinity to a given constant εi , lim

x→±∞ ε i (x) = εi . ( 52
)
(2) The construction of the localized solution starts at time t 1 such that the bar is entirely undamaged but the strain is just at the elastic limit:

ε 1 (x) = ε C , α 1 (x) = 0.
(3) At time step i > 1, by a suitable choice of the coordinates, damage is localized in an interval I i = (-D i , D i ), centered at 0 and the size of which, 2D i , has to be determined.

(4) In the (possibly empty) part J i ⊂ I i the bar is fully damaged : α i = 1.

(5) In the remaining part of I i , the bar is damaging but remains partially damaged: 0 ≤ α i-1 < α i < 1.

Figure 5: Sketch of the damage field in the considered non-homogenous solution.

Optimality conditions

Equipped with these assumptions and all the first order optimality conditions obtained in Subsection 5.1, let us formulate the conditions that the non-homogeneous solution must satisfy in each part of the bar at the time step i > 1.

(i) All along the bar, by the equilibrium condition, the stress is constant: σ i (x) = σ i . (ii) In the undamaged part of the R I i , the stress-strain relation, the damage criterion and the conditions at infinity read as

     E 0 (ε i (x) -η 2 ε i (x)) = σ i , |ε i (x)| ≤ ε C , lim x→±∞ ε i (x) = εi .
Solving the first differential equation, we deduce that ε i (x) decades exponentially as x → ±∞ and that

σ i = E 0 εi . ( 53
)
Note that the damage criterion requires εi to be less than ε C , necessary condition that we adopt henceforth:

0 < εi < ε C . ( 54
)
(iii) In the fully damaged part J i , when it exists, E(1) = 0 and the stress-strain relation and the damage criterion read as

     -η 2 ε i (x) = εi , |ε i (x)| ≥ ε L , α i (x) = 1. (55)
(iv) In the partially damaged part I i \ J i , the stress-strain relation and the damage criterion read as

     E(α i (x))ε i (x) -η 2 ε i (x) = εi , E (α i (x))ε i (x) 2 = E (0)ε 2 C , 0 ≤ α i-1 (x) < α i (x) < 1. ( 56
)
Let us first deduce some continuity properties for the localized solution. Since ε i is continuous everywhere, α i is continuous in I i \J i . Moreover, |ε i | must be equal to ε C at the boundary points x = ±D i of I i , and |ε i | must be equal to ε L at the boundary of J i . Therefore α i (±D i ) = 0, α i = 1 at the boundary of J i and the damage field is continuous everywhere.

Let us now determine the non-homogeneous solution. Let us remark that, at any point x, the damage is related to the strain by

α i (x) = α H (ε i (x)), (57)
with α H given by [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF]. Inserting that relation into the stress-strain relation leads to the following second order differential equation in terms of the strain field:

E(α H (ε i (x)))ε i (x) -η 2 ε i (x) = εi , ∀x ∈ R (58)
with the conditions at infinity ε i (±∞) = εi . Accordingly, it suffices to solve this differential equation for the strain field to obtain the damage field and hence the parts I i and J i .

The determination of the solution of ( 58) is based on the two following properties the proofs of which are reported in AppendixB. Proposition 5.1 (First integral). For a given εi ∈ (0, ε C ), the differential equation (58) with the condition at infinity (52) for the strain field ε i admits the first integral

Qualitative properties of the solutions

η 2 ε i (x) 2 = F i (ε i (x)), ∀x ∈ R, ( 59
)
where

F i (ε) = 2W H (ε) E 0 -2ε i ε + ε2 i . ( 60
)
Moreover the function F i verifies the following properties, as illustrated in Figure 6:

P1. For |ε| ≤ ε C , F i (ε) = (ε -εi ) 2 . P2. There exists a unique ε * * i ∈ (ε C , ε L ) such that F i (ε * * i ) = 0. P3. There exists a unique ε * i > ε C such that F i (ε * i ) = 0. P4. The value of ε * i increases monotonically from ε C to ∞ when εi decreases from ε C to 0.
The first integral allows us to establish the most important qualitative properties of the strain and damage fields at a given time step. Those properties are given in the following Proposition and illustrated in Figure 8.

Proposition 5.2 (Qualitative properties of the non-homogenous solutions). At time step i > 1

where εi < ε C , there exists a unique strain field solution to (59) and the associated damage field is given by ( 14) and (57). This solution is symmetric with respect to the origin and enjoys the following properties:

(i) The strain and damage fields are maximal at x = 0, their maximum value being ε * i defined in Proposition 5.1 and α H (ε * i ), respectively. Then they are monotonically decreasing with respect to |x|.

(ii) The size D i of the damaged zone I i is given by

D i = η ε * i εC 1 F i (ε) dε. ( 61
)
(iii) Outside the damaged zone, the strain field is given by

ε i (x) = εi + (ε C -εi ) exp - |x| -D i η when |x| ≥ D i . ( 62
)
(iv) Inside the damaged zone I i , the strain ε i (x) at the point x is given by the implicit relation

D i -|x| = η εi(x) εC 1 F i (ε) dε when |x| < D i . ( 63
) (v)
The existence of a fully damaged zone J i depends both on ε L and εi . Specifically, one has (a) If ε L = +∞, then there exists no fully damaged point whatever the value of εi

J i = ∅;
(b) If ε L < +∞, then there exists fully damaged points provided that εi is small enough so that ε * i > ε L . In such a case the fully damaged zone is a closed interval:

J i = [-Di , Di ] with Di = η ε * i εL 1 F i (ε) dε. ( 64
) 0 Figure 7: Plot of the localised solution curve ηε i = ± F i (ε i ) in the (ε, ηε ) plane when ε * i > εL.
The color bar at the top of the plot gives the value of the damage field α i calculated as a function of ε i using [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF]. The qualitative properties of this curve are given in Proposition 5.1.

Example. We consider the (LS) model under its form [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. Using the definitions introduced previously, we get

ε C = 2w 1 kE 0 , ε L = k ε C . ( 65
)
Since ε L is finite, there exists a critical value ε M of the applied strain at infinity under which a fully damaged zone exists. That value is given by

ε M = k -k (k -1) ε C . ( 66
)
For a given applied strain at infinity εi ∈ (0, ε C ), using Proposition 5.2 leads after some calculations to the following expressions for the maximal strain and the size of the damaged zones:

• When ε M ≤ εi < ε C ,          ε * i = ε L - ε L ε M -1 (ε i -ε M ), D i = D c := 2η √ k -1 arctan √ k + √ k -1 , Di = 0, ( 67
)
where we can note that the size of the damaged zone is fixed, independent of εi .

• When 0 < εi < ε M ,                    ε * i = ε L + 1 2 ε L ε C ε M εi -1 (ε M -εi ) D i = D c -2η √ k -1 arcsin 1 -εi /ε M 2(1 -εi /ε C ) + Di , Di = η ε L ε C ε M εi -1 ε M εi -1 . ( 68
)
In both cases, the strain field in the undamaged zone is exponentially dependent on |x| and given by (62). In the partially damaged zone, since σ H (ε) depends linearly on ε,

σ H (ε) = σ c k -1 k - ε ε c ,
the differential equation (58) for ε i is linear and therefore both the strain field and the damage fields are arcs of sinusoid as one can see on Figure 8. In the fully damaged zone, when it exists, the strain field is an arc of parabola as one can check directly from (55).

Asymptotic behavior when the stress goes to 0.

For the applications it is relevant to study the behaviour of the solution when εi → 0 (or equivalently when σ i → 0). Indeed, the solution with vanishing stress is the one that can be associated to a crack in the Griffith fracture model.

Proposition 5.3 (Asymptotic behaviour for vanishing stress). εi → 0 the non-homogenous solutions show the following asymptotic behaviour:

• The maximum strain value grows to infinity:

lim εi→0 ε * i = +∞. ( 69
)
• The size of the damaged interval grows to infinity:

lim εi→0 D i = +∞. ( 70
)
• The dissipated energy also tends to infinity:

lim εi→0 G i = ∞, G i := Di -Di w 1 α i (x) dx. ( 71
)
Proof. In the case of the (LS) model, these properties are straightforward consequences of (68). They remain valid for the general family of SG models studied in the present section and the proof is essentially based on Propositions 5.1 and 5.2.

• The property (69) is nothing but P4 of Proposition 5.1.

• To prove (70), let us start from the expression (61) of

D i . When ε C ≤ ε ≤ ε * i , since 0 ≤ F i (ε) ≤ F i (ε * * i ) ≤ w 1 , one has 1/ F i (ε) ≥ 1/
√ w 1 and, since the domain of integration goes to infinity, one gets (70).

• When ε L < ∞, by the same argument, (64) gives lim εi→0 Di = +∞ and (71) follows.

The case ε L = ∞ is more subtle since we need to bound from below α i (x). From ( 61) and (63), one gets

|x| η = ε * i εi(x) dε F i (ε) ≥ ε * i -ε i (x) √ w 1 when |x| > D i .
Therefore, for a given x ∈ R, since lim εi→0 ε * i = +∞ and lim εi→0 D i = +∞, one gets lim εi→0 ε i (x) = +∞. But, since α i (x) = α H (ε i (x)), one gets also lim εi→0 α i (x) = 1 and hence (71).

Let us note that the last estimate on the damage field says that every point of the bar tends to a fully damaged state when the stress goes to 0.

SG models with generic strain-gradient stiffness

We extend here the results of the previous Section by considering the generic case where the strain-gradient stiffness can vary with the damage field, namely G = G(α). We focus on the case in which G(α) is monotonically decreasing5 with α and distinguish two fundamentally different behaviour whether G(1) > 0 or G(1) = 0. In particular, we study the global minimum of the energy and we show that it is either infinite or zero. The consequences of this fact are illustrated in the numerical solution of the traction test on a bar with finite length.

SG1: the generic case with

G(1) > 0.
case can be easily reduced to the case studied in Section 5 with a constant straingradient stiffness. Indeed, denoting by Ḡ := min{G(α), α ∈ [0, 1]}, we have

L/2 -L/2 G(α) (u ) 2 dx ≥ L/2 -L/2 (u ) 2 Ḡ dx, ∀u, ∀α. (72)
Hence, the total energy (39) of the model with constant strain-gradient stiffness G(α) = Ḡ can be used as a lower bound for the generic case with a monotone decreasing G(α). Using the results of Proposition 5.3, we conclude that, as for G(α) = Ḡ, creating a crack or obtaining a solution with vanishing stresses in a bar with infinite length requires an infinite amount of energy, even if G(α) is not constant but G(1) > 0.

6.2. SG0: the generic case with G(1) = 0.

The functional (39) to minimize is not convex and can admit multiple local and global minimizers. To characterise these minima, in Sections 3-4 we solved the first order optimality conditions [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]. Here we follow a different approach based on the direct method of the calculus of variations, which is particularly effective to determine the properties of global minimizers of the energy.

We consider the minimization problem for

E(u, α) = L/2 -L/2 E 0 2 E(α) (u ) 2 + η 2 G(α) (u ) 2 + w 1 α dx, (73)
with E(1) = 0, G(1) = 0, the boundary conditions u(±L/2) = ±t L/2 and the unilateral constraint α(x) ≥ 0. For any value of the end-displacement it is possible to define the following admissible test field (see Figure 9):

u δ (x) =    t L 2 1 + 2x δ + sin(2π x/δ) π -δ/2 ≤ x ≤ +δ/2, 0, otherwise (74)
α δ (x) = 1, -δ/2 ≤ x ≤ +δ/2, 0, otherwise (75)
for which

ε δ (x) = u δ (x) =    2 t L δ (cos(πx/δ)) 2 , -δ/2 ≤ x ≤ δ/2, 0, otherwise (76)
It is immediate to verify that u δ , ε δ , and ε δ are continuous functions of x, while α jumps in x = ±δ/2. Being ε δ (±δ/2) = 0, all the jumps conditions (47) are verified in x = ±δ/2. Hence (u δ , α δ ) is an admissible test field for (73), for which the energy ( 73) is evaluated to be

E(u δ , α δ ) = w 1 δ, lim δ→0 E(u δ , α δ ) = 0. ( 77
)
Since E(u, α) ≥ 0 for any admissible state, (u δ , α δ ) is a minimising sequence attaining the global minimum of the energy for δ → 0, for any loading t. We conclude that for E(1) = G(1) = 0, the global minimum of the energy is 0, and this global minimum can be obtained through solutions localizing the damage field in a zone of vanishing width, where all the deformations are concentrated, without any energy dissipation. Remark 2. If G(1) > 0, the energy of the sequence (74) is given by

E(u δ , α δ ) = w 1 δ + E 0 η 2 G(1) π 2 t 2 L 2 δ 3 . ( 78
)
This expression emphasizes the crucial role played by the nonlocal stiffness at full damage and shows that the reasoning above cannot be applied for the models where G(1) > 0.

The fact that the global minimum of the energy is 0 for any loading lead us conclude that the strain gradient regularization is not effective when G(1) = 0, even if G(α) > 0 for α < 1. Indeed, even if there could be still meta-stable states that are local minima of the energy (an example of these is the elastic solution for t < ε C ), there are not stable states (in the sense of global minimization) that can sustain a finite stress.

Numerical illustrations

To confirm the theoretical results we consider two specific classes of models, giving the cases of gradient stiffness decreasing to zero with the damage and of constant gradient stiffness, respectively:

(SG0)-models : G(α) = E(α) ⇒ G(1) = E(1) = 0, (79)
(SG1)-models :

G(α) = 1. (80)
Hence we will denote by NS-G0 the model obtained when choosing E(α) = G(α) = (1 -α) 2 and by NS-G1 when E(α) = (1 -α) 2 but G(α) = 1. Similar notations are used for the LS models as in [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. The numerical solution is obtained by using an augmented lagrangian formulation keeping as primal variables the displacement u, the strain ε, and the damage field α. The compatibility between the displacement and the strain is enforced by using a scalar Lagrange multiplier field τ . Hence, we use linear Lagrange finite elements for ε and α, quadratic Lagrange elements for the displacement u, and elements with piecewise constant interpolation functions for τ . The energy functional is minimised by finding at each time step a fixed point of an alternate minimization algorithm where we solve iteratively the linear system for u, α, τ at fixed α, and the complementary system (a variational inequality) for α at fixed u, α, τ . The code is implemented using the FEniCS library. We refer the reader to the sources 6 for further details about the implementation.

We consider a traction test on a bar which is completely clamped at the two end, by imposing to 0 the variations of the displacement and the strain fields. The loading is given by an imposed end-displacement u(L) = t. Figure 10 reports the numerical results obtained using the NS-G1 model. In particular, in Figure 10a the evolution of the damage and strain fields are respectively reported by solid and dashed lines for different values of the load parameter. Figure 10b plots the relevant energetic contributions in (39): the elastic energy (E 1 , dotted), the nonlocal elastic energy (E 2 , dashed) and the dissipated energy (G, black solid) are shown. The analogous results for the LS-G1 model are given in Figure 11.

These results are coherent with their analytical counterparts shown in Figure 8 and with Proposition 5.3: as the load increases the width of the damaged region progressively increases to cover the whole bar and the dissipated energy grows towards its maximal value (w 1 L ≡ 1 in the numerical simulations on a bar of finite length). Figure 12 reports the dissipated energy at the end the load process for models from both the families (SG0 and SG1) of Strain-Gradient models. The NS-G1 and LS-G1 models with constant strain-gradient stiffness show a limit value of the dissipated energy which corresponds to the full damage of the whole bar; this value is independent of the mesh size. Instead, the dissipation in the NS-G0 and LS-G0 models are strongly mesh-dependent. Being G(1) = 0 the regularization is not effective. As in the local models, see Figure 2, the damage localizes in few elements and the dissipated energy is actually determined by the mesh-size h and vanishes when h → 0. In general the numerical results of the NS-G0 and LS-G0 models are not reliable, being extremely sensitive to the mesh, imperfections, and the parameters of the numerical solvers. This point is further illustrated by the snapshots of the damage and strain fields reported in Figure 13.

Conclusions

In brittle fracture models the energy required to create a crack of unit length is a material property, the fracture toughness. It is finite, non-null, and independent on the structural size (at least for large enough structures). In damage models a crack can be defined as a solution with vanishing elastic energy, obtained after a softening phase. Depending on the adopted damage model, the internal energy dissipation required to create a crack may be null or meshdependent, constant, or infinite. However, establishing an energetic equivalence of a damage model with brittle fracture models requires to fulfil the following requirements: 23 P1 The dissipated energy is mesh-independent, for sufficiently refined mesh. P2 The dissipated energy is non-vanishing. P3 The dissipated energy is finite and independent of the structural size, for sufficiently large structures.

Table 1 resumes the fulfilling of these properties for the damage models considered.

LOC DG SG1 SG0 P1 P2 P3 In local damage model (LOC), the dissipated energy is proportional to the mesh size h, see Figure 2c, and therefore vanishing for h → 0. Regularization techniques are aimed to overcome this drawback. We have shown that the damage-gradient DG models can be effective to this aim and verify all the desidered properties P1-P3, see Figure 3c. In DG models the toughness has a well defined nonvanishing value independent of the mesh-size and the bar length, see (71), that can be used to the identify the fracture toughness with the material parameters of the damage model.

We have shown that a similar identification is not possible for strain-gradient (SG) models in the form (39).

If the nonlocal stiffness function G(α) is bounded from below by a strictly positive constant, we have found that the stress vanishes only when the whole bar is fully damaged, see NS-G1 and LS-G1 curves in Figure 12. As the toughness is evaluated to be w 1 L, SG1 models fail to satisfy condition P3. For a bar of infinite length, we have proved that an infinite amount of energy is needed to create a crack, see Proposition 5.3.

If the nonlocal stiffness vanishes when the damage reaches its maximum value, G(α = 1) = 0, the dissipated energy can be mesh-dependent as in local models, see Section 6.2 and NS-G0 curve in Figure 12 (SG0 models).

These results are valid under the choice of a quadratic cost of the strain-gradient fields. This choice was regarded as the simplest since leading, at given damage, to a linear elastic constitutive relation between the effective stress and the gradient of strain. If other subtler choices to penalize the strain-gradients, see for instance [START_REF] Forest | Nonlinear regularization operators as derived from the micromorphic approach to gradient elasticity, viscoplasticity and damage[END_REF][START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF], are possible and lead to satisfactory properties of the resulting regularized model is not known to the authors.

AppendixA. Ill-posedness of local damage model

Proof. By hypothesis there is a set S ⊂ [-L/2, L/2] =: Ω of non-vanishing measure where the damage criterion is satisfied as an equality:

E (u i , α i )(v, ϕ) = 0, ∀v, ∀ϕ ≥ 0 in S, ϕ = 0 in Ω/S. (A.1)
The perturbed energy E(u i + hv, α i + hϕ) can be expanded with respect to h as:

E(u i + hv, α i + hϕ) = E(u i , α i ) + h 2 2 E (u i , α i )(v, ϕ) + o(h 2 ).
The stability depends on the sign of the second derivative of the energy, namely

E (u i , α i )(v, ϕ) E 0 = Ω E i v 2 dx + S 2E i u i v ϕ + 1 2 E i (u i ) 2 ϕ 2 + ϕ 2 w i dx = Ω E i v 2 dx + S 2E i u i v ϕ + (w i -w i E i /E i ) ϕ 2 dx (A.2) with E i = E(α H (t i )), E i = E (α H (t i )), E i = E (α H (t i )), w i = w (α H (t i )) and w i = w (α H (t i ))
. The second derivative is minimized at given ϕ with respect to v; using the boundary conditions v(±L/2) = 0 we obtain

v (x) = C E i (x) - E i (x)u i (x) E i (x) ϕ(x), C := S ϕ(ξ)E i (ξ)u i (ξ) E i (ξ) dξ Ω 1 E i (ξ) dξ . (A.3)
Substituting in (A.2), the energy second derivative reads

E (u i , α i )(v, ϕ) E 0 = Ω C 2 E i (ξ) dξ + S [w i -w i (E i /E i -2E i /E i )] ϕ 2 (ξ) dξ. (A.4)
By virtue of the stress softening condition [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF], the following inequality holds

w i -w i (E i /E i -2E i /E i ) < 0.
The second integral in (A.4) is negative for every perturbation ϕ and is linearly proportional to the measure, say θ, of the support set of ϕ. The first integral in (A.4) is instead positive but it is a quadratic function of the measure θ. Hence, choosing a ϕ with a sufficiently small measure of its support will produce a negative value of the energy second derivative. Therefore the unstable character of every solution (u i , α i ) verifying (A.1), and the need for regularizing terms.

AppendixB. Proof of Propositions 5.1 and 5.2

Sketch of the proof of Proposition 5.1. Throughout the proof, ε denotes the strain and accordingly ε(x) = u (x). Multiplying (58) by ε (x), integrating with to respect to x, and using equations ( 16)-( 18) one obtains:

η 2 ε (x) 2 = 2W H (ε(x)) E 0 -2ε i ε(x) + c, (B.1)
where c is an integration constant. Evaluating (B.1) at infinity with the boundary condition (52) and using [START_REF] Marigo | Constitutive relations in plasticity, damage and fracture mechanics based on a work property[END_REF] give the constant c:

c = 2ε 2 i -2W H (ε i )/E 0 = ε2 i .
Inserting that expression of c into (B.1) leads to (59)-(60). The function ε → F i (ε) is continuously differentiable everywhere and is twice differentiable everywhere except at ε ∈ {±ε C , ±ε L }. Denoting the derivative with respect to ε by a prime and using ( 16)-( 18) lead to

F i (ε) = 2E(α H (ε))ε -2ε i . (B.2) (i) If |ε| ≤ ε C , then α H (ε) = 0, 2W H (ε) = E 0 ε 2 and P1 follows. (ii) If ε C < ε < ε L , then σ H (ε)
is a decreasing function of ε by virtue of the softening condition [START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF] and hence F i is also decreasing. Since

F i (ε C ) = 2(ε C -εi ) > 0 and F i (ε L ) = -ε i < 0, then there exists a unique ε * * i ∈ (ε C , ε L ) such that F i (ε * * i ) = 0. (iii) By virtue of the previous property, F i is increasing in the interval (ε C , ε * * i ). Since F i (ε C ) = (ε C -εi ) 2 > 0, one gets F i (ε * * i ) > 0. We have also shown that F i is decreasing in the interval (ε * * i , ε L ). Moreover, if ε L < +∞ and ε ≥ ε L , then α H (ε) = 1 and F i (ε) = -2ε i < 0. Therefore, F i is decreasing in the interval (ε * * i , +∞).
Under the assumption (21), one gets lim ε→∞ W H (ε) = w 1 < ∞. Therefore lim ε→∞ F i (ε) = -∞ and there exists a unique ε * i > ε C such that F i (ε * i ) = 0.

25

(iv) Let us consider F i (ε) as a function of (ε i , ε), say

F (ε i , ε) := 2W H (ε) E 0 -2ε i ε + ε2 i .
By construction, ε * i satisfies ε * i > ε C > εi and F (ε i , ε * i ) = 0. By the implicit function theorem, we get

dε * i dε i = - ∂F ∂ εi (ε i , ε * i ) ∂F ∂ε (ε i , ε * i )
.

Therefore, since ∂F ∂ εi (ε i , ε * i ) = 2(ε i -ε * i ) > 0 and since, as we have shown in the previous property, ∂F ∂ε (ε i , ε * i ) < 0, ε * i is increasing when εi is decreasing. We immediately deduce from property P1 that F (ε C , ε C ) = 0 and hence ε * i = ε C when εi = ε C . The proof of the asymptotic behavior of ε * i when εi → 0 is given in Proposition 5.3. All these properties can be easily retrieved by geometrically intepreting F i (ε) as the algebric sum of the shaded areas reported in Figure B.14. For instance, chosen a level of stress σ < σ C , the critical strain ε = ε * i is the one leading to the equality between the areas marked as ⊕ and .

Sketch of the proof of Proposition 5.2. The localized solution must obey

ηε i (x) = ± F i (ε i (x)), (B.3)
which, by virtue of the properties of the function F i for εi < ε C , corresponds in the phase plane (ε, ε ) to the unique curve shown Figure 7. Since the damage is maximal where the strain is maximal and since we assume that the damage is maximal at x = 0, the strain is maximal at x = 0 where it takes the value ε * i . Accordingly, ones has ηε (x) = -sign(x) F i (ε i (x))

from which we deduce that

|x| = η ε * i εi(x)
dε F i (ε) .

The damage criterion dictates that α = 0 for ε ≤ ε C and α > 0 for ε > ε C . Therefore the length of the damage zone is given by (61).

By virtue of the definition [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF] of α H (ε), the fully damaged zone, when it exists, corresponds to the points x where ε i (x) ≥ ε L . Therefore its length is given by (64). The remaining properties are easily established.

AppendixC. Bifurcation analysis of homogeneous states for SG models

This appendix is devoted to the study the bifurcation of homogeneous solutions for SG models. The definition of stability of a state that we adopt here corresponds to the condition of local minimality of the energy [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] introduced in Section 2. We merely give the main steps of the calculation and the reader interested by a more detailed presentation should refer to [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] where a comprehensive analysis of a similar case is performed.

In the case of a SG model with G(α) = 1, the energy of the state (u, α) is given by and its first and second directional derivatives respectively read as Let t i be the given time step (which corresponds to the prescribed average strain) and let us consider the associated homogenous solution u(x) = t i x, α(x) = α H (t i ) with α H (t i ) given by [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF].

E(u, α) = L/2 -L/2
E (u, α)(v, β) E 0 = L/2 -L/2 E(α) u v + η 2 u v + 1 2 E (α)u 2 + w 1 E 0 β dx, (C.
We only consider the case when ε C < t i < ε L so that 0 < α H (t i ) < 1. Let (v, β) be an admissible direction of perturbation of the homogeneous state at that time step: v must belong to C 0 = {v ∈ H 2 (-L/2, L/2) : v(±L/2) = 0} and β must belong to L ∞ (-L/2, L/2) so that (u + hv, α + hβ) belongs to C i × D i for h positive and small enough. The homogeneous state (u, α) is a local minimum of the energy only if there exists h > 0 such that the following inequality holds true:

E(u + hv, α + hβ) ≥ E(u, α), ∀h ∈ [0, h].
Accordingly, let us expand the perturbed energy E(u + hv, α + hβ) with respect to h up to the second order:

E(u + hv, α + hβ) = E(u, α) + hE (u, α)(v, β) + h 2 2 E (u, α)(v, β) + o(h 2 ).
A straightforward calculation shows that the first derivate of the energy at (u, α) vanishes in any admissible direction of perturbation, i.e. E (u, α)(v, β) = 0. Therefore the stability of (u, α) depends on the sign of the second derivative of the energy. That derivative in the direction of perturbation (v, β) can read as

E (u, α)(v, β) = E 0 L/2 -L/2 E i v 2 + η 2 v 2 + 2E i t i v β + 1 2 E i t 2 i β 2 dx, (C.4)
where E i = E(α H (t i )), E i = E (α H (t i )) and E i = E (α H (t i )). In order to find under which condition the second derivative is positive for all admissible directions of perturbations, let us first minimize it with respect to β at given v. We immediately get that the minimum is obtained when

β(x) = - 2E i t i E i v (x).
Inserting into (C.4), the second derivative becomes the following quadratic form of v

E (u, α)(v, β) = E i (v) := E 0 L/2 -L/2 η 2 v 2 - 2E i 2 E i -E i v 2 dx.
Therefore, since by virtue of the strain hardening and stress softening conditions, the following inequalities hold

E i > 0, 2E i 2 -E i E i > 0,
the quadratic form E i (v) is the difference of two positive terms. Then, using the following classical result (the verification of which is left to the reader)

min v∈C0 L/2 -L/2 v 2 dx L/2 -L/2 v 2 dx = π 2 L 2 ,
one immediately obtains that the second derivative of the energy is positive in any non trivial direction of perturbation if and only the following inequality holds

L η < π E (α H (t i )) 2E (α H (t i )) 2 -E(α H (t i ))E (α H (t i )) . (C.5)
This inequality constitutes the stability condition of the homogeneous solution at time t i .

Figure 1 :

 1 Figure 1: Homogeneous responses for the (LS), left, and (NS), right, models, for which εL = kεC and εL = ∞, respectively.

)

Figure 2 :

 2 Figure 2: Numerical results for the non-homogenous traction test with the local model characterized by E(α) = (1 -α) 2 E 0 and w(α) = w 1 α. Figure2ashows the damage and displacement fields after failure, emphasising the discretisation (markers correspond to nodes in the finite element mesh). Figure2breports the evolution of the energy with the loading t. Before failure (t 1.0), the elastic energy is quadratic in t (elastic response) and the dissipated energy is zero. After failure, the elastic energy vanishes and the dissipated energy is strongly dependent on the mesh size. The damage localizes in a zone of a one or two elements and the dissipated energy scales linearly with the mesh size, as shown in Figure2cfor both the (LS) and (NS) models in[START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. Here the dashed lines are the energy values for damage localized in one (w 1 h) or two elements (2w 1 h). Figure2dillustrates the initial stiffness used in the simulations: E 0 (x) = {1 for |x| ≥ d/2; 1 -0.05 (2 x/d)2 for |x| < d/2}.

  Fields after failure (DG-LS model) Energy vs loading (DG-LS model) (c) Dissipated energy vs mesh size

Figure 3 :

 3 Figure3: Damage-Gradient model: numerical results for the traction test for the (LS) and (NS) models in[START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. The results are for = 0.1 L and a small imperfection on the Young modulus in x = 0.6L. The width of the localisation zone and the dissipated energies converge to well-defined finite values, see (37), when decreasing the mesh size. The dissipated energies after failure do not depend on the loading. Hence it is possible to associate a well-definite constant dissipated energy to a crack and the model is energetic equivalent with the Griffith fracture model.

Figure 4 :

 4 Figure 4: Conditions for the local uniqueness of the homogenous response as a function of the loading (ε/εC) and the normalized bar length (L/η). Light gray indicates regions where the homogeneous response is locally unique, hence stable. Left: (LS) model; right: (NS) model.

Figure 6 :

 6 Figure 6: Qualitative behavior of the function F i (ε).

Figure 8 :

 8 Figure 8: Evolution of the strain and damage fields for the LS model with k = 2 (Left) and for the NS model (Right) when εi is decreasing from εC to 0.

Figure 9 :

 9 Figure 9: Minimising sequence (74) for the SG0 model with E(1) = G(1) = 0. Black solid line: displacement u δ ; dashed line: damage α δ ; gray solid line: strain ε δ .

Figure 10 :

 10 Figure 10: Strain-Gradient model NS-G1:numerical results for the traction test. The simulation is without imperfection for η/L = 0.2 and n = 100 elements.

Figure 11 :

 11 Figure 11: Strain-Gradient model LS-G1: numerical results for the traction test. The simulation is without imperfection for η/L = 0.2 and n = 300 elements.

Figure 12 :

 12 Figure12: Strain-Gradient models: dissipated energies vs. mesh size at the end of the loading process. Markers are data from numerical simulations. The dashed and continuous lines serve as reference giving the dissipated energies corresponding to the full damage of few elements or of the whole bar.

Figure 13 :

 13 Figure13: Strain-Gradient models SG0: the results of the numerical simulations show uncontrolled localization of the damage ans strain fields and are extremely sensitive to the mesh, imperfections and solver parameters. The figures show a snapshot of the damage and strain field at t = 2.0 for η/L = 0.2 and n = 50 elements.

Figure B. 14 :

 14 Figure B.14: Geometrical interpretation of the energetic contents defining the function Fσ(ε). Upper left: (WH(ε) + σ 2 /(2E 0 )). Upper right: (σε). Bottom: the difference of the terms above, namely (E 0 Fσ(ε)/2).

2 )E 2 -L/ 2 E 2 E

 2222 (u, α)(v, β) E 0 = L/(α)v 2 + η 2 v 2 + 2E (α)u v β + 1 (α)u 2 β 2 dx. (C.3)

Table 1 :

 1 Résumé of the main properties for the damage models under consideration.

If one regards the phase-field model as a damage model with an internal length, the damage model is richer the Griffith model and one should regard the damage model as the true physical model and the Griffith model as an approximation. This view is opposite to that of Gamma-Convergence, where the damage model is regarded as a regularized approximation of the Griffith model.

https://bitbucket.org/cmaurini/gd-sd-damage-models-codes

We refer the reader to[START_REF] Pham | The variational approach to damage: I. The foundations[END_REF][START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF] for a time-continuous formulation of the evolution problem. Here we prefer to adopt the time-discrete setting to reduce the theoretical technicalities.

We voluntary refrain here from giving a precise definition of the functional setting. The question is subtle because the local damage model with softening is intrinsically ill-posed. Let us simply note that the condition of having a finite value of the energy functional functional[START_REF] Pham | The variational approach to damage: II. The gradient damage models [Approche variationnelle de l'endommagement: II. Les modèles à gradient[END_REF] does not prevent the damage field to be discontinuous and the displacement to jump wherever α = 1.

While it is physically reasonable to assume that the damage cannot increase the local and nonlocal stiffnesses, the case of a generic G(α) can be considered without additional difficulties.

See the in src/sgdam.py in https://bitbucket.org/cmaurini/gd-sd-damage-models-codes