Kim Pham

Hanen Amor

Jean-Jacques Marigo

Corrado Maurini
email: corrado.maurini@upmc.fr

Gradient damage models and their use to approximate brittle fracture

Keywords: Fracture, energy methods, non-local damage, stability, variational inequalities, finite elements

In its numerical implementation, the variational approach to brittle fracture approximates the crack evolution in an elastic solid through the use of gradient damage models. In the present paper, we first formulate the quasi-static evolution problem for a general class of such damage models. Then, we introduce a stability criterion in terms of the positivity of the second derivative of the total energy under the unilateral constraint induced by the irreversibility of damage. These concepts are applied in the particular setting of a one-dimensional traction test. That allows us to construct homogeneous as well as localized damage solutions in a closed form and to illustrate the concepts of loss of stability, of scale effects, of damage localization, and of structural failure. Considering several specific constitutive models, stress-displacement curves, stability diagrams, and energy dissipation provide identification criteria for the relevant material parameters, such as limit stress and internal length. Finally, the one-dimensional analytical results are compared with the numerical solution of the evolution problem in a two dimensional setting.

Introduction

Damage theory aims at modelling progressive degradation and failure in engineering materials such as metal, concrete, or rocks. Some of these materials exhibit a stress softening behavior that may lead to structural failure through damage localization in narrow bands. Damage localization may be interpreted as a regularized description of cracks, i.e. surfaces of discontinuities of the displacement field.

Local continuum models are not suitable to correctly predict damage localization processes. From the mathematical point of view, softening local damage models are ill-posed, showing lack of ellipticity and allowing for solutions with damage localization in band of null thickness with vanishing energy dissipation [START_REF] Benallal | Bifurcation and localization in rate independent materials[END_REF]. In finite element numerical simulations, the mesh size rules the thickness of the localization bands and thereby the energy stored inside the crack [START_REF] Bažant | Continuum theory for strain-softening[END_REF]. This fact implies an unacceptable dependence of the results on the discretization size. On the other hand, experimental evidences show that the strength of a structure considerably varies towards its size. For example, in metal composites with given volume fraction, the smaller is the size of a metal particle, the stronger is its response [START_REF] Lloyd | Particle reinforced aluminium and magnesium matrix composites[END_REF]. Again local approaches are unable to highlight these size effects since no characteristic length-scale is involved in their formulation. To avoid these pathological localizations and introduce scale effects, several regularization techniques have been developed, such integral (Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] or gradient [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF][START_REF] Comi | Computational modelling of gradient-enhanced damage in quasi-brittle materials[END_REF] approaches. Their common point is the enhancement of the description of the damage distribution around a material point by giving additional information on its "neighborhood". A consequence is the necessary introduction of an additional material parameter, the internal length. Despite of the large literature on the subject, the discussion on the effectiveness of the several possible regularization techniques is still open [START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF]. The main issue is that even the regularized non-local strain-softening models continue to be ill-posed, showing complex phenomena, such as infinite number of linearly independent solutions and continuous bifurcations [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF][START_REF] Pijaudier-Cabot | Strain localization and bifurcation in a non-local continuum[END_REF], which strongly depend on the specific constitute model under consideration.

