Arthur Geromel Fischer

Jean-Jacques Marigo

Gradient damage models coupled with plasticity and their application to dynamic fragmentation

Keywords: Plasticity, Damage, Variational approach, Dynamical fracture, Fragmentation

Introduction

The initiation and propagation of is still an unresolved question in fracture mechanics. Several models have been studied in different contexts [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF]1,[START_REF] Donzé | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF][START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF], in quasi static and dynamics [START_REF] Ravi-Chandar | Dynamic fracture of nominally brittle materials[END_REF][START_REF] Christopher | Models for dynamic fracture based on griffith's criterion[END_REF], and accounting for different phenomena. The objective of this chapter is to explain the development of the so-called "gradient damage models" [START_REF] Amor | Gradient damage models and their use to approximate brittle fracture[END_REF] and its extension to ductile materials and dynamic loading. The main idea of these models is that we represent the by a and and we do not need to know a-priori its path. However, some attention must be payed to the mesh in order to avoid creating a preferential direction for [START_REF] Negri | The anisotropy introduced by the mesh in the finite element approximation of the mumford-shah functional[END_REF].

Local models have proven to be not suitable to correctly predict [START_REF] Geers | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF][START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. Softening local damage models allow in infinitely thin bands and, consequently, cracks with zero energy dissipation [START_REF] Marigo | Bifurcation and stability issues in gradient theories with softening[END_REF]. In finite elements simulations, this implies that the mesh size determines the size of the localization zones and the results will necessarily depend on the mesh used.

In this context, the localization is controlled by adding a nonlocal term as an integral [START_REF] Pijaudier | Nonlocal damage theory[END_REF][START_REF] Zdenek | Nonlocal continuum damage, localization instability and convergence[END_REF][START_REF] De Borst | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF] or a gradient [START_REF] Comi | A non-local model with tension and compression damage mechanisms[END_REF][START_REF] Lorentz | Gradient constitutive relations: numerical aspects and application to gradient damage[END_REF][START_REF] Lorentz | Convergence of a gradient damage model toward a cohesive zone model[END_REF]. The family of gradient damage models contain the gradient of damage weighed by a parameter called the "" [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF] in order to avoid a localization in a band of null thickness.

These models have been originally proposed for quasi-static , but have also been extended to [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF] and dynamic loading [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Michael | A phase-field description of dynamic brittle fracture[END_REF]21]. The main objective of this chapter is to explain the necessary changes of the original model, in order to account for and inertial effects.

We first present the construction of gradient damage models for brittle softening materials based on the principle of . We discuss the main hypothesis and the need for regularization. We then briefly talk about Von-Mises criterion [START_REF] Marigo | Plasticité et Rupture[END_REF] and how to take it into account [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF]. We conclude the model by adding . Once the model is complete, we will briefly talk about the numerical implementation and show a few examples and results using the FEniCS library [START_REF] Mardal | Automated Solution of Differential Equations by the Finite Element Method -The FeniCS Book[END_REF] and an industrial code.

Theoretical Aspects

Gradient Damage Models

We present here a simplified construction of for brittle where there are no other dissipation effects. We are going to consider the case of theory and an material. For a more detailed construction of these models, see [START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Pham | Construction et analyse de modèles d'endommagement à gradient[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]. For the proofs of Gamma-convergence, consult [START_REF] Braides | Γ-convergence for beginners[END_REF][START_REF] Dal-Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF].

We denote the by σ, the by ε, the by u and the rigidity tensor by E. The contracted product of two tensors a and b will be denoted by a:b. When working in a 1-D scenario, we are going to call the Young's Modulus simply E.

We recall that ε = 1 2 (∇u + ∇ T u) and σ = E:ε.

Construction of a Damage Model (non-regularized)

The objective of this section is to describe a family of that can be applied to different types of materials. We will discuss the qualitative properties of these models.

We will suppose that 1. Damage can be represented by a scalar α ∈ [0, 1]. When α=0 the material is healthy and when α=1 the material is completely broken. Other choices for α are possible, for instance α ∈ [0, ∞) or α being a tensor, but we decided to keep the model here as simple as possible.

2. The rigidity tensor E(α) is a function of α, the material becomes less rigid when α increases and E(α=1) = 0 (no rigidity left when the material is broken). It is important to notice that, for a fixed damage value, the stress-strain relation is supposed to be linear.

3. Damage is , that is, it can only grow in time.

We now need to specify the conditions for . For that, we are going to use an idea similar to [START_REF] Griffith | The phenomena of rupture and flows in solids[END_REF], based on the notion of elastic , in its [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF].

The can be written as

ψ(ε, α) = 1 2 ε:E(α):ε. (1)
For a fixed deformation, a small increase δα>0 of damage causes a loss of -∂ψ ∂α (ε, α)δα>0 in the elastic energy. We compare the variation of elastic energy to a threshold k(α). As in Griffith's model, the rate of energy restitution is always smaller or equal to a threshold value and the crack only propagates when we have an equality. For this family of damage models, the propagation criterion can be written as

- 1 2 ε:E (α):ε ≤ k(α), α = 0 if -1 2 ε:E (α):ε < k(α) α ≥ 0 if -1 2 ε:E (α):ε = k(α) (2)
where k(α) is a function of α representing the necessary energy restitution necessary for damage to evolve. So far, we made only a few hypothesis concerning E(α) and k(α). One of the strengths of such models is that the specific forms of E(α) can be chosen such as to describe a specific material's behavior.

Example 2.1. As a first example, we are going to consider a 1-D bar under traction, where damage increases uniformly in space. We consider the functions

E(α) = E 0 (1 -α) 2 and k(α) = k 0 , (3)
where k 0 is a constant and E 0 is the Young's modulus when the material has not yet suffered any damage. It is easy to see that the damage criterion can be written as

(1 -α)E 0 ε 2 ≤ k 0 , α = 0 if (1 -α)E 0 ε 2 < k 0 α ≥ 0 if (1 -α)E 0 ε 2 = k 0
If the initial bar is not damaged, the bar will not be damaged while E 0 ε 2 <k 0 . We can define a ε c and critical stress σ c by

ε c = k 0 E 0 and σ c = E 0 ε c = k 0 E 0 .
When ε>ε c we can use the criterion once again to find that (1 -α)E 0 ε 2 = k 0 and thus

α = 1 - k 0 E 0 ε 2 .
For this model, we can see that while ε≤ε c , stress increases linearly with the strain and the bar suffers no damage. When ε>ε c , the bar damages and this causes a stress-strain relation that is no longer linear.

Figure 1 shows the normalized σ = σ/σ c and damage α as a function of the normalized ε = ε/ε c .

Let w(α) be a function such that w (α) = k(α). We define the energy density by

W (ε, α) = ψ(ε, α) + w(α). (4)
We can write the stress as

σ = ∂W ∂ε (ε, α) (5)
and the damage evolution criterion as

∂W ∂α (ε, α)• α = 0, (6)
where each of the two factors is non-negative. Now let β ≥ 0 be a small increase of damage in time. We have that

∂W ∂α (ε, α)•( β -α) ≥ 0. ( 7
)
Consider a structure whose initial configuration is given by Ω ⊂ R n (n=1, 2, 3). Suppose we have a volume force f acting on the whole structure, an imposed displacement u 0 on ∂ u ⊂ ∂Ω and a normal stress T on ∂ T ⊂ ∂Ω. We also suppose that

∂ u ∂ T = ∅ and ∂ u ∂ T = ∂Ω.
The static equilibrium can be written as

     divσ + f = 0 on Ω u = u D on ∂ u σ•n = T on ∂ T . (8)
We fix a test function w such that w = 0 on ∂ u . Then

Ω divσ•w + f •w dΩ = 0 (9)
and Green's formula shows that

∂u (σ•n)•wdS 0 + ∂T T •wdS - Ω σ : ε(w)dΩ + Ω f • wdΩ = 0. ( 10
)
We define

C = {u : u=u D on ∂ u } C 0 = {w : w=0 on ∂ u } . ( 11
)
The static equilibrium problem consists of finding u ∈ C such that

Ω ∂W ∂ε (ε(u), α):ε(w)dΩ = Ω f • wdΩ + ∂T T •wdS, ∀w ∈ C 0 . ( 12
)
If we consider the evolution problem where the time is denoted by t, by integrating [START_REF] Marigo | Bifurcation and stability issues in gradient theories with softening[END_REF], we obtain the following problem: find α≥0 such that

Ω ∂W ∂α (ε, α)•( β -α)dΩ ≥ 0, ∀ β ≥ 0. ( 13
)
We define the total energy of the system by

E(u, α) = Ω W (ε(u), α)dΩ - Ω f • udΩ - ∂T T •udS. (14)
It is easy to see that the evolution problem, given by equations ( 12) and ( 13), is equivalent to finding u ∈ C and α≥0 such that for all v ∈ C and β≥0 we have

∂E ∂u (u, α)(v -u) ≥ 0 and ∂E ∂α (u, α)( β -α) ≥ 0. ( 15
)

Regularized Model

It is now a well-known fact that local softening damage models are not viable ( [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], [START_REF] Amor | Gradient damage models and their use to approximate brittle fracture[END_REF]) as they allow damage localization in infinitely thin bands.

Example 2.2. To illustrated the problem of damage localization, we are going to consider a 1-D bar of length L and a material such that E(α)=E 0 (1 -α) 2 and w(α)=w 1 α.

When in equilibrium, we know that σ(x) = σ (constant). We will show that for any 0 < θ < 1 fixed, we can construct a solution to the damage problem such that there is no damage in the interval (0, θL) and uniform damage in (θL, L).

In fact, for x ∈ (0, θL), we have ε(x) = σ/E 0 . For x ∈ (θL, L), the damage criterion can be written as

w 1 = - 1 2 E (α) σ 2 E(α) 2 = E 0 (1 -α) σ 2 E 2 0 (1 -α) 4 = σ 2 E 0 1 (1 -α) 3 . (16)
Therefore, damage in this interval is given by

α * = 1 -3 σ 2 w 1 E 0 . ( 17
)
The dissipated energy can be calculated

D = L 0 w(α)dx = L θL w 1 α * dx = w 1 α * (1 -θ)L (18)
This shows that we have a solution of the damage problem for any θ. We can see that damage can be localized in an infinitely thin band and if we take θ → 1, the dissipated energy D tends to zero.

In a finite elements code, the size of the damage band will be determined by the mesh size. This means that refining the mesh will produce different results and dissipated energies that can tend to zero.

To solve this problem, a regularizing term is used. The main idea is to add a characteristic length in order to penalize sharp damage profiles and solve the problem of localization in thin bands and fracture without energy dissipation.

One simple way of doing this is by adding a term that depends of the norm of the gradient of damage.

We will use an energy density of the form

W (ε, α, ∇α) = ψ(α, ε) + w(α) + 1 2 w 1 2 ∇α•∇α, (19)
where is the and w 1 >0 is a normalization constant.

In the previous section, when describing the model, we first proposed an evolution law based on the energy restitution rate. We then expressed the static equilibrium and damage evolution by a principle of . For this new energy density, we are going to use directly the principle of minimum energy to obtain an evolution law, instead of manually proposing it.

We have

σ = E(α):ε = ∂W ∂ε (ε, α). ( 20
)
We define the dissipated energy by

D(α) = Ω w(α) + 1 2 w 1 2 ∇α•∇αdΩ (21)
and redefine the total energy

E(u, α) = Ω W (ε(u), α, ∇α)dΩ - Ω f • udΩ - ∂T T •udS. (22)
The evolution problem consists of finding u ∈ C and α≥0 such that

DE(u, α)(v -u, β -α) ≥ 0, ∀v ∈ C, ∀ β ≥ 0. ( 23
)
Example 2.3. We will show an example where we have a 1-D bar that breaks only in a small region. For more details, the reader is referred to [START_REF] Marigo | From the onset of damage to rupture: construction of responses with damage localization for a general class of gradient damage models[END_REF]. We are will assume that fracture occurs in the interval [x 0 -∆, x 0 + ∆], where x 0 is the center of the damage profile (assumed to be symmetric) and the value of ∆ is unknown.

Since the whole region is damaged, α satisfies the damage criterion in the interval:

1 2 E (α)ε 2 + w (α) -w 1 2 α = 0. ( 24
)
If we write S(α) = 1/E(α), then (omitting α) S = -1/E 2 E . Thus E = -E 2 S and

- 1 2 S (α)σ 2 + w (α) -w 1 2 α = 0. ( 25
)
We multiply this expression by α and integrate to obtain

- 1 2 S(α)σ 2 + w(α) -w 1 2 α 2 2 = constant = - σ 2 2E(0) . ( 26
)
We define the function H by

H(α) := σ 2 1 w 1 E(0) - S(α) w 1 + 2w(α) w 1 . ( 27
)
Then 2 α 2 = H(α). ( 28
)
The maximum value α max of α is such that H(α max ) = 0. For a given α max , we can find the value of the stress σ on the bar:

σ = 2w(α max ) S(α max ) -1/E(0) . ( 29
)
We can then find the damage profile using the relation

x(α) = x 0 αmax α H(β) dβ. ( 30
)
The different damage profiles are shown in Figure 2.3. Consider now that E(α)=(1 -α) 2 and w(α)=w 1 α. When the material in completely broken, we have σ=0 and H(α) = 2α. We then solve 2 α 2 = 2α, that is

dα √ α = √ 2dx, (31)
which, considering that α(

x 0 ) = 1 gives us √ 2 (1 - √ α) = x -x 0 (32)
We thus find for x ∈ (x 0 -∆, x 0 + ∆)

α(x) = 1 - |x -x 0 | √ 2 2 and ∆ = √ 2 . ( 33
)
The dissipated energy can now be calculated:

D = x 0 +∆ x 0 -∆ w 1 (α(x)) + 1 2 w 1 2 α 2 dx = 2 √ 2 0 2w 1 1 - y √ 2 2 dy. (34)
Thus

D = 4 √ 2 3 w 1 . (35)

Numerical Implementation of Damage

In the previous subsection, we described the model for brittle damage evolution. For a given set of boundary conditions, we want to minimize the total energy E(u, α) with respect to u and α, with the constraint α≥0 (irreversibility condition).

We will consider a problem discrete in time, that is, we divide the time in instants t 0 , t 1 , t 2 , ..., with t i < t i+1 , and study the evolution of the system. At an instant t i , we have a displacement u i and damage u i .

In the , u i and α i must minimize the E of the system. We recall that α i (x) ∈ [0, 1], for all x ∈ Ω, and if the damage at an iteration i is α i , then the condition is equivalent to α i+1 (x) ≥ α i (x), for all x ∈ Ω.

We want to E for u and α at the same time. The functional E may not be convex for the pair (u, α). It is, however, convex for each variable. That's why we use an [START_REF] Marigo | Numerical experiments in revisited brittle fracture[END_REF] procedure, that is, we alternately E in u for a fixed α and then, for a fixed u, minimize E in α, until we obtain convergence for the pair (u, α).

This method only guarantees that the solution is a stationary point or, in some cases, a local minimum of E. It does not ensure the of the solution.

To calculate (u i , α i ), we have to know the values of u i-1 and α i-1 . At each time iteration i, we have to:

(1) Update boundary conditions.

(2) Set (u (i,0) , α (i,0) ) := (u i-1 , α i-1 ).

(3) Iteration j ≥ 1:

(3.1) Solve u (i,j) := arg min u E(u, α (i,j-1) ); (3.2) Solve α (i,j) := arg min α i-1 ≤α≤1
E(u (i,j) , α);

(3.3) Stop when α (i,j) -α (i,j-1) is sufficiently small.

(4) We then set (u i , α i ) := (u (i,j) , α (i,j) ).

Note: at each minimization, u and α must respect boundary conditions. One important feature of this is that the energy decreases at each iteration:

E(u (i,j) , α (i,j) ) ≤ E(u (i,j-1) , α (i,j-1) ). (36)
There is no guarantee that the will converge in the general case and the choice of convergence criterion may influence the results. There are different aspects we can consider when choosing the convergence criterion:

• the convergence criterion depends on the variables or the total energy;

• relative tolerance or absolute tolerance;

• norm used: we consider only the difference between elements, the sum all over the structure, whether the gradient of the variables is important;

Damage Coupled with Plasticity

The family of models we have developed so far cannot account residual strains. In this section, we want to extend the damage models described in section 2.1.2 to materials.

The main idea is to define the total energy of the body in terms of the displacement field, the damage, the plastic strain and the cumulated plastic and study the of the system through the minimization with respect to each variable.

Section 2.1.2 was a review of the used to describe damage evolution for a quasi-static loading. We start this section by doing the same for . We then present and discuss an energy expression that can couple these two phenomena and show one example of material behavior.

Perfect Plasticity Model

We will denote the plastic strain by ε p . The stress is now given by

σ = E:(ε -ε p ). ( 37
)
In the general case, σ is admissible if it satisfies f (σ)≤0, where the function f depends on the criterion used. The evolution law is given by the relation

εp •f (σ) = 0. ( 38
)
In the case of yield criterion, we have in 1-D

f (σ) = |σ| -σ Y (39)
and in 3-D

f (σ) = 3 2 s:s -σ Y , (40)
where s := σ -Trσ 3 I is the and σ Y is the and is a material constant. We define the cumulated plastic strain from zero to the instant t as

pt = p0 + t 0 εp dτ (41)
and the energy density for a elasto-plastic material as

W 1D (ε, ε p , p) = 1 2 E(ε -ε p ) 2 + σ Y p (42)
in 1-D and as

W 3D (ε, ε p , p) = 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y p (43) in 3-D.
We consider a t 0 , t 1 , ... and we can calculate the cumulated plasticity by

pi = pi-1 + (ε p ) i -(ε p ) i-1 (44) 2.2.2. Plasticity Evolution in 1-D We define the function f (ε, p) = 1 2 Ep 2 -Eεp + σ Y |p -(ε p ) i-1 |. ( 45
)
It is clear the the minimization of f in p is equivalent to the minimization of

W 1D (ε(u), ε p ) in ε p .
The function f is strictly convex in p and is differentiable everywhere except in p = (ε p ) i-1 . As a consequence, f has one unique minimum.

We use two auxiliary results:

Proposition 2.4. For a given ε, set σ * = E(ε -(ε p ) i-1
). The value p that minimizes f (ε, p) can be characterized by:

(

1) If |σ * | ≤ σ Y , then the minimum is attained in (ε p ) i-1 . (2) If |σ * | > σ Y , then the minimum is attained at a point such that ∂f ∂p (ε, p) = 0.
Proof. We write p = (ε p ) i-1 + e. Then

f (ε, p) = f (ε, (ε p ) i-1 ) + 1 2 Ee 2 -σ * e + σ Y |e|. ( 46
) (1) If |σ * | ≤ σ Y , then σ * e ≤ σ Y |e| and f (ε, p) ≥ f (ε, (ε p ) i-1 ) + 1 2 Ee 2 .
Hence, the minimum is attained when e = 0, that is, when p = (ε

p ) i-1 . (2) If |σ * | > σ Y , we put e = hσ * /|σ * |, with h > 0. Then f (p) = f (ε, (ε p ) i-1 ) + 1 2 Eh 2 -σ * h + σ Y h. ( 47
) If h is small enough, then f (ε, p) < f (ε, (ε p ) i-1
). Since f is regular everywhere except e = 0, we must have ∂f ∂p (ε, p) = 0.

Proposition 2.5. In the evolution problem, we set

σ * = E(ε i -(ε p ) i-1 ). The minimization of W in ε p is equivalent to: (1) If |σ * | ≤ σ Y : (ε p ) i = (ε p ) i-1 . ( 48
) (2) If |σ * | > σ Y : (ε p ) i = (ε p ) i-1 + 1 - σ Y |σ * | ε i -(ε p ) i-1 (49)
and

E(ε i -(ε p ) i ) = σ Y . (50)
Proof. We have already proved (1) in proposition (2.4).

To prove (2), again by proposition (2.4), we have to find p such that ∂f ∂p (ε i , p) = 0. We notice that for e = 0 and |δe| < |e|, we have

|e + δe| = |e| + δe e |e| . ( 51
) Then ∂f ∂p (ε i , p) = Ee -σ * + σ Y e |e| = 0. ( 52
)
Hence,

E(p -(ε p ) i-1 ) -E(ε i -(ε p ) i-1 ) + σ Y e |e| = 0. ( 53
)
Rearranging the terms,

E(ε (i,j) -p) = σ Y e |e| . ( 54
e |e| = σ * -σ |σ * -σ| = σ * |σ * | . ( 56
)
Finally, by ( 53),

e = 1 E σ * -σ Y e |e| = 1 E σ * -σ Y σ * |σ * | (57) and (ε p ) i := p = (ε p ) i-1 + e = (ε p ) i-1 + 1 -σ Y σ Y |σ * | (ε i -(ε p ) i-1 ). ( 58
)
This is an elastic prediction -plastic correction procedure: we calculate the current strain and stress based on the previous time instant assuming that the material is (elastic prediction). If the stress is inside the , that is, |σ * | < σ Y , we keep it and the does not change. On the other hand, if the stress is not in the elastic domain, we update the ().

This behavior is shown in Figure 3. We can see the normalized strain ε:=εE 0 /σ Y , normalized stress and normalized plastic strain.

Plasticity Evolution in 3-D

In this section, we are going to write the same results as in the previous section, but now to a problem in 3-D.

We consider the usual assumption of plastic incompressibility:

tr(ε p ) = 0 (59) and 3 2 s:s ≤ σ Y , (60)
where the deviatoric stress tensor s is given by

s = σ - tr(σ) 3 I. ( 61
)
The total energy density is written as

W (ε, ε p ) = 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y ε p -ε p , (62)
where e = √ e:e.

We now define e as the deviatoric part of ε and since tr(ε p ) = 0, the minimization of W in ε p is equivalent to min

p : tr(p)=0 f (p), for every point in Ω ( 64
)
where

f (ε, p) := µp:p -2µe:p + 2 3 σ Y p -(ε p ) i-1 . ( 65
)
(The Lamé's coefficients are denoted by λ and µ.)

We set σ * = E:(ε i -(ε p ) i-1 ) and its deviatoric part is given by s * = 2µ(e -(ε p ) i-1 ).

The following propositions are the 3-D equivalents of the auxiliary results in section 2.2.3:

Proposition 2.6. The value p that minimizes f (ε, p) can be characterized by:

(1) if s * ≤ 2 3 σ Y , then the minimum is attained in (ε p ) i-1 ;

(2) if s * > 2 3 σ Y , then the minimum is attained at a point such that ∂f ∂p (ε, p) = 0.

Proof. We write p = (ε p ) i-1 + δ. Then

f (ε, p) := µp:p -2µe:p + 2 3 σ Y p -(ε p ) i-1 + f ((ε p ) i-1 ) -µ(ε p ) i-1 :(ε p ) i-1 + 2µe:(ε p ) i-1 = f ((ε p ) i-1 ) + µδ:δ + 2µp:(ε p ) i-1 -2µ(ε p ) i-1 :(ε p ) i-1 -2µe:δ + 2 3 σ Y δ = f ((ε p ) i-1 ) + µδ:δ + 2µδ:(ε p ) i-1 -2µe:δ + 2 3 σ Y δ = f ((ε p ) i-1 ) + µδ:δ + 2 3 σ Y δ -s * :δ. ( 66
) (1) If s * ≤ 2 3 σ Y , then f (ε, p) ≥ f (ε, (ε p ) i-1 ) + µδ:δ. Hence, the minimum is attained when δ = 0. (2) If s * > 2 3 σ Y , we put δ = hs * / s * , with h > 0. Then f (p) = f ((ε p ) i-1 ) + 2 3 σ Y h -s * h + µh 2 . ( 67
)
If h is small enough, then f (p) < f ((ε p ) i-1 ). Since f is regular everywhere except δ = 0, we must have ∂f ∂p f (p) = 0.

Proposition 2.7. The minimization of W in ε p is equivalent to:

(

1) If s * ≤ 2 3 σ Y : (ε p ) i = (ε p ) i-1 . ( 68
) (2) If s * > 2 3 σ Y : (ε p ) i = (ε p ) i-1 + 1 - 2 3 σ Y |s * | e i -(ε p ) i-1 . ( 69
)
Proof. The proof of (1) follows directly from the last proposition.

To prove (2), we have to find p such that ∂f ∂p (ε, p) = 0.

We derive f and apply it to a tensor δ:

∂f ∂p (ε, p):δ = 2µ(p -e i ):δ + 2 3 σ Y p -(ε p ) i-1 p -(ε p ) i-1 :δ = 0. (70) If σ i = E:(ε i -p) (71)
and s i is its part, we must have

s i = 2 3 σ Y p -(ε p ) i-1 p -(ε p ) i-1 . ( 72
)
It is clear that

s i = 2 3 σ Y . ( 73
)
We note that

s i = s * + 2µ((ε p ) i-1 -p) (74)
and, by equation ( 72), s i and s * have the same direction. Since we know s * , we obtain

s i = 2 3 σ Y s * s * . (75)
Finally, applying this to (72) and (74),

p -(ε p ) i-1 = s i p -p i-1 2 3 σ Y = s i s i -s * 2µ 2 3 σ Y = s * s i -s * 2µ s * = 1 - 2 3 σ Y s * s * 2µ = 1 - 2 3 σ Y s * (e i -(ε p ) i-1 ). ( 76
)
We conclude by taking (ε p ) i = p.

Damage-Plasticity Coupling

In this section, in order to construct a family of models that account for and damage, instead of proposing the for each variable, we work directly with a suitable form of and, by minimizing this energy, we deduce the . For simplicity, we remove volume forces from our calculations In section 2.1.2, is was postulated a total energy for brittle damage:

E brittle (u, α) = Ω ψ(α, ε(u)) + w(α) + 1 2 w 1 2 ∇α•∇α dΩ. ( 77
)
We recall that the evolution of the system for loading can be obtained minimizing this energy with respect to u and α. A perturbation in the direction u gives us the and a perturbation in α gives us the .

In section 2.2.1, we showed that the evolution of the plasticity minimizes the energy

E 1D plast (ε, ε p ) = Ω 1 2 E(ε -ε p ) 2 + σ Y p dΩ (78)
in 1-D, and

E 3D plast (ε, ε p ) = Ω 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y p dΩ (79)
in 3-D. By examining perturbations in ε and ε p , obtain the static equilibrium and the , respectively. As we can see, the problems for damage and plasticity are similar in the sense that the quasistatic evolution in both cases in found after minimizing the total energy. For the , we are going to use an energy form that is, in a way, a combination of the damage energy and the plastic energy. For that, we are going to assume that the now depends on the damage, that is, σ Y =σ Y (α).

We define the the following 1-D and 3-D energies for the damage-plasticity (DP) coupling:

E 1D DP (ε, ε p , p, α) = Ω 1 2 E(α)(ε -ε p ) 2 + σ Y (α)p + w(α) + 1 2 w 1 2 α 2 dΩ ( 80
)
and

E 3D DP (ε, ε p , p, α) = Ω 1 2 (ε -ε p ):E(α):(ε -ε p ) + 2 3 σ Y (α)p + w(α) + 1 2 w 1 2 |∇α| 2 dΩ. (81)
To obtain the evolution criteria, we the with respect to 3 variables (u, ε p and α):

• The minimization of the displacement gives us the static-equilibrium:

divσ = 0 , where σ = E(α):(ε -ε p ). (82)
• The minimization of the gives us

3 2 s:s ≤ σ Y (α) and εp • 3 2 s:s -σ Y (α) = 0. ( 83
)
• The minimization of α gives us the new damage criterion (after taking the derivative with respect to α and integrating by parts). In 1-D:

1 2 E(α )(ε -ε p ) 2 + σ Y (α)p + w(α ) -w 1 2 α ≥ 0 (84) In 3-D: 1 2 (ε -ε p ):E (α):(ε -ε p ) + 2 3 σ Y (α)p + w (α) -w 1 2 ∆α ≥ 0 ( 85
)
We also have α = 0 when we have a strict inequality.

Example 2.8. Consider a bar given by Ω=[0, L] under traction, where the displacement at the extremities are controlled. We want to calculate the evolution of damage and plastic strain for the homogeneous case. We consider σ 0 Y < √ w 1 E 0 . We take the functions

E(α) = E 0 (1 -α) 2 , w(α) = w 1 α , σ Y (α) = σ 0 Y (1 -α) 2 . ( 86
)
Since we are assuming uniformity in space, we have to calculate the scalars σ, ε, ε p and α.

We have 3 different stages: • elastic phase: it is easy to see that while ε < σ 0 Y /E 0 , then σ < σ Y (α)=σ 0 Y and there is no change in the plastic strain. Since there is no plastic strain, the damage criterion is the same for brittle materials and we see that the bar does not suffer any damage.

• plastic phase: if ε > σ 0 Y /E 0 , then plastic strain evolves. In a pure traction test, the plastic strain and the cumulated strain are the same and we must have

E 0 (ε-ε p )=σ 0 Y . Thus ε p =ε-σ 0 Y /E 0 . The damage criterion becomes -(1 -α) (σ 0 Y ) 2 E 0 -2(1 -α)σ 0 Y p + w 1 ≥ 0. ( 87
)
It is easy to see that for α=0, we have a strict inequality while

ε p < w 1 2σ 0 Y - σ 0 Y 2E 0 .
• damage-plastic phase: the plasticity continues to evolve and the criterion gives us ε p =ε-σ 0 Y /E 0 . The damage criterion is now

-(1 -α) (σ 0 Y ) 2 E 0 -2(1 -α)σ 0 Y p + w 1 = 0. ( 88
)
We can thus find

α = σ 0 Y E 0 + 2p -w 1 σ 0 Y σ 0 Y E 0 + 2p . ( 89
)
Figure 4 shows these three phases. We see the normalized (in function of damage threshold) σ = σ/σ c and strain ε = σ/ε c . We can clearly identify the three phases in the stress curve.

Dynamic Gradient Damage

To formulate the evolution of the dynamic system, we are going to use the principle of least action, as in [START_REF] Li | Gradient Damage Modeling of Dynamic Brittle Fracture[END_REF].

Suppose we have a mechanic system Ω whose displacement is u and stress is σ(u). At each instant t ∈ [t 1 , t 2 ] we impose a displacement u D (t) on ∂ u ⊂ ∂Ω and a normal stress T (t) on ∂ T ⊂ ∂Ω. We also suppose that ∂ u ∂ T = ∅ and ∂ u ∂ T = ∂Ω. We have the following equations:

     ρü = divσ + f on Ω u = u D (t) on ∂ u σ•n = T (t) on ∂ T . (90)
We fix a test function w such that w(x, t)=0 on ∂ u for all t ∈ [t 1 , t 2 ] and w(t=t 1 ) = w(t=t 2 ) = 0 on Ω. Then

Ω ρü•wdΩ = Ω divσ(u)•w + f •w dΩ (91)
and Green's formula shows that

Ω ρü • wdΩ = ∂u (σ•n)•wdA 0 + ∂T T •wdA - Ω σ(u):ε(w)dΩ + Ω f •wdΩ. ( 92
)
We integrate this equation between instants t 1 and t 2 , and after an integration by parts, we obtain

Ω ρ u•wdΩ t 2 t 1 - t 2 t 1 Ω ρ u• ẇdΩ dt = t 2 t 1 ∂T T •wdA dt - t 2 t 1 Ω σ(u):ε(w)dΩ dt + t 2 t 1 Ω f •wdΩ dt. (93)
We define the kinetic energy of the system

K( u) = Ω 1 2 ρ u 2 dΩ (94)
and the potential energy

P(u) = 1 2 Ω σ(u):ε(u)dΩ - Ω f •udΩ - ∂ T T •udA. ( 95
)
Applying the boundary conditions of w in t 1 and t 2 , we have

t 2 t 1 ∂P ∂u w - ∂K ∂ u ẇ dt = 0, ∀w. (96)
We have thus shown that the problem (90) implies equation (96). It is easy to see that equation (96) can be obtained by searching for stationary points of an action functional defined by

S(u, u) = t 2 t 1 P(u(t)) -K( u(t))dt. ( 97
)
This motivates us to construct a dynamic gradient model by defining a suitable form of the action functional. Instead of using a purely σ:ε, we are going to use the energies defined by equations ( 80) and (81) with the terms containing the plastic strain and energy dissipated by the damage process. We recall that they were written as

E 1D DP (ε, ε p , p, α) = Ω 1 2 E(α)(ε -ε p ) 2 + σ Y (α)p + w(α) + 1 2 w 1 2 α 2 dΩ and E 3D DP (ε, ε p , p, α) = Ω 1 2 (ε -ε p ):E(α):(ε -ε p ) + 2 3 σ Y (α)p + w(α) + 1 2 w 1 2 |∇α| 2 dΩ.
We take the external loads into account and define a as

P DP (u, ε p , α) = E DP (ε(u), ε p , p, α) - Ω f • udΩ - ∂ T T • udA. ( 98
)
We define the new by

L DP (u, u, ε p , p, α, t) = P DP (u(t), ε p (t), p, α(t)) -K( u(t)) (99)
and the by

S DP (u, u, ε p , p, α) = t 2 t 1 L DP (u, u, ε p , p, α, t)dt. (100)
We define the admissible displacement space C and admissible damage space D by

C = {u : u(t)=u 0 (t) on ∂ u } D = {α ∈ [0, 1] : α ≥ 0 on Ω} (101)
In order to preserve the of damage and plasticity, instead of searching for stationary points, we will now only consider the unilateral minimal condition of the action, that is, we search an displacement u∈C, damage α∈D and ε p such that S DP (u, u, ε p , p, α) ≤ S DP (w, ẇ, p, p -ε p + p, β)

for any w∈C, β∈D and p.

In particular, if we take β=α and p=ε p , we must have

∂S DS ∂u (w -u) + ∂S DS ∂ u ( ẇ -u) = 0 (103)
and, by following the previous calculations in reverse order, we find the problem given by (90). We now set w=u and p=ε p to study the . If at an instant t the damage is α t then we define the admissible damage D t taking α t and the irreversibility condition into account:

D t = {β : β ≥ 0 and β ≥ α t on Ω}. ( 104
)
For every

β∈D t ∂S DS ∂α (β -α) ≥ 0. ( 105
)
From this, it is easy to see that we obtain the same damage criterion for dynamic configurations and loading:

∂E DP ∂α (u, ε p , p, α)•(β -α) ≥ 0. (106)
Finally, the plastic evolution is obtained by taking w=u and β=α. Then, for any p, we must have S DP (u, u, ε p , p, α) ≤ S DP (u, u, p, p -ε p + p, α), which is the same criterion used in the quasi-static case, that is, for any p

E DP (u, ε p , p, α) ≤ E DP (u, p, p -ε p + p, α). ( 107
)
The whole set of equations can now be written:

• Dynamic evolution:      ρü = divσ + f on Ω u = u D (t) on ∂ u σ•n = T (t) on ∂ T . ( 108
)
• Damage evolution: for any β≥0 admissible, we have

∂E DP ∂α (u, ε p , p, α)•(β -α) ≥ 0. ( 109
)
• Evolution of : for any p, we have

E DP (u, ε p , p, α) ≤ E DP (u, p, p -ε p + p, α). ( 110
)

Numerical Implementation

To calculate the evolution of the system, we have to solve three sets of equation:

1. Dynamic equation. This differential equation gives us the of the system, which can be integrated to obtain the velocity and the displacement. It is the only equation where there is a dependence in time and is the main issue.

2. Damage evolution. This is a partial differential equation that must be solved globally (because of the gradient of damage) each time.

3. Plastic evolution. Since this is a local problem, it must be solved on each element independently.

The can be represented by the instants t 0 < t 1 < t 2 , ... and we use standard finite elements discretization for spacial representation. The displacement and damage fields consist of polynomials of degree one (P 1 elements) and the , total strain and plastic strain are constant on each element.

To solve the dynamic equation, we will use a slightly modified version of the () scheme, since it is easy to implement, the total energy is conserved and every term can be obtained explicitly, resulting in fast calculations. The down side is that the condition forces the time ∆t to be small and decreases when we refine the mesh [START_REF] Allaire | Analyse numérique et optimisation: Une introductionà la modélisation mathématique et à la simulation numérique[END_REF].

Once we have a new displacement field u i at the instant t i , search for α and ε p that solve equations ( 109) and (110) simultaneously. For that, we use an between and damage. We notice that when E(•) and σ Y (•) depend on the same way of α, only one iteration is needed.

To study the evolution of the plastic strain, we can use the algorithms described in sections 2.2.2 and 2.2.3. It is, however, important to notice that σ Y is still a constant in the damage problem, since it depends only on α and α is fixed during this step.

We consider the displacement u i-1 , the velocity v i-1 , the plastic strain (ε p ) i-1 and damage α i-1 at the instant t i-1 to be known.

(1) Calculate the acceleration a i :

ρa i-1 = divσ = div E(α i-1 ):(ε(u i-1 ) -(ε p ) i-1 ) .
(2) update the displacement:

u i = u i-1 + ∆tv i-1 + ∆t 2 2 a i-1 .
(3) Set ((ε p ) (i,0) , α (i,0) ) := ((ε p ) (i-1) , α i-1 ).

(4) Iteration j ≥ 1:

(4.1) Using sections 2.2.2 and 2.2.3, solve

(ε p ) (i,j) := arg min p E DP (u i , p, pi-1 + p -(ε p ) i-1 , α (i,j-1) ).
(4.2) Solve α (i,j) := arg min

α i-1 ≤α≤1 E DP (u i , (ε p ) (i,j) , pi-1 + (ε p ) (i,j) -(ε p ) i-1 , α).
(4.3) Stop when (ε p ) (i,j) -(ε p ) (i,j-1) and α (i,j) -α (i,j-1) are sufficiently small.

(5) We then set ((ε p ) i , α i ) := ((ε p ) (i,j) , α (i,j) ) and pi =p + (ε p ) (i,j) -(ε p ) i-1 .

(6) Calculate the acceleration a i+1 using

ρa i+1 = divσ = div E(α i+1 ):(ε(u i+1 ) -(ε p ) i+1 ) .
(7) Update the velocity:

v i+1 = v i + ∆t 2 (a i + a i+1 ).

Applications

In this first set of examples, we want to show the influence of the plasticity and of the dynamics.

1-D Fracture

We first study a 1-D bar under traction. We have already studied this behavior for a quasi-static loading and we were able to break the bar, where damage is localized in one region. The bar was then split into two bars, and each bar had zero stress.

In , however, we have to take the into consideration. When the bar breaks, the bar is not is static equilibrium and there are that propagate in the bar and cause the bar to continue breaking, even after the first .

Another interesting phenomenon is the importance of the waves and vibrations. When a wave reaches one of the extremities, the wave is reflected and a compression stress becomes a traction stress and vice-versa. When two waves interpose each other, the resulting wave can have a greater amplitude than the initial waves.

We show here a simple examples that illustrates this behavior. Damage evolution will only be considered when the stress is positive and the material is assumed to be .

We will study the shock between two bars. The first bar, on the left, has size L and initial speed v 0 > 0. This bar hits a second bar, of size 2L, which is initially at rest. The two bars are made of the same material and have the same thickness.

A then propagates between through the bars causing a crack in the middle of the larger bar.

In our model, we consider a single bar of length 3L. The interval [0, L) has initial speed v 0 and (L, 3L] has zero initial speed. The stress is supposed to be zero at the extremities. The can be seen in Figure 5.

Material Behavior

In this section, we present a small and far from extensive list of material behaviors that we can obtain only by changing how the function E(•), w(•) and σ Y (•) depend on α. In these simulations, we suppose that the system evolves uniformly in space. In Figure 6, we have E(α)=(1 -α) 2 and w(α)=α and we don't have plastic strain. We can clearly see an and then a phase where damage evolves. By taking into account the (Figure 7), we see that we have now three phases (elastic, plastic with no damage and plastic with damage). It is important to notice that, for both models, the stress is maximal before the beginning of the damage phase and then it decreases until it reaches zero.

For this next set of models, where we take w(α)=α 2 , we see that the behavior changes. In Figure 8, we see the evolution of . There is no longer an elastic phase and, as strain increases, both damage and the stress increase, even though the relation stress-strain is no longer linear because of .

Many other evolution laws could be created by taking, for instance, a different polynomial degree for the previous expressions or by combining them.

Dimensionless Parameters

In this first example, we are going to study a bar made of a material under traction. The objective here is to show that the brittle damage model in question depends only on two .

We consider a bar Ω = [0, L]. We write E(α) = E 0 a(α), where a(α=0)=1 and a(α= 1)=0.

The study of the consists in defining an energy and finding its minimum with respect to α: The dynamic equation in 1-D can be written as

E(u, α) = Ω 1 2 E 0 a(α) ε(u) 2 + w 1 α + 1 2 w 1 2 (α ) 2 dΩ. (111)
ρü = σ = E 0 a(α)ε(u) (112)
and we impose a displacement on the extremities:

u(x=0, t) = 0 u(x=L, t) = ε0 Lt. (113)
For this model, we have the following parameters: L, E 0 , w 1 , , ρ and ε0 . The first step is to reduce the number of parameters of the problem. Since we are interested in finding the minimizer of E, it is clear that we can redefine E as

E(u, α) = Ω 1 2 E 0 w 1 a(α) ε(u) 2 + α + 1 2 2 (α ) 2 dΩ. (114)
It is clear that the dynamic equation depends only of E 0 ρ . Therefore, we are only interested in 5 values: L, E 0 w 1 , E 0 ρ , , and ε0 We will change the scale of our variables in order to remove 3 parameters from our problem. We first write x = 1 L x and t = T t, for some constant T > 0 that we will specify later.

If x ∈ Ω, then x ∈ [0, 1]. The imposed displacements are now u(x=0, t) = 0 u(x=1, t) = ε0 Lt. ( 115
) If f (x, t) is a function of x and t, we define f (x, t) := f (x, t). (116)
We derive it to obtain

df (x, t) dx = 1 L d f (x, t) dx and df (x, t) dt = 1 T d f ( t, t) d t . (117)
Suppose we have a constant U 0 > 0. We define ũ(x, t) :

= 1 U 0 u(x, t) α(x, t) := α(x, t). (118)
Thus

E(ũ, α) = Ω 1 2 E 0 U 2 0 w 1 L 2 a(α) dũ dx 2 + α + 1 2 2 L 2 dα dx 2 dΩ (119) and 1 T 2 d 2 ũ d t2 = E 0 ρL 2 d dx a(α) dũ dx . ( 120
)
Since we only assumed T and U 0 were two positive constants, we can now fix them. We set U 0 = L w 1 E 0 and T = L ρ E 0 . We also define ˜ := L and ε0 := ε0

T L U 0 . We have E(ũ, α) = Ω 1 2 a(α) dũ dx 2 + α + 1 2 ˜ 2 dα dx 2 dΩ (121)
and the dynamics of the system is

d 2 ũ d t2 = d dx a(α) dũ dx . (122)
Considering x and t as the space and time variables (and removing the tilde from out notation), we obtain the dimensionless problem in Ω = [0, 1] :

• The α minimizes the energy E taking into account the condition, where E is given by

E(u, α) = Ω 1 2 a(α) ε(u) 2 + α + 1 2 2 α 2 dΩ. (123)
• The time evolution of the displacement u is given by

ü = a(α)ε(u) , (124)
under the imposed boundary conditions

u(x=0, t) = 0 u(x=1, t) = ε0 t. ( 125
)
Example 4.1. Suppose we want to study the fracture of a bar made of steel. We suppose this bar is 10 cm long and is being stretched with a constant speed of 100m/s. For this material, we have a density ρ=8000kg/m 3 , modulus of elasticity E 0 =210GPa, fracture toughness K IC =50MPa•m 1/2 and an ultimate tensile strength σ c =1000MPa.

When using this damage gradient model in a quasi-static scenario, we know that

w 1 = σ 2 C E 0 = (1000 • 10 6 ) 2 210 • 10 9 Pa = 4.7MPa.
The fracture energy G c rate can be found:

G c = K 2 IC E 0 = (50 • 10 6 ) 2 210 • 10 9 = 11.7kPa • m.
The energy dissipated by the fracturing process, given by equation [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], can be written as

G c = σ C 4 √ 2 3
.

We obtain = 1.2 • 10 -5 m and ˜ = /L = 1.2 • 10 -4 . The values of T and U 0 are

T = L ρ E 0 = 0.1 8000 210 • 10 9 s = 2 • 10 -5 s U 0 = L w 1 E 0 = 0.1 4.7 • 10 6 210 • 10 9 m = 4.7 • 10 -4 m.
The deformation can be now be found:

ε = du dx = U 0 L dũ dx = 4.73 • 10 -3 dũ dx .
The dimensionless deformation speed is

ε0 = ε0 T L U 0 = 4.3 • 10 -3 ε0 = 0.43.
This bar can be simulated using our model with only two parameters ( ˜ = 1.2 • 10 -4 and ε0 = 0.43).

When analyzing the results, one must keep in mind that a time of 1 in the simulation is equivalent to 2 • 10 -5 s. In the same way, a deformation of 1 in the simulation is equivalent to ε = 4.73 • 10 -3 in the real bar.

1-D Period Bar

We are interested in obtaining the number of fragments of a ring under . Instead of working with a ring in a 3-D scenario, we consider a bar [0, L] and the following periodic conditions in the strain ε and damage α:

ε(x+L, t) = ε(x, t), x ∈ R α(x+L, t) = α(x, t), x ∈ R, (126)
for every t ∈ R.

We also suppose that we start our study with a completely healthy bar under uniform (in space) strain rate ε0 . At the instant t = 0, we have

ε(x, 0) = ε0 , x ∈ [0, L] α(x, 0) = 0, x ∈ [0, L]. (127)
For simplicity, we assume that the initial strain is zero, that is, ε(x, 0) = 0, x ∈ [0, L]. The ring of initial perimeter L suffers uniform . At an instant t, the perimeter of the ring is L+ ε0 Lt. The displacement u at the extremities and ε must satisfy

L 0 ε(x, t)dx = u(L, t) -u(0, t) = ε0 Lt. ( 128
)
It is clear that ε(x, t) = ε0 t satisfies the above equations.

At each time instant t, we write

ε(x, t) = ε * (x, t) + ε0 t. (129)
The variable ε * is the difference between the real and the uniform strain. It is easy to see that ε * is in x and

L 0 ε * (x, t) = 0. ( 130
)
We define the function u * as the difference between the real displacement and the uniform displacement, that is,

u * (x, t) = u(x, t) -ε0 xt. (131)
By differentiating the above equation, we find that (u * ) = ε * . It is clear that that u * is periodic. The stress can be written as

σ(x, t) = A(α)((u * ) + ε0 t). (132)
The dynamic equation is ü * (x, t) = ü(x, t) = σ .

Finally, we write the total energy as

E(u * , α) = L 0 1 2 A(α)((u * ) + ε0 t) 2 + w(α) + 1 2 w 1 2 (α ) 2 . ( 134
)
The problem consists of finding two variables u * and α satisfying the dynamic equation and dE dα (u * , α)β = 0 for every β admissible.

Influence of Each Parameter on Damage

We now want to investigate the influence of the parameters in the in the bar. As we have seen, the gradient damage model can be described by two parameters: the and the deformation speed ε0 .

There is, however, another parameter we need to pay attention to: the size of elements ∆x. Even though it is a purely numerical parameter, it is has a great importance in the number of fragments we are able to obtain. Therefore, for each pair ( , ε0 ) we want to study, we run several simulations with different values of ∆x and see if they number of fragments converge for ∆x → 0. We emphasize that we are interested in the convergence of the number of the fragments, and not in the convergence of u and α. Since we have a periodic problem, a translation of (u, α) would be a different numerical result, but the same physical result.

Numerical simulations show that there is only a small difference between the results if the mesh is fine enough. The seems to converge as ∆x → 0. In figure 11, we can see the influence of the mesh and the number of elements used. We recall that we used standard P1 elements. We note that for 1/ ≤ 1000, a simulation using 2000 elements (N elem = 2000) is accurate. For 1/ ≤ 2000, 5000 elements are enough. This holds true for other values of and N elem . Hence, we will consider that the results are accurate if 1/ ≤ 1 3 N elem or, equivalently, ≥ 3/N elem . For other values of ε0 , this relation also seems to be valid.

Again from figure 11, we can see that the number of fragments increases linearly with 1/ for ε0 = 0.5. This linear behavior holds true for every other value of ε0 between 10 -4 and 10 2 tested.

We can see the influence of ε0 in figure 12. We make ε0 vary between 10 -4 and 10 2 . The change from less than 500 to over 1000.

We show in figure 13 the results in a log-log scale. When ε0 ≥ 1, the points appear to be in a straight line.

We conclude with two remarks. Firstly, we can see that the behavior is not monotone. For close values of ε0 , we see that there is an oscillation in the number of fragments. But nevertheless, we see a clear tendency for the number of fragments to increase as the strain rate increases.

Secondly, we remark that a small initial perturbation does not change the in the end. There is, however, some changes in how they appear and how many of these develop until total failure, but we emphasize that the number of cracks in the end is the same if the perturbation is sufficiently small.

Cylinder under Internal Pressure

The final application concerns the study of the fragmentation of a cylinder under a strong internal pressure. We want to know when and how the material breaks. In this case, we expect the cylinder to fragment into multiple parts.

The cylinder has an internal radius of R i =1.0 and external radius of R e =1.25. We impose an internal pressure of 1.0 and we assume it to be constant in space and time. The following material We first study a brittle cylinder. We consider the dependency of the rigidity tensor with respect to α as E(α)=E 0 (1 -α) 2 , which we have already studied (Figure 1). The damage coefficients used in the energy are w 1 =8•10 -3 and the characteristic length =10 -2 . For these parameters, the critical stress is σ c = √ E 0 w 1 =8.9, which is a much larger than the pressure imposed. However, the internal pressure causes waves of stress to propagate and, by combining with each other, they reach the stress threshold and cause fracture on the cylinder.

The problem in question is axisymmetric and we could expect a radial symmetry of the results. In practice, a small perturbation is enough to make damage localize. One could argue that by changing the perturbation used, we could obtain different results. This is true; however, we are interested at the moment of fragmentation and the number of fragments, which don't seem to change when the perturbation is sufficiently small.

We see the damage profile in Figure 14. We can see the cracks in red and the healthy material in blue. There are two main aspects that we see in this figure: the first one is that we have several damage bands and that they are somewhat evenly distributed. The second is the direction of the cracks. We see that the cracks follow the radial direction almost everywhere.

Next, we consider a cylinder. We consider once again E(α)=E 0 (1-α) 2 and, now, σ Y (α)=σ 0 Y (1α) 2 (see Figure 4). We take the initial yield stress as σ 0 Y =5.0, such that it is greater than the applied pressure, but smaller than the critical damage stress. This way, we obtain an elastic phase, a phase where plastic strain occurs but not damage, and a phase with damage and plastic strain evolution.

The damage profile in Figure 15 is different from the one obtained for brittle materials. We can see once again that we have damage bands that are evenly distributed, but their direction has changed. As expected for ductile materials, when combining damage and plasticity, the cracks evolve in a 45 • angle.

The last question we want to discuss is the number of fragments. We recall that the characteristic length is proportional to the thickness of the damage band for a quasi-static loading. For a dynamic problem this is no longer true, but it is still the most important parameter for determining the number of cracks. If we decrease the characteristic length, since each crack takes less place, we should obtain more cracks.

To test this hypothesis, we run a simulation with a characteristic length 10 times smaller ( =10 -3 ). The result, Figure 16, is very clear. We can see that each crack is now thinner and we have many more cracks. We also notice that there are now many cracks that cross each other.

Conclusion

In this work, we have briefly explained the hypotheses considered in the construction of gradient damage models for brittle softening materials for a infinitely slow loading, based on the principle of . We then explained the necessary changes to take plastic strains and inertia into account.

We have shown that these models are very flexible, and only by changing how the rigidity and yield stress depend on the damage, we can model material behaviors that are very different.

We then concluded by some applications, showing the influence of each parameter of the model and by studying the fragmentation of a cylinder under initial pressure.

[1] Brodbeck D. Rafey R. A. Abraham, F. F. and W. E. Rudge. Physical review letters. Instability dynamics of fracture: a computer simulation investigation, 73(2):272, 1994.

Figure 1 :

 1 Figure 1: Damage (dashed red) and normalized stress (blue) for a 1-D bar subject to traction, according to the model given in Example 2.1.

Figure 2 :

 2 Figure 2: Damage profiles for w1=1 and E(α) = (1 -α) 2 depending on the maximal value of α, given by equation 30.

)

  If we write σ = E(ε i -p), then, by taking the absolute values, we obtain |σ| = |σ Y |. We can write e = σ * -σ E . (55) Since we are working on the case |σ| = σ Y < |σ * |, we have

Figure 3 :

 3 Figure 3: Damage (dashed green) and normalized stress (blue).

Figure 4 :

 4 Figure 4: Damage (dashed red), normalized stress (blue) and plastic strain (green dots).

Figure 5 :

 5 Figure 5: Damage profile after failure in a 1-D bar as described in subsection 4.1.

Figure 6 :

 6 Figure 6: Evolution of stress and damage for a 1-D brittle material for the case E(α)=(1 -α) 2 and w(α)=α; plastic strains neglected.

Figure 7 :

 7 Figure 7: Evolution of stress and damage for a 1-D ductile material (using Von-Mises criterion) for the case E(α)=σY (α)=(1 -α) 2 and w(α)=α.

Figure 8 :

 8 Figure 8: Evolution of stress and damage for a 1-D brittle material for the case E(α)=(1 -α) 2 and w(α)=α 2 ; plastic strains neglected.

Figure 9 :

 9 Figure 9: Evolution of stress and damage for a 1-D ductile material (using Von-Mises criterion) for the case E(α)=σY (α)=(1 -α) 2 and w(α)=α 2 .

Figure 10 :

 10 Figure 10: The displacement u * is at the top left and the damage profile is at the bottom right. We can also see the velocity distribution and the stress σ.

Figure 11 :

 11 Figure 11: We see that the number of cracks is almost proportional to 1/ if the mesh is fine enough. In here, ε0 = 0.5

Figure 12 :

 12 Figure 12: Influence of ε0 on the number of cracks for =2 • 10 -4 .

Figure 13 :

 13 Figure 13: Influence of ε0 for a fixed value of (log-log).

Figure 14 :

 14 Figure 14: Damage profile for a brittle cylinder under internal pressure. We consider E(α)=E0(1-α) 2 and w(α)=w1α.

Figure 15 :

 15 Figure 15: Damage profile for a ductile cylinder under internal pressure. We consider E(α)=E0(1-α) 2 , σY (α)=σ 0 Y (1α) 2 and w(α)=w1α.