Non-local damage models appear also as an elliptic approximation of the variational fracture mechanics problem. The variational approach of brittle fracture recast the evolution problem for the cracked state of a body as a minimality principle for an energy functional sum of the elastic energy and the energy dissipated to create the crack [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. On the basis of the results of the mathematical theory of image segmentation and free-discontinuity problems [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF], [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] approximate the minima of this energy functional through the minimization of a regularized elliptic functional that may be mechanically interpreted as the energy of a gradient damage model with an internal length. Mathematical results based on Gamma-Convergence theory show that when the internal length of a large class of gradient damage models tends to zero, the global minima of the damage energy functional tend towards the global minima of the energy functional of Griffith brittle fracture [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]. The same is true for the corresponding quasi-static evolutions ruled by a global minimality principle [START_REF] Giacomini | Ambrosio-Tortorelli approximation of quasi-static evolution of brittle fractures[END_REF]. However, the behavior of the regularized models with respect to local minima and their relation with the Griffith theory of fracture is less clear.

This work studies the solutions of non-local damage models with a gradient term on the scalar damage variable. It constitutes essentially an extension of the work of [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF]. One of its aim is to establish some properties with respect to a local minimality condition of the solutions of the regularized models used in variational theory of brittle fracture. It adopts a variational approach to equilibrium and stability based on an unilateral energy minimality principle under the irreversibility condition on the damage field [START_REF] Nguyen | Stability and Nonlinear Solid Mechanics[END_REF]. The framework is common to the work of Del [START_REF] Del Piero | Elastic bars with cohesive energy[END_REF] on the equilibrium and stability of elastic bar with cohesive energies.

The paper is organized as follows. Section 2 is devoted to the formulation of gradient damage models in a general context. In Section 3, we consider the homogeneous solutions of the evolution problem for the 1D traction of a bar and study their stability. We compare several models on the basis of two fundamental global properties: the force-displacement curve and the stability diagram for homogenous solutions. This analysis provides useful criteria for the selection of the constitutive laws and the identification of the internal length. In Section 4, we construct solutions with damage localization. Section 5 describes the numerical strategy that may be adopted to solve the damage evolution problem in a general setting. The numerical results obtained for the 2D version of the traction test are compared to the analytical ones derived in the 1D setting.

The summation convention on repeated indices is implicitly adopted in the sequel. The vectors and second order tensors are indicated by boldface letters, like u and σ for the displacement field and the stress field. Their components are denoted by italic letters, like u i and σ ij . The third or fourth order tensors as well as their components are indicated by a sans serif letter, like A or A ijkl for the stiffness tensor. Such tensors are considered as linear maps applying on vectors or second order tensors and the application is denoted without dots, like Aε whose ij-component is A ijkl ε kl . The inner product between two vectors or two tensors of the same order is indicated by a dot, like a • b which stands for a i b i or σ • ε for σ ij ε ij .

Variational formulation of isotropic damage models

Constitutive assumptions

We consider a gradient damage model in which the damage variable α is a real number growing from 0 to 1, α = 0 being the undamaged state and α = 1 being the fully damaged state. The behavior of the material is characterized by the state function W ℓ which gives the energy density at each point x. It depends on the local strain ε(u)(x) (u and ε(u) denotes respectively the displacement and the linearized strain), the local damage value α(x) and the local gradient ∇α(x) of the damage field at x. Specifically, we assume that W ℓ takes the following form

W ℓ (ε(u), α, ∇α) := 1 2 A(α)ε(u) • ε(u) + w(α) + 1 2 w 1 ℓ 2 ∇α • ∇α (1)
where w 1 = w(1) represents the specific fracture energy of the brittle damage material defined in [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF] (see [START_REF] Ramtani | Orthotropic behavior of concrete with directional aspects: modelling and experiments[END_REF] for its experimental determination), A(α) the rigidity of the material in the damage state α. The function w(α) can be interpreted as the density of the energy dissipated by the material during a homogeneous damage process (such that ∇α(x) = 0), where the damage variable of the material point grows from 0 to α.

The last term in the right hand side of (1) is the non local part of the energy which plays a regularizing role by limiting the possibilities of localization of the damage field. For reasons of physical dimension, it involves a material characteristic length ℓ that will control the thickness of the damage localization zones. The strain energy density of the underlying local model is obtained by setting ℓ = 0. The local model with the strain energy W 0 (ε(u), α) is not well suited for the modeling of softening material behavior, because, as explicitly shown on various examples in Pham and Marigo (2009b), localization can be obtained without any cost of energy.

The qualitative properties of the damage model strongly depend on the properties of the stiffness function α → A(α), the dissipation function α → w(α), the compliance function α → S(α) = A -1 (α) and their derivatives. From now on, we will adopt the following hypotheses that characterize the material behavior:

Hypothesis 1 (Constitutive assumptions). The functions α → A(α), α → S(α) = A -1 (α) and α → w(α) are continuously differentiable on [0, 1) and fulfill the following conditions Positive elasticity : A(α) > 0, S(α) > 0 with A(1) = 0.

(2a)

Decreasing stiffness : A ′ (α) < 0, S ′ (α) > 0. (2b) Dissipation : w(0) = 0, w(α = 0) > 0, w ′ (α) ≥ 0. ( 2c
)
Remark 1. The present model is a brittle damage model without any residual strains and plasticity-like effects. Including plastic behaviors in our formulation is not a major difficulty. It would require to change the form of the strain energy by replacing the total strains by the elastic strains in the elastic energy part, and by adding a plastic dissipated energy (possibly damage dependent). More generally, as far as rate independent behaviors are concerned, one can develop a variational approach like the present one. Such formulations are a generalization of the concept of Standard Materials, see [START_REF] Nguyen | Stability and Nonlinear Solid Mechanics[END_REF], [START_REF] Mielke | Evolution of rate-independent systems[END_REF] and Pham and Marigo (2010a). On the other hand, the modeling of rate dependent behaviors by such a variational approach is less developed, because of fundamental theoretical difficulties. Specifically, the difference between rate independent and rate dependent behaviors is that the dissipated energy (or more precisely, the dissipation potential in the language of Generalized Standard Materials) is homogeneous of degree 1 in the former case, but not in the latter. This lack of one homogeneity induces serious difficulties from a theoretical viewpoint, cf [START_REF] Bourdin | The variational approach to fracture[END_REF].

Energy functional

We consider a homogeneous body Ω ⊂ R n made of the non-local damaging material characterized by the state function (1). Let us denote respectively by u t , F t , and f t the imposed displacement on the boundary part ∂Ω U , the surface forces on the complementary boundary part ∂Ω F , and the volume forces over Ω. We define the space of kinematically admissible displacement fields at time t as:

C(U t ) = v ∈ (H 1 (Ω)) n : v = u t on ∂ U Ω , (3)
H 1 (Ω) being the usual Sobolev space of functions defined on Ω which are square integrable and have square integrable first derivative. The set of admissible damage fields is the following convex subset of the Sobolev space H 1 (Ω)

D 1 = β ∈ H 1 (Ω) : β(x) ∈ [0, 1] for almost all x (4)
Then for any admissible pair (u, α) in C(U t ) × D 1 at time t, we define the total energy of the bar

P t (u, α) := E(u, α) - Ω f t • u dΩ - ∂ F Ω F t • u dΓ (5)
where E(u, α) stands for the total strain energy (see Clapeyron's theorem below for a justification of this terminology)

E(u, α) = Ω W ℓ (ε(u)(x), α(x), ∇α(x)) dΩ. ( 6
)
Remark 2. In the definition (4) of the admissible damage fields, the damage is left free on the boundary of Ω, without specifying any boundary condition of the Dirichlet type. The complementary Neumann boundary conditions will automatically appear through the variational approach, as it is shown in the following equation (12d). Other types of boundary conditions are possible, but, for sake of conciseness, this is the only one considered in the analytical study. The numerical simulations of Section 4 will illustrate the case where the damage is set to zero on the boundary.

2.3. First-order optimality conditions: evolution problem, equilibrium, and damage criterion

In the variational approach, the quasi-static evolution problem for the displacement and the damage fields is formulated as a first-order unilateral minimality condition on the functional (5) under the irreversibility condition that the damage can only increase, that is α ≥ 0 (here and henceforth the superimposed dot denotes the derivative with respect to t). Considering an initial undamaged state at t = 0, that leads to the following variational evolution problem (see [START_REF] Marigo | Constitutive relations in plasticity, damage and fracture mechanics based on a work property[END_REF]; [START_REF] Nguyen | Bifurcation and postbifurcation analysis in plasticity and brittle fracture[END_REF]; [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF]):

For each t > 0, find (u t , α t ) in C(u t ) × D 1 such that ( ut , αt ) ∈ C( ut ) × D and for all (v, β) ∈ C( Ut ) × D, DP(u t , α t )(v -ut , β -αt ) ≥ 0, (7)
with the initial condition α 0 (x) = 0. In (7), D denotes the convex cone of positive damage rate

D = β ∈ H 1 (Ω) : β(x) ≥ 0 for almost all x (8)
whereas DP(u, α)(v, β) denotes the Gâteaux derivative of P at (u, α) in the direction (v, β):

DP(u, α)(v, β) = Ω A(α)ε(u) • ε(u) dΩ + + Ω 1 2 A ′ (α)ε(u) • ε(u) + w ′ (α) βdΩ + Ω w 1 ℓ 2 ∇α • ∇β dΩ. (9)
From this global variational formulation, we can deduce the standard local formulation of the damage model by integration by parts and classical localization arguments.

The equilibrium equation for the stress σ t = A(α)ε(u t ) is obtained by testing the variational inequality (7) for β = αt and v ∈ C( ut ). This gives the standard equilibrium equations

div σ t + f t = 0 in Ω, σ t • n = F t on ∂ F Ω. ( 10
)
The damage problem is obtained testing (7) for arbitrary β in the convex cone D with v = ut . That leads to the variational inequality governing the evolution of the damage

Ω w 1 ℓ 2 ∇α t • ∇β dΩ + Ω 1 2 A ′ (α t )ε(u t ) • ε(u t )β dΩ + Ω w ′ (α t )β dΩ ≥ 0 (11)
which must hold for all β ∈ D, the inequality becoming an equality when β = αt . After an integration by parts and using classical tools of the calculus of variations, we find the strong formulation for the damage evolution problem in the form of the Kuhn-Tucker conditions for unilateral constrained variational problems: The two local problems (10) and ( 12) are implicitly linked and must be satisfied simultaneously to get a solution of the evolution problem (7).

Irreversibility : αt ≥ 0 on Ω (12a) Damage criterion : 1 2 A ′ (α t )ε(u t ) • ε(u t ) + w ′ (α t ) -w 1 ℓ 2 ∆α t ≥ 0 on Ω (12b) Energy balance : αt 1 2 A ′ (α t )ε(u t ) • ε(u t ) + w ′ (α t ) -w 1 ℓ 2 ∆α t =

The stability criterion

In softening materials, when the yield criterion of damage is reached somewhere in the body, a multiplicity of solutions can appear for the same loading history [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF]Pham and Marigo, 2009a). This possible loss of uniqueness, which depends in particular on the size of the body, requires a criterion to select the stable states. We adopt here an energetic point of view by defining the stable states at a given time t as those which are Unilateral Local Minima of the total energy P t . Specifically, for a given couple (u, α) of admissible fields in C t × D 1 , the stability condition can be written formally as

∃r > 0, ∀(v, β) ∈ C(0) × D : (v, β) = 1, ∀h ∈ [0, r] P t (u, α) ≤ P t (u + hv, α + hβ) (13)
where • stands for the H 1 norm.

Let us emphasize that the stability of a state is only tested in the positive direction β of damage. This restriction is due to the irreversibility of damage. Indeed, even if there exists in the neighborhood of the tested state (u, α) another state (u * , α * ) with a smaller energy but with α * (x) < α(x) somewhere in the body, because of the irreversibility of the damage, it is not possible that the body evolves to this state, see also [START_REF] Nguyen | Bifurcation and postbifurcation analysis in plasticity and brittle fracture[END_REF][START_REF] Nguyen | Stability and Nonlinear Solid Mechanics[END_REF].

In practice, we test the stability of any state (u, α) (not necessarily associated with an evolution path) by doing a Taylor expansion of the total energy (26) for a given loading around this state in the admissible direction (v,

β) ∈ C(0) × D P(u + hv, α + hβ) = P(u, α) + h DP(u, α)(v, β) + h 2 2 D 2 P(u, α)(v, β) + o(h 2 ). ( 14
)
Suppose that the displacement field u is the one satisfying the equilibrium (10) at the damage state α. Then, we can then distinguish two cases:

(i) if the damage criterion (12b) at the state (u, α) is a strict inequality everywhere in the structure, then, for any non-zero damage field direction β ∈ D, the first derivative ( 9) is strictly positive and the state (u, α) is stable;

(ii) if there exists a non-zero measure subset Ω D (α) ⊂ Ω over which the damage criterion (12b) is an equality, then the first derivative of the energy is not sufficient to conclude and the stability of the state is given by the sign of the second derivative of the energy over the convex subset D(α) of D,

D(α) = {β ∈ D, Supp(β) ⊂ Ω D (α)} , (15)
see [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF], Pham and Marigo (2009b), Pham and Marigo (2010b).

Strain hardening and stress softening

For a homogeneous damage distribution α, let us respectively define the elastic domain in the strain and stress spaces as

A(α) = ε ∈ Sym, -A ′ (α)ε • ε ≤ 2w ′ (α) , (16)
A * (α) = σ ∈ Sym, S ′ (α)σ • σ ≤ 2w ′ (α) , ( 17
)
where Sym is the space of symmetric tensors. Depending on the behavior of the quadratic forms defining these domains, the material is said to be -strain hardening when α → (-A ′ (α)/w ′ (α)) is decreasing with respect to α, i.e.:

w ′ (α)A ′′ (α) -w ′′ (α)A ′ (α) > 0. ( 18
)
This means that the domain of admissible strains in the elastic regime is increasing for increasing damage.

-stress hardening (resp. softening) when α → (S ′ (α)/w ′ (α)) is decreasing (resp. increasing) with respect to α, i.e.:

w ′ (α)S ′′ (α) -w ′′ (α)S ′ (α) < (resp. >)0, (19)
then the domain of admissible stress in the elastic regime is increasing (resp. decreasing) for increasing damage.

The importance of these definitions will appear in the study of the solutions of the variational problem.

A Clapeyron's theorem for the evolution problem

In linear elasticity, Clapeyron's theorem states that the total strain energy of a body in equilibrium under a prescribed loading, is equal to one-half of the work done by the external forces in the displacement solution field. In the case of our damage problem, we can formulate an analog Clapeyron's theorem. Indeed we have the following property.

Property 1. Let (u t , α t ) be a displacement-damage path parametrized by a increasing loading parameter t ∈ [0, 1] with (u 0 , α 0 ) = (0, 0) and (u 1 , α 1 ) = (ũ, α), and satisfying the evolution problem (7) at each time. Then the total strain energy of its final state (ũ, α) is given by

E(ũ, α) = Ω 1 2 w 1 ℓ 2 ∇α • ∇α dΩ + Ω 1 2 A(α)ε(ũ) • ε(ũ) dΩ + Ω w(α) dΩ = 1 0 P e ( ut )dt, ( 20
)
where P e ( ut ) is the power of the external forces in the displacement rate ut

P e ( ut ) = Ω f t • ut dΩ + ∂ F Ω F t • ut dΓ + ∂ U Ω σ t n • ud t dΓ. ( 21
)
This means that the total energy at time t is equal to the work of the external forces up to t. Moreover, since the work of the external forces is equal to the strain work up to time t, this latter one is a state function and equal to the total strain energy:

1 0 Ω σ t • ε( ut ) dΩ dt = E(ũ, α). ( 22
)
The proof is detailed in the Appendix.

Approximation of variational brittle fracture

In fracture mechanics approaches, material failure is modeled by the nucleation and propagation of surfaces of discontinuity of the displacement field. Considering a solid Ω ⊂ R n , the cracks, say Γ, are surfaces of dimension n -1. The quasi-static problem of fracture mechanics requires to determine the evolution of the displacement field u and the crack set Γ as a function of the loading. The variational approach to brittle fracture formulates this evolution problem as a global minimization condition on the following Griffith energy functional, under an irreversibility condition for the crack set Γ [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]:

F(u, Γ) = Ω\Γ 1 2 A 0 ε(u) • ε(u)dΩ + G c H n-1 (Γ) ( 23
)
The energy is the sum of the elastic energy stored in the cracked body Ω \ Γ and the surface energy required to create the crack, H n-1 being the surface measure of Γ, i.e. the total crack length for n = 2 or the total crack surface for n = 3. The constant G c is the material toughness, i.e. the energy required to create a unit crack surface according to Griffith theory.

The functional ( 23) has a close analogy with the damage energy functional in the form (6). Both of them are the sum of an elastic energy term and a dissipated energy term, the dissipated energy being a volume integral for the damage model and a surface integral for the fracture model. Indeed, the damage functional is a regularized version of the Griffith functional. Gamma-Convergence results [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF] show that, under some constitutive requirements, the global minimum of the damage functional (6) converges toward the global minimum of the Griffith functional (23) when the internal length ℓ goes to zero, provided that the dissipated energy function w(α) is related to the toughness G c and the internal length ℓ by

G c = 2ℓ √ 2 1 0 w 1 w(β)dβ. ( 24
)
3. Application to the 1D tension test

The 1D problem

In the following, we illustrate the solution of the evolution problem for the traction of a one-dimensional bar made of a homogeneous material with stress softening. The end x = 0 is blocked while the end x = L has a displacement U t = t L imposed by a hard device (see Figure 1). Hence, the admissible displacement field u must satisfy the boundary conditions

u(x = 0) = 0, u(x = L) = U t , U 0 = 0. ( 25
)
U t The total energy (5) reads now as

P(u, α) = L 0 1 2 E(α)u ′2 dx + w(α) + 1 2 w 1 ℓ 2 α ′2 dx, (26)
where E(α) is the one-dimensional axial stiffness and (•) ′ = ∂(•)/∂x. The equilibrium equation ( 10) reads as:

σ ′ t (x) = 0, σ t (x) = E(α t (x))u ′ t (x), ∀x ∈ (0, L). (27)
Then, the stress σ t along the bar is necessarily constant. Using the boundary conditions ( 25), we get

σ t = U t L 0 S(α t (x))dx , (28)
with S(α) = 1/E(α). Then the damage problem (12) becomes

Irreversibility : αt ≥ 0, α 0 = 0, (29a)
Damage criterion :

-w 1 ℓ 2 α ′′ t + 1 2 E ′ (α t )u ′ t 2 + w ′ (α t ) ≥ 0, ( 29b
)
Energy balance : αt -w 1 ℓ 2 α ′′ t + 1 2 E ′ (α t )u ′ t 2 + w ′ (α t ) = 0, (29c)
Boundary conditions :

α ′ t (0) ≥ 0, α ′ t (L) ≥ 0. ( 29d
)
Depending on the ratio L/ℓ between the length of the bar and the internal length, this problem can admit a unique solution or, on the contrary, an infinite number of solutions. In the latter case, the different solutions depend on the zone where the damage criterion (29b) and the boundary conditions (29d) become an equality. The following sections are devoted to the analysis of the problem (29).

Homogeneous solutions

We consider only the case of strain hardening materials, in the sense of equation ( 18). We focus here on the study of homogeneous solutions (tx, α t ) of the evolution problem, i.e. solutions for which the damage and the strain fields are uniform throughout the structure. Hence, for the given boundary conditions, the solution of the elastic problem (27) takes the form:

u ′ t (x) = t, u t (x) = t x, σ t = t E(α t ). ( 30
)
The damage criterion is also uniform in the bar and reads as

t 2 2 E ′ (α t ) + w ′ (α t ) ≥ 0, αt t 2 2 E ′ (α t ) + w ′ (α t ) = 0. ( 31
)
Since the bar is undamaged at the beginning of the loading (α 0 = 0), the solution of the evolution problem is characterized by an elastic phase and a damaging phase, as detailed below.

Elastic phase

For t increasing from 0, the damage criterion (29b) is a strict inequality if the end displacement U t is smaller than the elastic limit

U e = L - 2w ′ (0) E ′ (0) , with σ e = L 2w ′ (0) S ′ (0) , (32)
The stress value σ e can be interpreted as the yield stress. If w ′ (0) > 0, then U e > 0 and the material has a genuine elastic phase: as long as U t lies in the interval [0, U e ], the material remains sound (α t = 0) with a stiffness E 0 = E(0) and the stress is given by σ t = E 0 U t /L. On the other hand, if w ′ (0) = 0, then the damage criterion becomes an equality at the onset of the loading and we cannot observe an elastic phase.

Damaging phase

For U t ≥ U e the damage criterion (31) becomes an equality and the damage can grow. Hence, (31) gives the following implicit relation between the prescribed displacement U t and the associated homogeneous damage α t :

U t L = - 2w ′ (α t ) E ′ (α t ) . ( 33
)
Owing to the strain hardening assumption, (-w ′ /E ′ ) is a monotonic increasing function. Thus, there is a unique solution α t for a given U t . In other words, the strain hardening hypothesis rules out snap-back phenomena during the evolution in time of homogeneous solutions. The corresponding stress is

σ t = E(α t )t = 2w ′ (α t ) S ′ (α t ) . ( 34
)
Equation ( 34) clearly shows that the stress decrease (resp. grows) with the displacement if the material is stress softening (resp. stress hardening) at the damage level α t .

Remark 3 (Peak Stress). From the constitutive assumptions (2), w ′ and S ′ are positive continuous functions on α ∈ [0, 1). Hence, we can define the peak stress in homogeneous solutions as

σ M = sup α∈[0,1) 2w ′ (α) S ′ (α) . ( 35
)
When σ M < +∞, then σ M is the maximal stress that the material can sustain. In particular, in the case of stress softening laws, the peak stress ( 35) is attained for α = 0.

Stability and role of the internal length

Although the homogeneous strain-damage solution (tx, α t ) t≥0 verifies the set of the equilibrium equations and the damage criterion, it could be not observable during a tensile test. Indeed, some states (tx, α t ) may be unstable. A stability analysis is then necessary to select interesting solutions from the physical point of view. It consists in studying the sign of the first and second derivative of the total energy at the state (tx, α t ) (see Section 2.4). These derivatives read as:

DP(tx, α t )(v, β) = L 0 1 2 E ′ (α t )t 2 + w ′ (α t ) βdx (36a) D 2 P(tx, α t )(v, β) = L 0 w 1 ℓ 2 β ′2 + 2E ′ (α t )tv ′ β + E ′ (αt) 2 v ′2 + w ′′ (α t ) β 2 dx. (36b)
As discussed in section 2.4, the states of the elastic phase are stable. Indeed, the inequality is strict in (31-a) and since the first derivative of the energy is positive in each direction β = 0, these states are stable. Beyond the elastic phase, i.e. for U t ≥ U e , the damage criterion (31-a) becomes an equality and the first derivative of the energy vanishes. The second derivative reads as

D 2 P(tx, α t )(v, β) = L 0 w 1 ℓ 2 β ′2 dx+ L 0 E(α t ) v ′ + E ′ (α t ) E(α t ) tβ 2 dx - L 0 1 2 S ′′ (α t )σ 2 t -w ′′ (α) β 2 dx. ( 37
)
The first two terms are positive. The last term of (37) is positive if and only if the behavior is with stress hardening (see ( 19)). Indeed, combining (34) with ( 19) leads to

0 < w ′′ (α t )S ′ (α t ) -w ′ (α t )S ′′ (α t ) = -S ′ (α t ) 1 2 S ′′ (α t )σ 2 t -w ′′ (α t ) , (38)
where S ′ (α t ) > 0 by virtue of (2b). In this case, we deduce that the state is necessarily stable.

On the other hand, if the material is stress softening, the last term of (37) is negative and may induce an instability. In this case, the sign of the second derivative may be assessed through the study of the following Rayleigh ratio

R t (v, β) = L 0 w 1 ℓ 2 β ′2 dx + L 0 E(α t ) v ′ + E ′ (αt) E(αt) tβ 2 dx L 0 1 2 S ′′ (α t )σ 2 t -w ′′ (α t ) β 2 dx . ( 39
)
Specifically, a sufficient (resp. necessary) condition for stability is that

ρ = min (v,β)∈C 0 ×D R t (v, β) > (resp. ≥)1. ( 40
)
After some calculations which are not reproduced here, the infimum of the Rayleigh ratio ( 39) is given by

ρ = min E(α t )S ′ (α t ) 2 σ 2 t ; π 2 w 1 ℓ 2 L 2 E(α t ) 2 S ′ (α t ) 4 σ 4 t 1/3 1 2 S ′′ (α t )σ 2 t -w ′′ (α t ) . ( 41
)
We may finally summarize the results on the stability of homogeneous states in the following proposition.

Property 2 (Stability of homogeneous states). In the elastic phase, U t < U e , the homogeneous state of strain-damage (U t x/L, 0) is stable. For U t ≥ U e , if the material has a stress hardening behavior (19), then the state is stable. In the case of stress softening behavior (19), the state (U t x/L, α t ) with α t given by (33) is stable if and only if the length of the bar L satisfies the inequality

L 2 ℓ 2 ≤ π 2 w 1 E(α t ) 2 S ′ (α t ) 4 σ 4 t 1 2 S ′′ (α t )σ 2 t -w ′′ (α t ) 3 , (42)
where σ t = E(α t )U t /L is the stress at the equilibrium.

Remark 4. The Property 2 points out a size effect in the stability result due to the presence of the internal length. Indeed, the value α t of the homogeneous damage given by (33) does not depend on the length of the bar. Accordingly, small or large bars give the same stress-strain diagram for homogeneous responses. However, according to (42), for large ratio L/ℓ, the homogeneous state will be unstable and will not be observable during an experiment. On the contrary, small bars allow homogeneous stable states.

To illustrate the previous stability analysis, we consider three examples of damage laws.

Example 1: a model with an elastic phase

We consider the following damage law

E(α) = E 0 (1 -α) 2 , w(α) = w 1 α, (43)
This law satisfies the Hypothesis (2a)-(2c) and both the strain hardening (18) and stress softening (19) conditions for any α. Therefore the evolution problem admits a unique homogeneous strain-damage solution. Since w ′ (0) > 0, there exists an elastic phase in the evolution problem (Section 3.2.1). Hence, using ( 32) and ( 35), the displacement and the stress at the elastic limit read

σ e = σ M = w 1 E 0 , U e = w 1 E 0 L = σ M E 0 L. (44)
Using ( 33) with ( 43) we deduce the value of the damage at each time

α t = max 0, 1 - U e U t 2 . ( 45
)
The associated stress-strain relation reads

σ t =        σ M U t U e if U t ≤ U e , σ M U e U t 3 otherwise. ( 46
)
We notice that the peak stress σ M and the yield stress σ e given by ( 32) are identical for this law. Figure 2(a) reports the stress in the bar versus the normalized end displacement U t /U e , which is proportional to the homogeneous strain t = U t /L. After the elastic phase, the stress decreases asymptotically to 0, a distinctive feature of stress softening. Let us now discuss the stability of the homogeneous solutions. For U t < U e the state (tx, 0) is stable. For U t ≥ U e , inserting (45)-( 46) into (42), we deduce that the homogeneous state (tx, α t ) is stable if and only if the length of the bar L satisfies the following condition:

0 1 2 0 1 U t U e Σ t Σ M (a) Stress vs strain response 0 1 2 3 0 Λ c U t U e L Stable Unstable (b) Stability diagram
L ℓ ≤ λ c U e U t , with λ c = 4π 3 √ 3 . ( 47
)
Figure 2(b) resumes these conditions on a stability diagram in the U t /U e -L/ℓ plane. Depending on the length of the bar, two different behaviors occur:

• For L > λ c ℓ, the homogeneous state is unstable for any U t ≥ U e and a damage localization necessarily arises at the end of the elastic phase;

• For L < λ c ℓ, after the elastic phase (U t > U e ), there exists a finite interval [U e , λ c U e ℓ/L) of the prescribed displacement for which the homogeneous state is still stable. However, the homogeneous state becomes unstable for U t ≥ λ c U e ℓ/L.

Example 2: a model without an elastic phase

We consider the following damage law

E(α) = E 0 (1 -α) 2 , w(α) = w 1 α 2 , (48)
This law satisfies the Hypothesis (2a)-(2c) and the strain hardening condition (18) for any α. However the stress softening condition ( 19) is ensured only for α ≥ 1/4. Moreover since w ′ (0) = 0, this damage law does not have an elastic phase (σ e = 0). Hence, using (35), the peak stress is

σ M = 3 √ 3 8 √ 2 w 1 E 0 . (49)
From ( 33), we deduce that the evolution of the homogeneous damage level is given by

α t = U 2 t U 2 t + 3U 2 M with U M = 16σ M 9E 0 L. ( 50
)
Therefore the relation between the stress of the homogeneous solution and the prescribed displacement at x = L is given by

σ t = E 0 9U 4 M (U 2 t + 3U 2 M ) 2 U t L . (51)
The material becomes stress softening for U t ≥ U M . The peak stress is σ M and is reached for U t = U M . Figure 3(a) shows the stress in the bar versus the homogeneous strain.

In this case, the stability of the state is ruled by second derivative of the energy. Indeed, there is not an elastic phase, and the first derivative of the potential energy is always zero. During the stress hardening phase (U t < U M ), the homogeneous state is stable no matter the length of the bar. Beyond U M , the bar is in a stress softening regime and we can apply Property 2: the state is stable if the length of the bar verifies the following condition

0 1 2 3 0 1 U t U M Σ t Σ M (a) Stress vs strain response 0 1 2 3 0 10 U t U M L Stable Unstable (b) Stability diagram
L ℓ ≤ π √ 3 4 U 2 t /U 2 M (U 2 t /U 2 M -1) 3/2 . ( 52
)
The stability diagram is plotted on Figure 3(b). Beyond the stress hardening regime (U t ≥ U M ), whatever the choice of the length ratio L/ℓ, size effects rule the stability of the solution. For sufficient large displacements (the critical value depending on the bar length), the homogeneous state becomes unstable and a localization arises somewhere in the bar.

Example 3: a family of models with the same homogeneous strain-stress response

We consider the following family of damage models indexed by the parameter p > 0:

E(α) = E 0 (1 -α) p , w(α) = w 1 1 -(1 -α) p/2 . ( 53
)
This case is a generalization of the law (43) which is recovered for p = 2. It satisfies the Hypothesis (2a)-(2c) and both the strain hardening (18) and stress softening (19) conditions for any α ∈ [0, 1) and any p > 0. Since w ′ (0) > 0, the damage evolution contains an elastic phase (Section 3.2.1) and there exists a unique homogeneous strain-damage solution. Hence, using (32) and ( 35), the displacement and the stress at the elastic limit read Using ( 33) with ( 53) we deduce the value of the damage at each time

σ e = σ M = w 1 E 0 , U e = w 1 E 0 L = σ M E 0 L. (54)
α t = max 0, 1 - U e U t 4/p . (55)
The corresponding stress is

σ t =        σ M U t U e if U t ≤ U e , σ M U e U t 3
otherwise.

(56)

Therefore, this family of damage laws, indexed by p, leads to the same strain-stress diagram whatever the value of p.

This example points out that a force-displacement diagram obtained from a tensile test is not sufficient to identify the damage model. Indeed, the force-displacement response only involves one combination of the two state functions (E(α) and w(α)). At least another experimental curve is necessary. We claim that the stability analysis provides the missing information. Indeed, even if the homogeneous (elastic) states are stable for U < U e , for U ≥ U e the stability of an homogeneous state depends both on the length of the bar and on p. Specifically, the homogeneous state (tx, α t ) is stable when U t ≥ U e if and only if

L ℓ ≤ λ (p) c U e U t 4/p-1 , with λ (p) c = 8π 3 √ 3p . ( 57
)
Therefore the stability diagram depends on p, cf Figure 4:

• When p < 4, for sufficient large displacements (the critical value depending on the bar length) the homogeneous state becomes unstable.

• When p > 4, if L > λ p c ℓ, all homogeneous states (U t x/L, α t ) with U t in the interval [U e , U (p) c ) being U (p) c = λ (p) c ℓ L p 4 -p U e .
are unstable. However, the homogeneous states are stable again for

U t ≥ U (p) c .
This example underlines the relevance of the stability properties in the identification the damage laws E(α) and w(α).

Localized solutions

When a homogeneous state looses its stability, the damage field localizes on zones whose size is controlled by the internal length ℓ. We study here the properties of solutions with damage localized by focusing on the class of stress-softening material with an elastic phase. We only indicate the main lines of the construction of the localized solutions. The interested reader should refer to Pham and Marigo (2009a) for the proofs and the details of the calculations.

Optimal damage profile

Let σ ∈ (0, σ M ) be the equilibrium stress of the bar and S = (x 0 -D, x 0 + D) be a putative localization zone, where the thickness D has to be determined and x 0 is an arbitrary point of the bar. To construct the damage profile in the localization zone, we suppose that the damage criterion ( 29b) is an equality only on S and that the damage is zero on the remaining part of the bar. Accordingly, the damage field α satisfies

-σ 2 S ′ (α) + 2w ′ (α) -2w 1 ℓ 2 α ′′ = 0 on S, α = 0 on (0, L)/S. ( 58
)
Since α and α ′ must be continuous1 at x 0 ± D, we have

α(x 0 ± D) = α ′ (x 0 ± D) = 0. ( 59
)
Multiplying ( 58) by α ′ and integrating with respect to x, we obtain the first integral where C is a constant. Evaluating (60) at x 0 ± D, the conditions ( 59) and Hypothesis 1 give C = -σ 2 /E 0 . Hence the first integral may be written in the form

-σ 2 S(α) + 2w(α) -w 1 ℓ 2 α ′2 = C in S, ( 60
)
0 x 0 2 x 0 x 0 2 L 0 1 U t U e Α
ℓ 2 α ′ (x) 2 = H(σ, α(x)) in S, (61)
where

H(σ, β) := 2 w(β) w 1 - σ 2 w 1 E 0 (E 0 S(β) -1) , with β ∈ [0, 1). ( 62
)
The localized damage profile is the cusp-shaped curve reported in Figure 5. We define ᾱ(σ) as the maximal value of the damage along the bar at the given stress value σ, which is attained at x 0 , the center of the localization zone. This maximal value depends only on σ. When σ decreases from σ M to 0, ᾱ(σ) increases from 0 to 1.

The size D 0 of the localization zone is deduced from (61) by integration. It is a function of the stress σ:

D(σ) = ℓ ᾱ(σ) 0 dβ H(σ, β) . ( 63
)
For the assumed constitutive behavior, the integral above is well defined and always exists.

The position x 0 of the center of the localized solution can be chosen arbitrarily in the interval [D(σ), L -D(σ)]. We finally deduce from (61) that, in the localization zone, the damage field is given by the following implicit relation between x and α:

|x -x 0 | = ℓ ᾱ(σ) α dβ H(σ, β) . ( 64
)
It is easy to see that the damage field is symmetric with respect to the center of the localization zone, decreasing continuously from ᾱ(σ) at the center to 0 at the boundary (see Figure 5).

Energy of localized solutions

The total strain energy of the bar is given by

E = L 0 1 2 w 1 ℓ 2 α ′2 + 1 2 S(α)σ 2 + w(α) dx. ( 65
)
Assuming that a single localization zone exists and that it is surrounded by an elastic zone, we deduce from the first integral (60) that the total strain energy reads as

E = 1 2 σ 2 E 0 L + L 0 2w(α(x)) dx. ( 66
)
Since α is equal to 0 on (0, L)/S and w(0) = 0, the integral in ( 66) can be reduced to the subdomain S. Finally, using the symmetry of the profile α over S and performing the change of variable y = α(x), we obtain

E = 1 2 σ 2 E 0 L + 4ℓ ᾱ(σ) 0 w(β)dβ H(σ, β) . (67)

Cracks as localized damage zones with zero stress

In our damage context, a crack can be viewed as a damage localization whose corresponding equilibrium stress σ is 0. In this limit case, we deduce from previous analysis that

ᾱ(0) = 1, ℓ 2 α ′ (x) 2 = 2 w(α) w 1 , D 0 := D(0) = ℓ 1 0 w 1 2 w(β) dβ. (68)
As far as the total strain energy (67) is concerned, the first term vanishes, whereas the second gives

G c = 2 ℓ 1 0 2 w 1 w(β) dβ. (69)
This energy corresponds to the energy dissipated in the creation of a single crack and can be identified with the fracture toughness of the material.

Fracture, damage and internal length

The analytical results of this section show that the properties of the homogeneous and localized solutions of the gradient damage models with stress softening are characterized by several parameters, as E 0 , σ M , G c , D 0 , and w 1 . These parameters may be regarded as material constants to identify by carefully designed experimental procedures. The Young modulus E 0 and the peak stress σ M are easily obtained through classical uniaxial tensile tests. The fracture toughness G c is another classical parameter found in material databases. In a uniaxial traction test, it may be calculated as the ratio of the total dissipated energy when breaking the bar in two parts and the bar cross sectional area. With field measurements (displacement, thermal) and analysis of cracked specimens, it is also possible to estimate the thickness 2D 0 of the localization zone around the crack. The specific fracture energy w 1 corresponds to the total dissipated energy to damage up to failure a bar in a homogeneous traction test2 , divided by the volume of the specimen. Differently from the fracture toughness G c , which is a measure of the energy dissipated in a localized fracture mode, w 1 is a measure of the energy dissipated in homogeneous damage processes.

On the other hand, the gradient damage models introduced in Section 2 are characterized by two constitutive functions, the damaged stiffness α → E(α) and the density of dissipated energy α → w(α), and an additional material parameter, the internal length ℓ.

The present analytical study of the homogeneous and localized solutions of the gradient damage models rigorously establishes several explicit relations between the material and the model parameters. Because of their relevance, they are resumed below (recall also that E 0 = E(0) and w 1 = w(1)):

σ M = sup α∈[0,1) 2w ′ (α) S ′ (α) ; (70a)
G c = 2 ℓ 1 0 2 w 1 w(β) dβ; (70b)
D 0 = ℓ 1 0 w 1 2 w(β) dβ. (70c)
The equations above include a precise relation between the fracture toughness G c defined in brittle fracture models [START_REF] Bourdin | The variational approach to fracture[END_REF]), the peak stress σ M , the specific fracture energy w 1 , and the internal length defined in damage models [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF]). In this sense, they establish an explicit correspondence between brittle fracture and localized solutions of damage models. These properties are closely related to the Gamma-convergence results available for regularized approaches to brittle fracture à la Ambrosio-Tortorelli (see [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF].

Example

We consider the model characterized by the constitutive functions (43). A straight calculus of ( 68)-( 69) gives

G c = 4 √ 2 3 w 1 ℓ, D 0 = √ 2ℓ, σ M = w 1 E 0 . (71)
Therefore, the knowledge of the fracture toughness G c , the yield elastic stress σ M , and the Young modulus E 0 , allows us to compute the value of the internal length ℓ (for given constitutive functions). In addition, the damage profile in the zone where the bar is broken is

α(x) = 1 - |x -x 0 | √ 2ℓ 2 , x ∈ [x 0 -D 0 , x 0 + D 0 ], (72)
α being zero elsewhere (see Figure 5).

For concrete, a brittle material, the typical values for the material parameters are (see e.g. [START_REF] Bazant | Measurement of characteristic length of nonlocal continuum[END_REF] and [START_REF] Comi | Fracture energy based bi-dissipative damage model for concrete[END_REF]):

E 0 = 29 GPa, σ M = 4.5 MPa, G c = 70 N/m. ( 73
)
Using the relations (71), we find 2D 0 = 106 mm, ℓ = 38 mm,

w 1 = 698 N/m 3 . ( 74
)
The value of the thickness crack band in uniaxial tension test is fairly reasonable and can be compared to the phenomenological law used in the literature as a function of average aggregate size d a . As a rule of thumb, it is commonly considered that a typical value for 2D 0 is approximatively 3d a , where, for concrete, d a ∼ 10 ÷ 40mm.

Numerical implementation and examples

This section reports on the numerical solution of the damage evolution problem (7) for the 2D version of the uniaxial traction problem. We consider the damage model of Example 1 in Section 3.2.4, where:

A(α) = E 0 (1 -α) 2 , w(α) = w 1 α. (75)
There are two main reasons to select this model: (i) the model has a non-vanishing elastic phase; (ii) the associated functional P(u, α) turns out to be quadratic in the damage variable, a valuable feature for its numerical minimization.

Time-discrete evolution and solution algorithm

Section 2.3 formulates the variational evolution problem with a continuous time variable t. For the numerical work, the time line is discretized in a set of N + 1 time instants {t i } N i=0 with a uniform time step ∆t. The time-discrete counterpart of (7) states that, given the displacement and the damage field (u i-1 , α i-1 ) at time step t i-1 , the solution at time step t i is obtained by solving the following bound-constrained minimization problem

inf{P t i (u, α) : u ∈ C(t i ), α ∈ D i }, (76)
where

D i = {α ∈ H 1 (Ω) : α(x) ≥ α i-1 a.e.}.
The unilateral constraint α(x) ≥ α i-1 is the time-discrete version of the irreversibility of damage.

The functional to be minimized at each time step, P t (u, α), is not convex in the pair (u, α). However, it is convex (and even quadratic for the model ( 75)) in each of the two variables u and α, once the other is fixed. Accordingly, the adopted numerical strategy consists in a series of alternate minimizations, that is minimization with respect to u at a fixed α and then minimization with respect to α at a fixed u, until convergence. Considering the displacement loading (3) with u t (x) = t i u e (x) imposed on a part of the boundary ∂ u Ω, the corresponding algorithm reads as follows:

• Initialization : Set (u (0) , α (0) ) := (u i-1 , α i-1 )

• Iteration p :

1. Compute, under the constraint u(x) = t i u e (x) on ∂ u Ω,

u (p) := arg min u P t i (u, α (p-1) ) (77) 2. Compute, under the constraint α i-1 ≤ α ≤ 1 on Ω, α (p) := arg min α P t i (u (p) , α) (78)
• End : Repeat until α (p) -α (p-1) ∞ ≤ δ 1 . Set (u i , α i ) := (u (p) , α (p) ).
The problem is discretized in space with standard finite elements. The finite element code is implemented in Matlab using compiled fortran and C functions for matrix assembly. Standard linear solvers are used to solve the displacement minimization (77). The damage minimization (78) with the irreversibility constraint is solved by the large-scale bound-constrained quadratic optimization solver included in the Matlab Optimization Toolbox. With a similar implementation we solved problem with up to 200 000 degrees of freedom. A much more efficient, massively parallelized, numerical implementation relying on High Performance Computing tools is being developed [START_REF] Bourdin | Numerical implementation of the variational formulation of quasi-static brittle fracture[END_REF].

Remark The numerical algorithm, and even the numerical code itself, is the same as the one used in the simulation of crack evolution in the variational approach to brittle fracture and it was firstly developed by [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and [START_REF] Bourdin | Numerical implementation of the variational formulation of quasi-static brittle fracture[END_REF]. Adopting the regularized version of brittle fracture proposed by [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functional via Gamma-Convergence[END_REF], the work of Bourdin et al. implements the damage model of Equation ( 48), which does not have a purely elastic phase. Moreover, it avoids the unilateral constrained minimization of equation ( 78) by imposing irreversibility only on fully developed cracks as bilateral constraints on the damage variable. For further details on the numerical applications of the variational approach to brittle fracture and its extensions, the reader may refer to Del Piero et al. ( 2007), [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF], and [START_REF] Lancioni | The variational approach to fracture mechanics. A practical application to the French Panthéon in Paris[END_REF].

Discretization and mesh independency

To obtain a discretized formulation of the minimization problem, we use standard linear triangular finite elements. In two-dimensional elasticity, each element has three nodal degrees of freedom: the two components of the displacement vector and α. Differing from most of numerical approaches to damage, the damage field is here treated as a nodal variable. The main goal of the gradient damage term is to render the numerical results independent of the mesh. However, to obtain this independency, the finite element meshes must be kept unstructured, uniform and fine enough.

Indeed, as pointed out by [START_REF] Negri | The anisotropy introduced by the mesh in the finite element approximation of the Mumford-Shah functional[END_REF], structured meshes cause a misestimation of the surface energy by introducing anisotropic effects which can favor some particular damage orientations. Similarly, non-uniform meshes can induce artificial inhomogeneities and favor the nucleation of cracks (in the sense of damage localizations) in refined zone. As far as the choice of the mesh size is concerned, an accurate estimation of the discretized dissipated energy requires the element size, say h, to be (much) smaller than the thickness of the localization zone (say, D 0 in our 1D analysis). Numerical experiments show that setting h ∼ E 0 ℓ/σ M still gives qualitatively reasonable results. However, the computed dissipated energy will be overestimated, see [START_REF] Bourdin | The variational approach to fracture[END_REF] and [START_REF] Bourdin | Numerical implementation of the variational formulation of quasi-static brittle fracture[END_REF] for details.

Even when these three conditions on the meshes (unstructured, uniform, fine enough) are fulfilled, it is not ensured that the numerical solution is independent of the mesh. Indeed, there is none general result of uniqueness and, more, there exist many situations where the number of solutions are infinite (see for example [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF]). In such a case, the numerical solution will depend on the mesh (and also on the algorithm), because small imperfections (due to physics and/or numerical errors) may greatly affect the final result. A key issue is to introduce in the numerical treatment the concept of stability of state in order to select the physically meaningful solutions 3 . This is a very fundamental but difficult task to which future works will be devoted. However, this type of mesh dependency is fundamentally different from the one obtained using non regularized models. For gradient damage models, it is closely related to the imperfection sensitivity of the physical system. On the other hand, for local damage models with stress-softening, mesh dependency is a pure consequence of the mathematical inconsistency of the modeling approach.

Numerical results for the 2D uniaxial traction test

We solve numerically the time-discrete evolution problem for the 2D version of the uniaxial traction test studied analytically in the previous sections. We consider a bar in plane-stress condition with an aspect ratio L/H = 10. Figure 6 illustrates the geometry and the boundary conditions. We consider the damage model ( 75), which has a finite elastic yield stress σ M reached for the end-displacement U e = L σ M /E 0 . The numerical simulations refer to dimensionless variables. Scaling the displacement with U e , the only relevant dimensionless parameter affecting the result of the simulation is the ratio:

λ = L ℓ . (79)
Following the stability diagram in figure 2(b), we distinguish short bars with L < λ c ℓ and long bars with L > λ c ℓ, being λ c = 4π/3 √ 3 ≃ 2.41. In the case of short bars, the homogeneous solution remains stable beyond the elastic limit U e . On the contrary, in the case of long bars, the homogeneous solution becomes unstable at the elastic limit U = U e , in a state with null damage field. We perform numerical simulations for long bars with L = 2λ c ℓ and short bar with L = λ c ℓ/2. In the case of a bar made of concrete with an internal length ℓ = 38 mm, as estimated in Section 3.3.5, this would correspond to bars of length L = 2λ c ℓ ≈ 180 mm and L = λ c ℓ/2 ≈ 46 mm. It would be almost impossible to perform experimentally the test on the short bar. This is the reason why [START_REF] Ramtani | Orthotropic behavior of concrete with directional aspects: modelling and experiments[END_REF] resorted to the PIED test (including aluminum reinforcements) to perform experiments with homogenous damage fields. The time step for loading increment is set to ∆T = 0.01. A uniform finite element mesh with a typical element size h = 0.01 is used. 4.4. Case of a long bar: L = 2 λ c ℓ For t < 1 all the material points are in an elastic regime. The alternate minimization algorithm gives a quasi-static evolution characterized by an elastic response (α = 0) with uniform deformation. For t > 1 we know that the elastic solution is unstable and the numerical algorithm gives a localized solution. Figure 7 shows the localized solution obtained for two types of boundary conditions on the damage field: (BC1) without Dirichlet boundary conditions on the its convergence toward unstable solutions. damage field, (BC2) with α = 0 on x = ±L/2. In the first case, the damage localization appears at the boundary. In the second case, the boundary conditions force the localization to appear inside the bar. Its precise location is energetically indifferent. Because of mesh imperfections and numerical errors, this location may vary from one simulation to another even if the initial data are identical (see the comments on mesh dependency in Section 4.2). As expected the response is essentially 1D, the field distribution being independent of y. The damage field profile versus x coincides to that obtained analytically in Section 3.3 (see Figure 5), where only the boundary conditions (BC1) were considered. In this case, only one half of the damage profile is sufficient to break the bar. It may be indifferently located at the left or the right end of the specimen. Its actual location in a numerical (and physical!) tests is ruled by imperfections and numerical errors. This type of solutions requires only one half of the energy of the full damage profile. Hence, it is a better candidate for energy minimization. Figure 8 reports the elastic and dissipated energies versus time, and the global force versus the imposed displacement at the end of the bar for the boundary conditions (BC1). Using Clayperon's theorem (20), the global force is calculated as being the time derivative of the total energy. The force-displacement diagram is linear until t = 1, then the force drops to zero. The dashed line of Figure 8(b) represents the homogeneous solution of equation ( 46). This homogeneous solution is stable only until the elastic limit t = 1. For t < 1, the elastic energy is quadratic in t and the dissipated energy is zero. For t > 1 the elastic energy vanishes and the dissipated energy is equal to G c H/2, in agreement with the analytical result of equation ( 69). For the boundary conditions (BC2), the diagrams are identical, except that the final dissipated energy is G c H. The numerical results show that the response of short bars is characterized by three phases: (i) elastic response with null damage for t < 1; (ii) progressive damage, with homogeneous damage distribution increasing with the loading for 1 < t < t * ; (iii) failure at a finite loading, say t * Figure 9 compares the time-evolution of the damage field distribution along the central 10 reports the energy diagram and the force-displacement relationship for L = λ c ℓ/2. For this case complete failure is obtained at t * = 1.65. The dissipated energy at the end of the test is close to G c H/2. The presence of the phase of progressive homogeneous damage is in agreement with the analytical result of Section 3.2.4, which predicts that homogeneous states are stable beyond the elastic limit when L < λ c ℓ. In particular the numerical force-displacement diagram in Figure 10(b) coincides with the one in figure 2(a) for t < t * . However, according to the analytical results of equation ( 47), the homogeneous solution is expected to be stable until a limit loading t s = λ c ℓ/L = 2. The numerical solution switches from the homogeneous state to a non-homogeneous one with vanishing residual stress at the smaller loading t * = 1.65. Applying the bifurcation analysis of [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF], it may be shown that t * is close to the loading for which a bifurcation from the homogeneous response occurs. (But the proof of that is beyond the scope of this paper.) Hence we may argue that, due to numerical errors and mesh inhomogeneities, the adopted numerical algorithm leaves the stable homogeneous solution branch as soon as a bifurcated non-homogeneous solution is available. A precise assessment of this point and of the bifurcation and stability loss phenomena in the numerical response is currently under investigation. 46), coinciding with the numerical response for t < t * = 1.65.

Concluding remarks

We have shown in this paper that gradient damage models are good candidates to approximate the fracture of materials. Inserting them in a variational setting allows to formulate the non local damage evolution problem in a very elegant and concise manner. Moreover, the concept of stability of states rises immediately from this variational approach. Thanks to the presence a characteristic material length, we can account both for size effects and for the stability of the homogeneous states for sufficiently small bodies. The determination of the size beyond which an homogeneous state is no more stable is one of the main step in the analysis of the properties of a given model. That gives an additional information to the stress-strain curve which is fundamental to identify the state functions entering in the model. Even if the major part of the results have been obtained here in a one-dimensional context, all these concepts can be easily extended to a full three-dimensional setting. These theoretical extensions will be the aim of future works. From a numerical viewpoint, computations of very complex geometry are already available and the interested reader can refer to [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], [START_REF] Bourdin | The variational approach to fracture[END_REF] or [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF].

Figure 1 :

 1 Figure 1: 1D traction problem with imposed end-displacement.

Figure 2 :

 2 Figure 2: Properties of the homogeneous solutions for example 1.

Figure 3 :

 3 Figure 3: Properties of the homogeneous solutions for example 2.

Figure 4 :

 4 Figure 4: Stability diagram for Example 3 with p = 2, 4, 8.

Figure 5 :

 5 Figure 5: The damage profile α(x) centered at x0 for σ = 0 when considering the damage law (43), see equation (72).

Figure 6 :

 6 Figure 6: 2D uniaxial traction test with imposed displacements.

4. 5 .

 5 Case of a short bar: L = λ c ℓ/2

Figure 7 :

 7 Figure7: Damage field distribution in a long bar (L = 2 λcℓ) for t > 1. We report the solutions obtained for the damage field free at the boundary (BC1) and the damage field imposed to be zero on x = ±L/2 (BC2).

Figure 8 :

 8 Figure 8: Numerical results for the 2D traction of a long bar (L = 2 λcℓ) with damage free at the boundary (BC1). The dashed line in the force-displacement diagram reports the 1D homogeneous solution in equation (46), coinciding with the numerical response for t < 1.

  Figure

  Figure 10: Numerical results for the 2D traction of a short bar (L = 2 λcℓ) with damage free at the boundary (BC1). The dashed line in force-displacement diagram reports the 1D homogeneous solution in equation (46), coinciding with the numerical response for t < t * = 1.65.

The continuity of α ′ at x0 ± D is obtained as a first order optimality condition on P(u, α).

[START_REF] Ramtani | Orthotropic behavior of concrete with directional aspects: modelling and experiments[END_REF] designed a specific experimental apparatus to perform homogenous uniaxial test preventing damage localization, the so-called PIED (Pour Identifier l'Endommagement Diffus) test. It allows for the experimental measure of a homogeneous stress-strain field of concrete by sticking aluminium plates on the lateral faces of the specimen. The aluminium reinforcements stabilize the diffuse damaged response (microcracks) and prevents the appearance of transverse macrocracks.

The present alternate minimization algorithm does not explicitly test stability. However, we rarely observed

Acknowledgements

The present work has been developed in the framework of the French research network MoMaS (Modélisations Mathématiques et Simulations numériques liées aux problèmes de gestion des déchets nucléaires). Its financial support is gratefully acknowledged.

Appendix: Proof of the Clayperon Theorem (20)

As the couple (u t , α t ) is solution of the evolution problem, it verifies in particular the balance energy (12c) at each time t

which after integrating by parts and using the boundary conditions

reads

Inserting into the time derivative of the total energy, we obtain

Using the equilibrium equations ( 10) and the displacement boundary conditions, the right hand side of (83) becomes

Integrating in time the above relation and using the initial condition α 0 (x) = 0 lead to equation (20). The equality ( 22) is a simple consequence of (84).