The main objective of this work was the study of the fragmentation of a metallic shell. This thesis is divided into four parts: construction of a damage model, numerical implementation, calibration of the model parameters using experimental data and analytical works.

In this work, we considered a model that couples the standard gradient damage models with plasticity and dynamics. Using the energy and the action of the system, we could obtain all the equations necessary to describe the dynamic ductile model: the equations of dynamics, the plasticity criterion and the damage criterion. We then detail the numerical implementation of these models.

Some qualitative behaviours are then obtained, such as the number and the direction of cracks, and the convergence to the quasi-static model.

In order to better understand the influence of the parameters, we studied the problem analytically. By studying the amplitude of the perturbations, we describe how to obtain an analytic approximation for the number of cracks in a ring under expansion.

In order to run realistic simulations, it is needed to calibrate the material parameters. We focus here on a simple case of brittle materials. The experimental data were obtained in a series of shockless spalling tests performed by the CEA.

We also study other forms of regularization, now applied to the plastic strain, avoiding localization in zero-thickness bands. We considered using the dissipative properties of the temperature field to regularize the model. Finally, we conclude with plasticity models where we add a term depending on the gradient of the plastic strain (gradient plasticity models).

iii En premier lieu, j'aimerais remercier mon directeur de thèse, Jean-Jacques Marigo, d'avoir suivi mon travail durant ces trois années. J'ai apprécié sa gentillesse et sa pédagogie, ainsi que ses grandes connaissances en mécanique.

Je suis reconnaissant à mes rapporteurs d'avoir lu avec attention mon manuscrit et à mon jury de thèse d'avoir fait le déplacement pour ma soutenance. Merci à eux pour leurs questions et remarques très pertinentes.

Un grand merci au CEA Saclay et au CEA Gramat qui ont permis la réalisation et le bon déroulement de cette thèse en finançant mes recherches ; en particulier, Daniel Guilbaud et Gilles Damamme pour leur encadrement, ainsi que Julien Grunenwald, Skander El Mai et Jean-Lin Dequiedt qui m'ont suivi pendant ces années.

Résumé

Cette thèse porte sur l'étude de la fragmentation d'enveloppes métalliques avec des applications dans le domaine militaire. L'enveloppe est mise en expansion par la détonation d'explosifs et la très forte pression (quelques centaines de kilo-bars) ainsi générée. L'état de contrainte induit dans le matériau va conduire à sa fragmentation et à la génération d'un très grand nombre d'éclats. Le principal objectif de cette thèse est de prévoir le nombre, la forme et la distribution massique de ces fragments.

Cette thèse est divisée en cinq parties : la construction d'un modèle d'endommagement, l'implémentation numérique, des études analytiques, la calibration des paramètres du modèle en utilisant des données expérimentales, et des travaux analytiques.

Tout d'abord, nous avons considéré des modèles qui couplent les modèles d'endommagement classiques avec la plasticité et la dynamique. En utilisant l'énergie et l'action du système, nous avons obtenu toutes les équations qui décrivent le modèle dynamique et ductile : l'équation de la dynamique, le critère de plasticité et le critère d'endommagement.

Nous avons ensuite détaillé l'implémentation numérique de ces modèles. Deux codes ont été utilisés : la bibliothèque d'éléments finis FENICS et le logiciel EUROPLEXUS. Dans un premier instant, nous avons implémenté les modèles d'endommagement avec la bibliothèque FENICS pour des tests initiaux, en particulier pour des problèmes unidimensionnels. Ensuite un des modèles de fracture ductile a été implémenté dans le code industriel EUROPLEXUS, avec lequel nous avons fait des simulations des problèmes tridimensionnels.

En ce qui concerne la performance du code, le problème d'endommagement peut être écrit comme un problème linéaire, où la matrice en question est symétrique définie-positive. Par conséquent, nous avons pu utiliser la méthode du gradient-conjugué, déjà implémentée dans la bibliothèque PETsC, et qui marche très bien dans les codes parallélisés. v Dans un premier instant, des résultats qualitatifs ont pu être obtenus, comme le nombre et la direction des fissures, ainsi qu'une étude de la convergence vers le modèle quasi-statique.

La principale application est l'explosion d'un cylindre métallique à cause d'une forte pression intérieure. De façon surprenante, le problème d'un cylindre avec un chargement radial perd sa symétrie et nous obtenons plusieurs fissures inclinés qui se croisent. Une première question que se pose c'est de comprendre pourquoi ces zones de localisation de l'endommagement apparaissent.

Afin de mieux comprendre l'influence de chaque paramètre du modèle, nous avons fait des études analytiques. Le problème du cylindre a été simplifié en un anneau, qui peut être vu comme une barre avec des conditions aux limites périodiques. A partir de l'observation de l'amplitude des perturbations, nous avons pu décrire comment obtenir une approximation analytique du nombre de fissures pour l'anneau en expansion.

Cependant, pour être capable de simuler des problèmes réalistes, il est nécessaire de calibrer les paramètres du modèle. Nous nous sommes intéressés plus particulièrement au problème d'écaillage de matériaux fragiles (céramiques). A partir des données expérimentales obtenues par une série d'expériences réalisée par le CEA, nous avons pu calibrer les paramètres de notre modèle pour avoir une bonne approximation de l'énergie dissipée par le processus de rupture.

Des travaux complémentaires ont également été réalisés concernant l'écaillage et la modélisation de la striction. Afin d'empêcher la localisation de la déformation plastique dans des bandes d'épaisseur nulle, d'autres formes de régularisation ont été étudiées, comme par exemple, l'utilisation des propriétés dissipatives du champ de température. Enfin, nous avons conclu ce travail en proposant des modèles de plasticité où l'énergie dépend aussi du gradient de la déformation plastique (modèles de plasticité à gradient). D'une forme générale, les travaux effectués pendant cette thèse ont aidé à mieux comprendre l'évolution de l'endommagement dans un contexte dynamique. Sur le plan numérique, ces modèles marchent, peuvent être parallélisés et donnent des bonnes directions de fissures. La fragmentation d'un cylindre sous pression a été étudiée en 1D et 3D et l'influence de chaque paramètre du problème a pu être identifiée. Comme continuation de ce travail, nous avons encore deux grandes questions théoriques : la convergence vers le modèle quasi-statique et l'épaississement des régions endommagées.

vi Introduction

The initiation and propagation of cracks is still an unresolved question in fracture mechanics. Several models have been studied in different contexts (Barenblatt [9], Abraham and Rudge [START_REF] Brodbeck | Physical review letters. Instability dynamics of fracture: a computer simulation investigation[END_REF], Hentz and Daudeville [START_REF] Donzé | Discrete element modelling of concrete submitted to dynamic loading at high strain rates[END_REF], Hakim and Karma [START_REF] Hakim | Laws of crack motion and phase-field models of fracture[END_REF]), in quasi static and dynamics (Ravi-Chandar [START_REF] Ravi-Chandar | Dynamic fracture of nominally brittle materials[END_REF], Larsen [START_REF] Christopher | Models for dynamic fracture based on griffith's criterion[END_REF]), and accounting for different phenomena. The main objective of this thesis is to explain the development of the so-called "gradient damage models" (Pham [START_REF] Amor | Gradient damage models and their use to approximate brittle fracture[END_REF]) and its extension to ductile materials under a dynamic loading.

The main application behind this thesis is that of a metallic shell that expands due to a strong internal pressure, until it fragments. Several models have been proposed to estimate the number and size of the resulting fragments. These models focus mostly on the one-dimensional expanding ring and use statistical arguments or presence of micro-voids, as in Mott and Linfoot [START_REF] Mott | A theory of fragmentation[END_REF], Grady [START_REF] Grady | Fragmentation of Rings and Shells: The Legacy of N.F. Mott[END_REF]. Our approach differs from the previous ones in the sense that we consider a homogeneous and sound material, and no random phenomenon is considered.

The idea behind the models used in damage mechanics is that we can represent the crack by a damage field. No a-priori hypothesis of its path is made.

As we will see in a simple example, local models are not capable of correctly predicting damage evolution (Peerlings and Brekelmans [START_REF] Geers | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF], Pham et al. [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]). Softening local damage models allow damage localization in infinitely thin bands and, consequently, cracks with zero energy dissipation (Benallal [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. Moreover, some attention must be payed to the mesh in order to avoid creating a preferential direction for crack propagation (Negri [38]).

In this context, the problem of localization is solved by adding a nonlocal Contents term, such as an integral (Pijaudier-Cabot and P. Bazant [START_REF] Pijaudier | Nonlocal damage theory[END_REF], P. Bazant and Pijaudier-Cabot [START_REF] Zdenek | Nonlocal continuum damage, localization instability and convergence[END_REF], Peerlings and de Vree [START_REF] De Borst | Gradient enhanced damage for quasi-brittle materials[END_REF], Lorentz and Andrieux [START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF]) or a gradient (Comi [16], Lorentz and Benallal [START_REF] Lorentz | Gradient constitutive relations: numerical aspects and application to gradient damage[END_REF], Lorentz et al. [START_REF] Lorentz | Convergence of a gradient damage model toward a cohesive zone model[END_REF]) of the damage or the strain. The family of gradient damage models contain the gradient of damage weighed by a parameter called the "characteristic length" (Pham et al. [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with Bibliography 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 brittle damage evolution, but have also been extended to ductile (Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], Ambati et al. [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], Miehe et al. [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-elasticplastic solids[END_REF]) and dynamic loading (Bourdin et al. [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF], Borden et al. [START_REF] Michael | A phase-field description of dynamic brittle fracture[END_REF], Li et al. [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]). In this thesis, we explain the necessary changes to the original model, in order to take plastic deformations and inertial effects into account.

In the first chapter, we briefly present the construction of gradient damage models for brittle softening materials based on the principle of minimum energy. We discuss the main hypothesis and the need for regularization. We then talk about the Von-Mises plasticity criterion, how to write it using a principle of minimum energy and how to couple plasticity and damage by using a suitable form of energy, as done by Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF]. We conclude the model by removing the hypothesis of static equilibrium at each instant and adding inertial effects. We follow the same methodology of Li et al. [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]: we write the Lagrangian and the action of the system, and find the equations of dynamics, along with the criterion of damage and plasticity evolution, by using the principle of least action.

In Chapter 2, once the model is complete, we detail how the evolution of damage and plasticity is calculated numerically, and the schemes used for the temporal integration. In a first stage, we consider the standard Lagrangian discretization using P1 elements. We then show some examples to validate our implementation, test the convergence rate in function of the mesh size and the time-step and have a first insight in how the parameters of the problem affect the results. We conclude this chapter by detailing the implementation of the discontinuous Galerkin (DG) methods for quasi-static and dynamic damage-plasticity problems. The FEniCS library (Logg [START_REF] Mardal | Automated Solution of Differential Equations by the Finite Element Method -The FeniCS Book[END_REF]) and the industrial code EUROPLEXUS are used.

In Chapter 3, we study the particular case of a cylinder in expansion. After fragmentation, we want to count the number of fragments obtained and determine how it depends on the parameters used. The problem is axially symmetric and, therefore, we should obtain an axially symmetric profile for the damage. Surprisingly, this is not what happens, as we obtain radial cracks somewhat evenly spaced. In order to understand what is causing the evolution of these cracks, we focus mostly on the one-dimensional case, that is, a ring. By studying the linearised system, we show that some modes of Chapter 1

Dynamic Gradient Damage Models

The objective of this first chapter is to explain the development of the so-called "gradient damage models" and its extension to ductile materials and dynamic loading. The main idea of these models is that a crack can be represented by a scalar (the damage field). No hypothesis are made a-priori of its path.

We explain how gradient damage theory deals with the question of damage localization in infinitely thin bands (and, consequently, cracks that dissipate no energy) by adding a term containing the gradient of damage weighed by a parameter called the "characteristic length". Roughly speaking, this constant is going to determine how thick the crack region is going to be. These models have been originally proposed for quasi-static brittle damage evolution, but have also been extended to ductile (Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], Ambati et al. [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], Miehe et al. [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-elasticplastic solids[END_REF]) and dynamic loading (Bourdin et al. [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF], Borden et al. [START_REF] Michael | A phase-field description of dynamic brittle fracture[END_REF], Li et al. [START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF]). The main objective of this chapter is to explain the necessary changes in order to account for both plastic deformations and inertial effects.

More precisely, we first present the construction of gradient damage models for brittle softening materials based on the principle of energy minimization. We discuss the main hypothesis and present one example in order to illustrate the need for regularization. We then briefly talk about Von-Mises plasticity criterion and how to take it into account. We conclude the model by moving from quasi-static to dynamic loadings.

Chapter 1 -Dynamic Gradient Damage Models

Gradient Damage Models

We present here a simplified construction of gradient damage models for brittle elastic materials when there are no other dissipative phenomena. We are going to consider the case of small strains theory and isotropic material.

For a more detailed construction of these models, see Marigo [START_REF] Marigo | L'endommagement et la Rupture : hier, aujourd'hui et demain[END_REF], Bourdin et al. [START_REF] Bourdin | The variational approach to fracture[END_REF], Pham [START_REF] Pham | Construction et analyse de modèles d'endommagement à gradient[END_REF], Pham et al. [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with Bibliography gradient damage models[END_REF]. For the proof of Gamma-convergence to Griffth's model (Griffith [21]), the reader is referred to Braides [START_REF] Braides | Γ-convergence for beginners[END_REF], Dal-Maso and Toader [START_REF] Dal-Maso | A model for the quasi-static growth of brittle fractures: Existence and approximation results[END_REF].

We denote the stress by σ, the strain by ε, the displacement by u and the rigidity tensor by E. When working in a 1D scenario, we are going to call the Young's modulus simply by E, the stress by σ and the strain by ε.

We recall that ε = 1 2 (∇u + ∇ T u). It is clear that we consider that the variables in question are regular enough so that trace and the energies are well-defined. Unless otherwise stated, the variables will be either in the L 2 (Ω) or the H 1 (Ω) spaces. The contracted product of two tensors a and b will be denoted by a:b and, for elastic materials, the stress can be written σ = E:ε.

Construction of a Damage Model (non-regularized)

In this section, we are going to describe a family of damage models that can be applied to different types of materials. We will discuss the qualitative properties of these models.

We begin the construction by making the following hypothesis:

1. Damage can be represented by a scalar α∈[0, 1]. When α=0 the material is sound and when α=1 the material is completely broken.

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

3. Damage is irreversible, that is, it can only grow in time ( α≥0).

We now need to specify under which circumstances damage increases. For that, we are going to use an idea similar to Griffith's criterion (Griffith [21]), based on the notion of elastic energy restitution, in its variational form (Francfort and Marigo [18]).

The elastic energy can be written as ψ(ε, α) = 1 2 ε:E(α):ε.

(1.1)

-Gradient Damage Models

For a fixed deformation, a small increase δα>0 of damage causes a loss of elastic energy equivalent to -∂ψ ∂α (ε, α)δα>0. 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(α) (1.2)
where k(α)≥0 is a function of α representing the necessary energy restitution necessary for damage to evolve. Let w(α) be a function such that w (α) = k(α). We define the energy density by W (ε, α) = ψ(ε, α) + w(α).

(1.3)

We can write the stress as

σ = ∂W ∂ε (ε, α) (1.4)
and the damage evolution criterion as ∂W ∂α (ε, α)• α = 0, (1.5) 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. (1.6)
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 in Ω u = u D on ∂ u σ•n = T on ∂ T . (1.7)
We fix a test function w such that w = 0 on ∂ u . Then

Ω divσ•w + f •w dΩ = 0 (1.8)
Chapter 1 -Dynamic Gradient Damage Models and Green's formula shows that

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

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

Ω ∂W ∂ε (ε(u), α):ε(w)dΩ = Ω f •wdΩ + ∂T T •wdS, ∀w ∈ C 0 . (1.11)
If we consider the evolution problem where the time is denoted by t, by integrating (1.6), we obtain the following problem: find α≥0 such that

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

E(u, α) = Ω W (ε(u), α)dΩ - Ω f •udΩ - ∂T T •udS. (1.13)
It is easy to see that the evolution problem, given by equations (1.11) and (1.12), is equivalent to finding u ∈ C and α≥0 such that

DE(u, α)(v -u, β -α) ≥ 0, ∀v ∈ C, ∀ β ≥ 0. (1.14)

Regularized Model

It is now a well-known fact that local softening damage models are not viable (Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], Pham [START_REF] Amor | Gradient damage models and their use to approximate brittle fracture[END_REF]) as they allow damage localization in infinitely thin bands. The example below illustrates this problem: Example 1.1.1. Consider a 1D bar represented by the interval [0, 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).

-Gradient Damage Models

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 .
(1.15)

Therefore, damage in this interval is given by

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

D = L 0 w(α)dx = L θL w 1 α * dx = w 1 α * (1 -θ)L (1.17)
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.

The main idea behind gradient damage models is to add to the energy a term that depends on the gradient of damage. This way, sharp damage profiles will dissipate an infinite amount of energy and will not be minimizers of this energy. This construction leads to the notion of a characteristic length of the damage problem. We will now use an energy density of the form

W (ε, α, ∇α) = ψ(α, ε) + w(α) + 1 2 w 1 2 ∇α•∇α, (1.18)
where is the characteristic length 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 minimum energy. 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 notice that for a homogeneous damage profile, we obtain the same damage criterion.

We have

σ = E(α):ε = ∂W ∂ε (ε, α). (1.19)

Chapter 1 -Dynamic Gradient Damage Models

We define the dissipated energy by

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

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

DE(u, α)(v -u, β -α) ≥ 0, ∀v ∈ C, ∀ β ≥ 0. (1.22)

Damage Coupled with Plasticity

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

For that, we first review the plasticity model that we are going to use, showing how it can be written as a problem of energy minimization. We then discuss the model obtained when writing an energy functional that contains both plasticity and damage dissipation terms.

In this thesis, only the Von-Mises criterion will be studied, even though only a few adaptations are needed if we want to consider other criteria. More details about the coupling gradient damage and plasticity can be seen in Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF], Tanne [START_REF] Tanne | Variational phase-field models from brittle to ductile fracture: nucleation and propagation[END_REF].

We finish this section by showing some examples of material behaviour that can be obtained using this approach.

Perfect Plasticity Model Unidimensional Model

In this section, we follow the approach of Marigo [START_REF] Marigo | Plasticité et Rupture[END_REF], Alessi et al. [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF]. We will denote the plastic strain by ε p . The total strain can be decomposed in an elastic part (a part that contributes to the stress) and a plastic part (a permanent strain). The stress-strain relation is now

σ = E(ε -ε p ).
(1.23)

-Damage Coupled with Plasticity

Figure 1.1: Damage (dashed green) and normalized stress (blue) for a generic ductile material.

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. (1.24)
We are interested in the Von-Mises yield criterion, where σ Y is the yield stress. In 1D, we have

f (σ) = |σ| -σ Y ≤ 0. (1.25)
This behaviour is shown in Figure 1.1. We can see the normalized strain, normalized stress and normalized plastic strain.

We define the cumulated plastic strain from zero to the instant t as

p(t) = t 0 | εp |dτ (1.26)
and the energy density for a elasto-plastic material as

W 1D (ε, ε p , p) = 1 2 E(ε -ε p ) 2 + σ Y p. (1.27) Proposition 1.2.1.
The 1D version of the Von-Mises criterion, written as

1. yield criterion: |σ| ≤ σ Y 2. flow rule:      εp ≥ 0 if σ = +σ Y εp = 0 if |σ| < σ Y εp ≤ 0 if σ = -σ Y (1.28)
Chapter 1 -Dynamic Gradient Damage Models is equivalent to 1. stability condition: for any p * , we have

W 1D (ε, ε p , p) ≤ W 1D (ε, p * , p + |p * -ε p |) 2. energy balance: Ẇ 1D (ε, ε p , p) = σ ε.
Proof. We first notice that for any p * , we have

W 1D (ε, p * , p+|p * -ε p |)-W 1D (ε, ε p , p) = 1 2 E(p * -ε p ) 2 -E(ε-ε p )(p * -ε p )+σ Y |p * -ε p |.
(1.29) If the yield criterion holds, it is easy to see that the stability condition is also true.

Conversely if the stability condition holds and by taking p * →ε p , we obtain

-E(ε -ε p )(p * -ε p ) + σ Y |p * -ε p | ≥ 0. (1.30)
If we divide it by p * -ε p and study the cases p * ≥ε p and p * ≤ε p , we obtain the yield criterion.

We take the derivative of W 1D :

Ẇ 1D (ε, ε p , p) = σ( ε -εp ) + σ Y | εp |. (1.31)
If the flow rule holds, it is easy to see that the the energy balance is also verified.

Finally, if the energy balance holds, then

σ εp = σ Y | εp |. (1.32)
Thus, if |σ|≤σ Y , then εp =0. Otherwise εp and σ have the same sign.

By using this, it is clear that the plasticity criterion can be written as the minimization an energy defined as the integral of W 1D .

Three-dimensional Model

In 3D, the stress-strain relation can be written as

σ = E:(ε -ε p ).
(1.33)

The Von-Mises criterion is now given by the function

f (σ) = 3 2 s:s -σ Y , (1.34)
where s:=σ -Trσ 3 I is the deviatoric stress. We also recall that this criterion imposes that Trε p =0.

We can define the energy density as

W 3D (ε, ε p , p) = 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y p (1.35)
and, by following the same steps described in 1D, we can show that calculating the evolution of plasticity in 3D is equivalent to minimizing this energy.

Damage-Plasticity Coupling

In this section, in order to construct a family of models that account for plasticity and damage, instead of proposing the evolution laws for each variable, we work directly with a suitable form of energy and, by minimizing this energy, we deduce the constitutive relations. For simplicity, we remove volume forces from our calculations

We recall that, in section 1.1.2, we obtained a total energy for brittle damage:

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

In section 1.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Ω (1.37)
in 1D, and

E 3D plast (ε, ε p ) = Ω 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y p dΩ (1.38)
in 3D. By examining perturbations in ε and ε p , obtain the static equilibrium and the plasticity criterion, respectively. As we can see, the problems of damage and plasticity are similar in the sense that the quasi-static evolution in both cases is found after minimizing the total energy. For the coupled problem, we are going to use an energy Chapter 1 -Dynamic Gradient Damage Models 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 yield stress now depends on the damage, that is, σ Y =σ Y (α).

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

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

E 3D DP (u, ε p , p, α) = Ω 1 2 (ε(u)-ε p ):E(α):(ε(u)-ε p )+ 2 3 σ Y (α)p+w(α)+ 1 2 w 1 2 |∇α| 2 dΩ.
(1.40) To obtain the quasi-static evolution criterion, we minimize the total energy with respect to all three variables (u, ε p and α):

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

(1.41)

• The minimization of the plastic strain gives us

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

1 2 E(α )(ε(u) -ε p ) 2 + σ Y (α)p + w(α ) -w 1 2 α ≥ 0 (1.43)
In 3D:

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

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

E(α) = E 0 (1 -α) 2 , w(α) = w 1 α , σ Y (α) = σ 0 Y (1 -α) 2 . (1.45)
Since we are assuming uniformity in space, we only 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=ε p =ε-σ 0 Y /E 0 . The damage criterion becomes -(1 -α) (σ 0 Y ) 2 E 0 -2(1 -α)σ 0 Y p + w 1 ≥ 0. (1.46)
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 plastic evolution 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. (1.47)
We can thus find

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

Material Behaviour

In order to illustrate the reach of such models, we show some examples of material behaviour that we can obtain only by changing how the function E(•), w(•) and σ Y (•) depend on α. The curves were obtained considering homogeneous damage, as in example 1.2.2.

In Figure 1.3, we have E(α)=(1 -α) 2 and w(α)=α and we don't have plastic strain. We can clearly see an elastic phase and then a phase where damage evolves. By taking into account the plastic evolution (Figure 1.4), 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 behaviour changes. In Figure 1.5, we see the evolution of brittle damage. 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 damage evolution.

-Material Behaviour

-Dynamic Damage Models

The list of models described above is, of course, far from extensive. Many other evolution laws could be created by taking, for instance, a different polynomial degree for the previous expressions or by combining them. It is important to notice that when we take functions E, w and σ Y that depend linearly or quadratically on α, the damage problem is linear, which is a easier to calculate, specially numerically. If E(α)/σ Y (α) is constant for every α, then the plasticity problem depends only on the strain, and not on the damage.

Dynamic Damage Models

To formulate the evolution of the dynamic system, we are going to use the principle of least action, as in Li [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 . (1.49)
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Ω (1.50)
and Green's formula shows that

Ω ρü•wdΩ = ∂u (σ•n)•wdA 0 + ∂T T •wdA - Ω σ(u):ε(w)dΩ + Ω f •wdΩ.
(1.51) 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.
(1.52)

Chapter 1 -Dynamic Gradient Damage Models

We define the kinetic energy of the system

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

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

t 2 t 1 ∂P ∂u w - ∂K ∂ u ẇ dt = 0, ∀w. (1.55)
We have thus shown that the problem (1.49) implies equation (1.55). It is easy to see that equation (1.55) 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.
(

This motivates us to construct a dynamic gradient model by defining a suitable form of the action functional. Instead of using a purely elastic energy σ:ε, we are going to use the energies defined by equations (1.39) and (1.40) with the terms containing the plastic strain and energy dissipated by the damage process. We remember 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 potential energy as We define the admissible displacement space C and admissible damage space D by

P DP (u, ε p , α) = E DP (ε(u), ε p , p, α) - Ω f •udΩ - ∂ T T •udA. ( 1
C = {u : u(t)=u 0 (t) on ∂ u } D = {α ∈ [0, 1] : α ≥ 0 on Ω} (1.60)
In order to preserve the irreversibility 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, β) (1.61)
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 (1.62)
and, by following the previous calculations in reverse order, we find the problem given by (1.49). We now set w=u and p=ε p to study the damage evolution. 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 Ω}.
(1.63)

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

∂E DP ∂α (u, ε p , p, α)•(β -α) ≥ 0. (1.65)
Finally, we look at the plastic evolution 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, α), Chapter 1 -Dynamic Gradient Damage Models 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, α).

(1.66)

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 .
(1.67)

• Damage evolution: for any β≥0 admissible, we have

∂E DP ∂α (u, ε p , p, α)•(β -α) ≥ 0.
(1.68)

• Evolution of plastic strain: for any p, we have

E DP (u, ε p , p, α) ≤ E DP (u, p, p -ε p + p, α). (1.69)
Chapter 2

Numerical Implementation and Validation

In this chapter, we detail the numerical implementation of the gradient damage model using the finite elements method. We consider a spatial discretization based on the standard Lagrange family of P 1 elements, unless otherwise stated. We first discuss the damage problem, calculations of the plastic strain and dynamics, showing the algorithms and numerical methods used for each separate problem. We then show some test cases to validate our implementation and we discuss some qualitative properties of the dynamic damage model.

For the time discretization, we consider the instants t i , with t i+1 =t i + ∆t. In 1D, the elements of the mesh have the same length ∆x. In 2D and 3D, we will specify whether we are using a structured mesh or an unstructured mesh.

We finish this chapter by detailing the implementation of the discontinuous Galerkin (DG) methods. We write the variational formulation associated to it and how the type of element used affects our results.

For the calculations, we used the FeniCS (Logg [START_REF] Mardal | Automated Solution of Differential Equations by the Finite Element Method -The FeniCS Book[END_REF]) library and the industrial code Europlexus [START_REF] Team | Europlexus: a computer program for the finite element simulation of fluid-structure systems under transient dynamic loading[END_REF].

Implementation of Plasticity

We present here the algorithm we used when calculating the evolution of the plastic strain. We will consider that the total deformation ε is known and fixed, and we are only interested in the evolution of ε p . Even though this Chapter 2 -Numerical Implementation and Validation algorithm is commonly used in solid mechanics, we considered important to detail it here.

The only important remark here is that this algorithm considers the yield stress σ Y to be constant. When coupling plasticity and damage, the only necessary change for the algorithm is to consider the yield stress as a function of damage.

Plasticity in 1D

We first discuss how the plasticity was implemented. We remember from section 1.2.1 that the evolution of the plasticity can be found by minimising W 1D defined by

W 1D (ε, ε p , p) = 1 2 E(ε -ε p ) 2 + σ Y p. (2.1)
It is important to notice that this is a local problem, that is, it can be solved independently in each element or Gauss point.

Suppose that the plastic strain is (ε p ) i . We define the auxiliary function

f (ε, p) = 1 2 Ep 2 -Eεp + σ Y |p -(ε p ) i-1 |. (2.2)
In the discrete problem, 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.1.1. 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|. (2.3) (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 .

-Implementation of Plasticity

f (p) = f (ε, (ε p ) i-1 ) + 1 2 Eh 2 -σ * h + σ Y h. (2.4)
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.1.2. 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 . (2.5) (2) If |σ * | > σ Y : (ε p ) i = (ε p ) i-1 + 1 - σ Y |σ * | ε i -(ε p ) i-1 (2.6)
and

E(ε i -(ε p ) i ) = σ Y . ( 2
(ε i , p) = Ee -σ * + σ Y e |e| = 0. (2.9)
Hence,

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

E(ε (i,j) -p) = σ Y e |e| . ( 2
) i := p = (ε p ) i-1 + e = (ε p ) i-1 + 1 -σ Y σ Y |σ * | (ε i -(ε p ) i-1 ). (2.15)
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 (elastic prediction). If the stress is inside the elastic domain, that is, |σ * |<σ Y , we keep it and the plasticity does not change. On the other hand, if the stress is not in the elastic domain, we update the plastic deformation (plastic correction).

Plasticity in 3D

In this section, we are going to write the same results as in the previous section, but now to a problem in 3D. We use the standard Von Mises criterion:

tr(ε p ) = 0 (2.16) and 3 2 s : s ≤ σ Y , (2.17)
where the deviatoric stress tensor s is given by s = σ -tr(σ) 3 I.

(2.18)

-Implementation of Plasticity

The discrete energy density at the instant i + 1 is written as

W (ε, ε p ) = 1 2 (ε -ε p ):E:(ε -ε p ) + 2 3 σ Y ε p -(ε p ) i + 2 3 σ Y pi
(ε i -(ε p ) i-1
) and its deviatoric part is given by s * =2µ(e -(ε p ) i-1 ).

The following propositions are the 3D equivalents of the auxiliary results in section 2.1.2: Proposition 2.1.3. 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 * :δ.
(2.23)

Chapter 2 -Numerical Implementation and Validation

(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 . (2.24) 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.1.4. The minimization of W in ε p is equivalent to:

(1) If s * ≤ 2 3 σ Y : (ε p ) i = (ε p ) i-1 . (2.25) (2) If s * > 2 3 σ Y : (ε p ) i = (ε p ) i-1 + 1 - 2 3 σ Y |s * | e i -(ε p ) i-1 . (2.26)
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. (2.27) If σ i = E:(ε i -p) (2.28)
and s i is its deviatoric part, we must have

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

s i = 2 3 σ Y .
(2.30)

-Implementation of Damage

We note that

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

s i = 2 3 σ Y s * s * . (2.32)
Finally, applying this to (2.29) and (2.31),

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
).

(2.33)

We conclude by taking (ε p ) i = p.

Implementation of Damage

The implementation of damage using the FEniCS library is straight forward: we define the total energy of the system and find its derivative with respect to α and in the direct β using the derivative(energy, alpha, beta) command.

We solve the resulting constrained minimisation problem using the class OptimisationProblem along with the NonlinearVariationalSolver.

The main advantage of this approach is that, once the code is implemented, studying the influence of the functions E(α), w(α) and σ Y (α) demands little effort in terms of programming.

We now detail the implementation of the damage problem in the industrial code EUROPLEXUS. We use the model

E(α)=(1 -α) 2 E 0 , w(α)=w 1 α and σ Y (α)=(1 -α) 2 σ 0
Y . The damage problem consists of finding α ∈ [α min , α max ] that minimizes the total energy, that is

Ω (α -1)ε el Eε el β + w 1 β + w 1 2 ∇α∇β + 2(α -1)σ Y pβ ≥ 0,∀β. (2.34)
We have to solve a linear system on the form Aα = b, where

A = (ε el Eε el + 2σ Y p) + w 1 2 ∇ T ∇ (2.35) and b = ε el Eε el + 2σ Y p -w 1 . (2.36)
The first step is the initialization of our variables and of the libraries used (PETSC and TAO).

Initialization:

• Create the table containing the value of α at each node.

• Store the determinant of the jacobian J 0 on each element at t = 0.

• Initialize PETSC and TAO.

• Assemble the constant matrix A 0 := w 1 2 ∇ T ∇J 0 .

• A 0 (i, j) = A 0 (j, i) = δ ij if the material of the node i cannot be damaged.

• Set the vectors α min = 0 and α max = 1

At each time step i + 1:

• update the value of α min :=α i ;

• assemble the matrix A 1 := ε el Eε el + 2σ Y p;

• assemble vector b and the matrix A = A 0 + A 1 ;

• find the vector α i+1 by solving Aα = b using TAO;

It is important to recall that the solution of this last linear system using TAO also takes into account the irreversibility condition. The GPCG solver (J. Moré and Toraldo [START_REF] Moré | On the solution of large quadratic programming problems with bound constraints[END_REF]) is used. Others solvers, such as the Scalable Nonlinear Equations Solvers (Balay et al. [START_REF] Balay | PETSc users manual[END_REF]) were tested, but with a less satisfying performance in our problems.

Dynamic Numerical Schemes

In this section, we describe the schemes used to solve the dynamic evolution. We first discuss each method for an elastic material and then we add its extension to damage and plasticity.

We are going to detail three schemes: the explicit Newmark scheme, an implicit variational scheme and the generalized midpoint rule scheme, studying the influence of vibrations, convergence rate, energy dissipation and stability for each one.

-Dynamic Numerical Schemes

Explicit Newmark Scheme

This method is used to solve second order (linear or non-linear) differential equations. It is commonly used in civil engineering for numerical evaluation of the dynamic response of structures. In particular, the industrial code we used for the large calculations uses this scheme, so we first investigate its properties and how to couple it with damage.

The finite elements method gives us a system of the form

M Ü (t) + KU (t) = f (t).
(2.37)

We will approximate U (t), U (t) and Ü (t) by the sequences

U i , U i and Ü i satisfying      U i+1 = U i + ∆t U i + (∆t) 2 2 Ü i U i+1 = U i + ∆t 2 ( Ü i + Ü i+1 ) M Ü i+1 = f (t i+1 ) -KU i+1 . (2.38)
It is easy to see that this implies that

(U i+1 -U i ) = (U i -U i-1 ) + ∆t( U i -U i-1 ) + (∆t) 2 2 ( Ü i -Ü i-1 ) U i -U i-1 = ∆t 2 ( Ü i-1 + Ü i ).
(2.39)

Therefore U i+1 -2U i + U i-1 = (∆t) 2 Ü i . (2.40)
We obtain the equivalent equation for the Newmark scheme

M U i+1 -2U i + U i-1 (∆t) 2 + KU i = f (t i ). (2.41)
We can see that this scheme is a second order scheme (in time). It is now a well-known fact that this scheme is stable if

max j λ j (∆t) 2 < 4, (2.42)
where λ j are the eigenvalues of the problem KU = λMU .

In particular, for a uniform 1D mesh with elements of length ∆x, we have that the eigenvalues of the laplacian are smaller than 4/(∆x) 2 . For a material of density ρ and Young's modulus E, we obtain thus the following CFL condition ∆t ∆x ≤ ρ E .

(2.43)

We now have to change the matrix K because of damage and the force f to include the effects of plasticity. In order to avoid confusion, we will use the Chapter 2 -Numerical Implementation and Validation notation u for the displacement, v for the velocity, a for the acceleration , α for the damage and ε p for the plastic strain. When written in the variational form, form a test function w, we have

Ω ρawdV + Ω E(α)ε(u)ε(w)dV = Ω E(α)ε p ε(w)dV (2.44)
which, can be written as

Ma + K(α)u = f (α, ε p ). (2.45)
We propose the following algorithm:

(1) Update boundary conditions.

( and plastic strain (ε p ) i+1 .

) Calculate u i+1 = u i + ∆tv i + (∆t) 2 2 a i . 2
(4) Find the acceleration

Ma i+1 = f (α i+1 , (ε p ) i+1 ) -K(α i+1 )u i+1 .
(5) Find the velocity v i+1 = v i + ∆t 2 (a i + a i+1 ). (6) Advance to time step i + 1.

As we'll see in the validation section, this method produces good results with an almost-constant energy. The main problem we face here is the apparition of vibrations.

We also notice that when the yield stress and rigidity tensor have the same dependency on the damage, the result of the plasticity problem is independent of the damage. For this reason, the damage and the plasticity problem will be solved only once, allowing us to gain a significant amount of time at teach iteration.

Implicit Variational Scheme

The objective of this subsection is to propose a simple, intuitive and easy to implement scheme that allows us to compare results using a different discretization.

-Dynamic Numerical Schemes

Contrary to the previous section, we now propose a method for the dynamic equation that should be solved, along with the problems of damage and plasticity, until convergence of all variables.

For this scheme, we use the same ideas of Bourdin et al. [START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF]. We use the following approximation for the second derivative of u:

üi ≈ u i -2u i-1 + u i-2 ∆t 2 . (2.46)
The problem of dynamics becomes finding the displacement u i+1 solution of

Ω (ε(u i+1 )-(ε p ) i+1 ):E(α i+1 ):ε(w)dV = ρ Ω u i+1 ∆t 2 •wdV -ρ Ω 2u i -u i-1 ∆t 2 •wdV ,
(2.47) for any test function w admissible.

We show that this scheme is dissipative for the elastic problem. In fact, the problem can be written as

Ω ε(u i+1 ):E:ε(w)dV = ρ Ω v i+1 -v i ∆t •wdV , (2.48)
where

v i = u i -u i-1 ∆t . ( 2

.49)

By taking w=v i+1 ∆t=u i+1 -u i , we find

Ω ε(u i+1 ):E: ε(u i+1 ) -ε(u i ) dV = ρ Ω (v i+1 -v i )•v i+1 dV , (2.50)
We remark that

(v i+1 -v i )•v i+1 = 1 2 v i+1 2 - 1 2 v i 2 + 1 2 v i+1 2 + 1 2 v i 2 -v i •v i+1 = 1 2 v i+1 2 - 1 2 v i 2 + 1 2 v i+1 -v i 2 .
(2.51)

This shows us that

( 1 2 ε(u i+1 ):E:ε(u i+1 ) + ρ 2 v i+1 2 ) -( 1 2 ε(u i ):E:ε(u i ) + ρ 2 v i 2 ) = - 1 2 (ε(u i+1 -ε(u i )):E:(ε(u i+1 -ε(u i )) - ρ 2 v i+1 -v i 2 .
(2.52)

Chapter 2 -Numerical Implementation and Validation

We see that the total energy between two time steps always decreases if there are no sources of energy to the system. From these calculations, we can expect this change of energy to be proportional to ∆t 2 at each time step.

In conclusion, this method is very dissipative (dissipation proportional to ∆t), and the numerical solution is more regular than the analytic solution (this becomes very evident when studying a problem containing shockwaves). This regularizing effect, even though undesirable in certain cases, can be helpful when analysing the behaviour of displacement and velocity near the cracks, as will be discussed later.

Even though we don't have stability problems, we still have to pay attention to the time step in order to control how much energy is dissipated numerically.

Generalized Midpoint Rule Scheme

In this section, we detail a 1D algorithm that can be seen as a generalization of the previously proposed variational scheme. This scheme and the calculations presented could easily be extended to 2D and 3D, but we chose to focus on the 1D case, as we wanted to try a third method to obtain the qualitative behaviour of our model.

We start by presenting the scheme and the calculations in Simo and Hughes [START_REF] Simo | Computational Inelasticity. Interdisciplinary applied mathematics[END_REF] for an elastic-plastic problem.

We fix η ∈ [0, 1], θ ∈ [0, 1] and δ ∈ [0, 1], and define

f (σ) = |σ| -σ Y . (2.53)
For a variable X, we define the notation

X θ = (1 -θ)X i + θX i+1 . (2.54)
We have the following equations for an elastic-plastic material:

               u i+1 -u i = ∆tv η Ω ρ (v i+1 -v i ) ∆t w + Ω σ θ w = 0, ∀w ε p i+1 -ε p i = p sign(σ θ ), p ≥ 0 f (σ δ ) ≤ 0 p • f (σ δ ) = 0. (2.55)
To obtain the stability, we take the test function

w = u i+1 -u i = ∆tv η . Therefore Ω ρ(v i+1 -v i )v η + Ω σ θ (u i+1 -u i ) = 0.
(2.56)

We have

(v i+1 -v i )v η = (v i+1 -v i ) (1-η)v i +ηv i+1 = ηv 2 i+1 -(1-η)v 2 i +(1-2η)v i v i+1 . (2.57) Since v i v i+1 = 1 2 -(v i+1 -v i ) 2 + v 2 i+1 + v 2 i , (2.58)
we have

(v i+1 -v i )v η = 1 2 v 2 i+1 -v 2 i -(1/2 -η)(v i+1 -v i ) 2 .
(2.59)

From this general expression, we can study the properties for each particular case.

Purely Elastic Case

When there is no plastic strain, we have

Ω ρ(v i+1 -v i )v η + Ω Eu θ (u i+1 -u i ) = 0.
(2.60)

We can then see that the difference of total energies (elastic energy plus kinetic energy) between two different instants is

Ω ρ 2 v 2 i+1 + 1 2 E(u i+1 ) 2 - Ω ρ 2 v 2 i - 1 2 E(u i ) 2 = Ω ( 1 2 -η)(v n+1 -v i ) 2 +( 1 2 -θ)E(u i+1 -u i ) 2 .
(2.61) From that, we obtain the the scheme is table if η≥1/2 and θ≥1/2 and non-dissipative if η=θ=1/2.

Elastic-Plastic Case

For an elastic-plastic material, the relation

Ω ρ(v n+1 -v n )v η + Ω Eσ θ (u n+1 -u n ) = 0 (2.62)
still holds. We have then

σ θ (u i+1 -u i ) = σ θ 1 E (σ i+1 -σ i ) + p sign(σ θ ) = 1 E σ θ (σ i+1 -σ i ) + p|σ θ |.
(2.63) On one hand, we have

σ θ (σ i+1 -σ i ) = 1 2 σ 2 i+1 -σ 2 i -(1/2 -θ)(σ i+1 -σ i ) 2 .
(2.64)

Chapter 2 -Numerical Implementation and Validation

On the other,

p|σ θ | = p•f (σ θ ) + pσ Y . (2.65) This implies that Ω ρ 2 v 2 i+1 + 1 2E (σ i+1 ) 2 - Ω ρ 2 v 2 i - 1 2E (σ i ) 2 + Ω pσ Y = Ω (1/2 -η)(v i+1 -v i ) 2 + (1/2 -θ) E (σ i+1 -σ i ) 2 -pf (σ θ ).
(2.66)

We use the fact that pf (σ δ )=0. Therefore, if θ=δ, then the last term disappears.

The scheme is stable and non dissipative if θ=η=δ=1/2.

Adding Damage

We now try to obtain similar results for the case where damage is considered. The difficulty appears from the fact that E is no longer constant between iterations. There is also a question of how to define σ θ from the other variables. For instance, we can take

     σ 1 θ = (1 -θ)σ i + θσ i+1 σ 2 θ = E(α θ )(u θ -ε p θ ) σ 3 θ = (1 -θ)E(α i ) + θE(α i+1 ) (u θ -ε p θ ) = E θ (u θ -ε p θ ).
(2.67)

As in the calculation of plasticity, we are free to choose γ∈[0, 1] and define at which moment between t i and t i+1 we want the damage criterion to be satisfied, leading to the following system of equations

                   u i+1 -u i = ∆tv η Ω ρ (v i+1 -v i ) ∆t w + Ω σ θ w = 0, ∀w ε p i+1 -ε p i = p sign(σ θ ), p ≥ 0 f (σ θ ) ≤ 0 p • f (σ θ ) = 0 α γ := argmin α≥α i E(u γ , ε p γ , pγ , α).
(2.68)

We propose then the following algorithm:

-Dynamic Numerical Schemes

(1) Solve the displacement problem:

u i+1 -u i = ∆tv 1/2 Ω ρ (v i+1 -v i ) ∆t w + Ω σ 1/2 w = 0, ∀w (2.69)
(2) Solve the plasticity problem at the instant θ=1/2.

(3) Solve the damage problem at the instant γ.

(4) Check the convergence of u i+1 , (ε p ) i+1 and α i+1 . Repeat if it did not converge, otherwise go to the next iteration.

From the previous calculations, it is natural to consider θ=1/2. Numerical experiments have shown that for all of these choices of σ θ or γ, the system will have the same behaviour when the time-step and mesh size are small enough.

From a theoretical point of view, we can show that by taking σ θ =E θ (u θ -ε p θ ), we obtain a scheme that is stable. In fact, as in the previous calculation, we take the test function

w= ∆t 2 (v i +v i+1 )=u i+1 -u i . We have Ω ρ 2 (v i+1 -v i )(v i + v i+1 ) + Ω σ 1/2 (u i+1 -u i ) = 0. (2.70)
It is clear that the first term is equivalent to the change in the kinetic energy. We will then focus on the second term. We write E i :=E(α i ) and obtain

σ 1/2 (u i+1 -u i ) = 1 4 (E i + E i+1 )(u i+1 -ε p i+1 + u i -ε p i ) • (u i+1 -ε p i+1 ) -(u i -ε p i ) + (ε p i+1 -ε p i ) = 1 4 (E i + E i+1 ) (u i+1 -ε p i+1 ) 2 -(u i -ε p i ) 2 + 1 4 (E i + E i+1 )(u i+1 -ε p i+1 + u i -ε p i )(ε p i+1 -ε p i ) = 1 2 E i+1 (u i+1 -ε p i+1 ) 2 - 1 2 E i (u i -ε p i ) 2 + 1 4 (E i -E i+1 ) (u i+1 -ε p i+1 ) 2 + (u i -ε p i ) 2 + + σ 1/2 (ε p i+1 -ε p i ).
(2.71)

If we write the dissipated energy as

D = Ω 1 4 (E i -E i+1 ) (u i+1 -ε p i+1 ) 2 + (u i -ε p i ) 2 + σ 1/2 (ε p i+1 -ε p i ), (2.72)
Chapter 2 -Numerical Implementation and Validation we obtain, as in the previous cases, that the difference between the kinetic plus the elastic energy between two instants changes of -D, that is,

Ω ρ 2 v 2 i+1 + Ω 1 2 E(α i+1 )(u i+1 -ε p i + 1) - Ω ρ 2 v 2 i + Ω 1 2 E(α i )(u i -ε p i ) +D = 0.
(2.73) To show that this scheme is stable, we only have to show that D≥0. In fact, as in the ductile case, we know that σ 1/2 (ε p i+1 -ε p i )≥0. Using the irreversibility of the damage, we know that E i ≥E i+1 , proving that D≥0 and, therefore, that the scheme is stable.

Numerical Verification

In this section, we show some test cases in order to verify our numerical implementation. One of the difficulties found when testing the model was the absence of simple analytic solutions that could be used to test all phenomena at the same time.

The first test was the homogeneous solution found in example 1.2.2. Due to the simplicity of this case, the results will not be shown in this section.

In order to test the convergence rate and dissipative properties of the proposed algorithms, we ran a series of simulations for fixed material parameters, and we compared the results for different numerical parameters.

The third test we used was the study of a rigid-plastic bar subjected to damage. For this case, it is possible to construct an analytic solution for the dynamic damage problem. The numerical implementation cannot be completely rigid, since the dynamics of the system depends on the elasticity. For this reason, we ran simulations using very large values for the rigidity.

Finally, as a last test, we propose an algorithm in order to obtain the evolution of damage for a quasi-static loading. The main difficulty of this case is that we want to impose the dissipated energy, and not the displacement. Since this problem is hard to be solved numerically, it will allow us to test if the damage profile can be found in a reliable way.

Rate of Convergence

In this section, we will detail the tests used to validate our numerical implementation. More precisely, we use these tests in order to verify rate of convergence, stability and energy conservation of each scheme.

Since we want to test the quadratic convergence in time of these schemes, we want a non-homogeneous problem such that the second derivative in time of the solution is continuous in time. When testing plasticity evolution and

-Numerical Verification

damage localization, we want to avoid border effects, so we want maximum stress to happen far from the extremities. Based on these considerations, we propose the problem below on the interval [0, 1] for the displacement u: initial conditions:

u(t=0) = 0 v(t=0) = -π cos(πx); (2.74)
boundary conditions:

u(x=0) = -sin(πt) u(x=1) = sin(πt). (2.75)
We then solve the equation ü=u numerically. For this 1D problem, it is clear that the analytic solution is u(x, y)= -sin(πt) cos(πx), but we prefer calculating the error of our simulations in a purely numeric way, as this approach is also suitable for the cases with damage and plasticity. The approach remains simple: if we want to test the error using a mesh and a time-step, we will first run a simulation where the time-step and the mesh are a lot finer than the ones in question and consider this result as the exact solution of the problem.

We first test the elastic Newmark scheme. When taking ∆t=∆x, we can verify the superconvergence property (Hughes [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]), where we obtain errors of order 10 -14 even for a coarse mesh.

For this first series of tests, we will take ∆t=0.5∆x in order to respect the stability condition and avoid the superconvergence case. For the Genezalized Midpoint Rule scheme, the mesh and the time-step are independent, and superconvergence does not occur. Since we are interested in the convergence in time, we take a mesh fine enough so that the errors caused by the discretization in space are negligible when compared to errors produced by the temporal discretization. The errors are obtained using the L 2 norm.

In order to test the convergence rate and dissipative properties of the proposed algorithms, we ran a series of simulations for fixed material parameters, and we compared the results for different numerical parameters. We can see in the quadratic convergence in Figure 2.1 for both schemes.

When comparing both schemes, it is important to keep in mind that the error in the Genezalized Midpoint Rule scheme comes exclusively from temporal discretization, whereas in the Newmark scheme the error originated from the spatial discretization is non-negligible, since the time-step and meshsize are of the same order. Since the Newmark scheme is explicit, it is the one who runs the fastest between the two, when using the same mesh and time-step. We also verified that the energy is conserved for both schemes. The next step is to add damage and plasticity to the schemes. We will use the same initial and boundary conditions given by equations 2.74 and 2.75.

For both schemes, adding damage to the system does not change the convergence rate. We consider E(α)=(1-α) 2 When considering plasticity (with or without damage), we obtain a clear loss of the quadratic convergence. Even if we try changing the instant be-

-Numerical Verification

tween t i and t i+1 when plasticity is calculated, we do not obtain a better convergence rate. One possible explanation for this is that the plasticity depends on the gradient of u, which cannot be obtained with enough precision with the type of finite elements that we are using.

In summary, the two schemes we implemented have a quadratic convergence rate for purely elastic and brittle damage problems, but a poorer convergence rate when plasticity is considered.

For the Newmark scheme, when there is only damage or only plasticity, or when the yield stress and rigidity tensor have the same dependency on the damage, we solve each problem only once, making this scheme very fast in each iteration.

For the Genezalized Midpoint Rule scheme, since we have to solve each problem a few times before converging, along with the fact that the dynamic problem is implicit, each iteration takes more time. This fact is, however, counterbalanced by the fact that there are no restrictions on the time-step.

In terms of stability, numeric experiments show that both schemes can be considered to conserve the total energy, but we have no theoretical results supporting our claim. For the Genezalized Midpoint Rule scheme, we proved that the scheme is always stable.

Rigid-Plastic Bar

We consider a bar of density ρ and [-L, L]. We suppose that its extremities are being pulled with constant speed ε0 L > 0. We also suppose the initial speed is uniform. We also suppose that the bar is not damaged at t=0 and the damage phase is imminent.

We first construct an analytic solution to this problem and then we compare it to the numerical results. We consider

σ P (α) = (1 -α)σ 0 P (2.76) and w(α) = w 1 α. (2.77) Dynamic equation: σ -ρ v = 0. (2.78) Compatibility equation: v -ε = 0. (2.79)
Since there is no elastic deformation and the movement is always in the same direction, we can consider that the accumulated plasticity is equal to the total deformation: ṗ = | ε|.

(2.80) The boundary conditions can be written as

     v(x = -L) = -ε0 L v(x = L) = ε0 L α(x = -L) = α(x = L) = 0 . (2.81)
The initial condition are

     v(x, t = 0) = ε0 x p(x, t = 0) = ε(x, t = 0) = ε 0 α(x, t = 0) = 0 . (2.82)
Since we are about to start the damage phase, we have

ε 0 = w 1 σ 0 P . (2.83)
We define the parameter

c 2 0 = σ 0 P ρε 0 . (2.84)
We search for a symmetric solution and we divide the bar in two zones: in zone (1), the bar plastifies and is damaged; in zone (2), the bar is damaged but there is no evolution in plasticity. The boundary between the two regions is expressed by X(t) or T (x) and x = X(t) ⇔ t = T (x).

Region (1):

By hypothesis, there is plastic and damage evolution in this region.

-Numerical Verification

By (2.76) and (2.78), ρ v = σ = -α σ 0 P and we obtain

α = - v c 2 0 ε 0 and α = -ε c 2 0 ε 0 . (2.85)
From the damage criterion, we have

0 = w 1 -σ 0 P ε -w 1 2 α . ( 2
α (x, t) = A(x) cosh c 0 t + B(x) sinh c 0 t (2.89)
We have α = 0 for t = 0, so A(x) = 0. Differentiating (2.86) in time, we find

v = ε = -ε 0 2 α . (2.90)
The initial conditions also give us

v 0 = ε0 = -ε 0 2 α (t = 0) = -c 0 ε 0 B(x). Therefore, α (x, t) = - ε0 c 0 ε 0 sinh c 0 t (2.91) and v = ε0 cosh c 0 t . (2.92)
Using again (2.86), ε = ε 0 (1 -2 α ) and

ε = ε 0 + ε0 c 0 sinh c 0 t (2.93) Region (2):
By hypothesis, there is only damage evolution in this zone. We only consider here the points x ≥ 0.

Since there is no plastic change, the speed of all points is the same:

v = ε0 L.
Differentiating (2.86) in time, we find α = 0, α = F (t) and α = F (t)(L -x).

(2.94)

Therefore, α can be written as

α = G(t)(L -x) + H(x).
(2.95)

Chapter 2 -Numerical Implementation and Validation

In the boundary:

In region (1), equation (2.92) gives us v 1 . In region (2), we know that v 2 = ε0 L. By continuity, we must have v 1 (X(t), t) = v 2 (X(t), t). Thus,

v 1 = ε0 cosh c 0 t (x -X(t)) + ε0 L. (2.96)
We know the (plastic) deformation in ( 1) is given by (2.93).

Since the plasticity in (2) does not change in time, we have

ε 2 (x, t) = ε 2 (x, T (x)) = ε 1 (x, T (x)) and
ε 2 = ε 0 + ε0 c 0 sinh c 0 T (x) (2.97)
Since region ( 2) is being damaged, (2.86) is valid and

2 α 2 = 1 - ε 2 ε 0 . (2.98)
By (2.97) and (2.98),

α 2 = - ε0 c 0 ε 0 sinh c 0 T (x) .
(2.99) Therefore α 2 = 0. We differentiate (2.91) in time to find

α 1 = - ε0 ε 0 2 cosh c 0 t . (2.100)
We take x → 0 in (2.96) to find

v 1 (0 + , t) = ε0 L -X(t) cosh c 0 t . (2.101)
We integrate (2.86) in order to obtain the jump in 0:

[u] + 2 ε 0 [α ] = 0 (2.102)
and, by symmetry,

α (0 + , t) = - u(0 + , t) 2 ε 0 . (2.103) Thus α 1 (0 + , t) = - v 1 (0 + , t) ε 0 2 = - ε0 ε 0 2 L -X(t) cosh c 0 t (2.104)

-Numerical Verification

and, since α 2 is constant,

α 2 = α 2 (L, t) = α 1 (0 + , t) + L 0 α = α 1 (0 + , t) + α 1 X(t) = - ε0 L ε 0 2 . (2.105) We imposed α(L, t) = 0, which implies α2 (x, t) = - ε0 L ε 0 2 (x -L).
(2.106)

We can find

α 1 (x, t) = α 1 (0 + , t)+ x 0 α 1 = - ε0 ε 0 2 L- ε0 ε 0 2 cosh c 0 t x-X(t) . (2.107)
We integrate the above expression and, by continuity of α in X(t),

α1 (x, t) = ε0 L ε 0 2 (L -x) - ε0 2ε 0 2 cosh c 0 t x -X(t) 2 .
(2.108)

We remember that v 1 is given by (2.96). We calculate

v1 = ε0 c 0 sinh c 0 t (x -X(t)) -ε0 cosh c 0 t Ẋ(t). (2.109)
We use (2.85) to obtain

α 1 = - v1 c 2 0 ε 0 = - ε0 2 ε 0 cosh c 0 t (x-X(t))+2 ε0 c 0 ε 0 sinh c 0 t Ẋ(t)+ ε0 c 2 0 ε 0 cosh c 0 t Ẍ(t).
(2.110) From (2.107) and (2.110), we get the differential equation

cosh c 0 t Ẍ(t) + 2c 0 sinh c 0 t Ẋ(t) = - c 2 0 L 2 . (2.111)
Since the whole bar is undergoing plastic deformation at t=0, we must have X(0)=L.

At t=0, σ is uniform and v1 =σ =0. Since all points continue to plastify, X does not change and Ẋ(0) = 0.

We can now solve (2.111):

X(t) = L cosh( c 0 t ) . (2.112)
We find the inverse of this function to obtain T (x):

T (x) = c 0 cosh -1 L x . ( 2

.113)

Chapter 2 -Numerical Implementation and Validation

Equation (2.108) becomes α1 (x, t) = ε0 L 2 ε 0 2 - ε0 2ε 0 2 cosh c 0 t x 2 - ε0 L 2 2ε 0 2 cosh c 0 t -1 . (2.114)
Since α 1 (t=0)=0, we can find

α 1 (x, t) = ε0 L 2 c 0 ε 0 - x 2 2L 2 sinh c 0 t + c 0 t -arctan tanh c 0 t 2 . (2.115)
We can also find σ 1 = (1 -α 1 )σ 0 P . We showed in (2.106) that α2 does not depend of t. Therefore

α 2 (x, t) = α 1 (x, T (x)) + α2 t -T (x) = α 1 (x, T (x)) + ε0 L ε 0 2 (L -x) t -T (x) .
(2.116) In region (2), σ 2 = v2 =0. We can find the stress in this region by

σ 2 (x, t) = σ 1 x, T (x) .
(2.117)

ANALYTIC SOLUTION

We have the two regions separated by:

X(t) = L cosh( c 0 t ) or T (x) = c 0 cosh -1 L x .

Region (1):

This region is defined by the pair (x, t) such that |x| < X(t)

α 1 (x, t) = ε0 L 2 c 0 ε 0 -x 2 2L 2 sinh c 0 t + c 0 t -arctan tanh c 0 t 2 σ 1 (x, t) = (1 -α 1 )σ 0 P ε 1 (x, t) = ε 0 + ε0 c 0 sinh c 0 t .

Region (2):

This region is defined by the pair (x, t) such that

X(t) < |x| ≤ L α 2 (x, t) = α 1 (x, T (x)) + ε0 L ε 0 2 (L -x) t -T (x) σ 2 (x, t) = σ 1 x, T (x) ε 2 (x, t) = ε 0 + ε0 c 0 sinh c 0 T (x) .

-Numerical Verification

Numerical solution

We solved this problem numerically and compared the numerical solution to the analytic one.

The numerical schemes implemented so far work only for elastic bars and, by hypothesis, we have a rigid-plastic bar. This issue was solved by setting the Young's modulus large enough so that the elastic deformation is negligible compared to the plastic deformation. We are going to consider the Young's modulus as E(α)=(1 -α) q E 0 , for q ≥ 0.

In the rigid bar, p = ε and at t = 0, we have ε = w 1 /σ 0 P . For the elastic bar, we must take into account this difference:

p(x, 0) = w 1 σ 0 P - qσ 0 P 2E 0 ε(x, 0) = p(x, 0) + σ 0 P E 0 .
(2.118)

Numerical results

We obtained good results numerically and when we refine the mesh and increase E 0 , the results are also more precise.

Chapter 2 -Numerical Implementation and Validation

Final remarks

• Numerically, we obtained the most accurate results for q=1, that is,

E(α)=(1 -α)E 0 .
• Increasing the value of E 0 , we increase the amplitude of oscillations for σ. We can solve this problem by refining the mesh.

• We also needed to refine the mesh in order to find a smooth profile for X(t).

• The stability condition for the Newmark scheme requires ∆t ≤ ∆x √ E 0 . Since we need big values for E 0 and we need a fine mesh, the time step is small and the simulations take a long time to run.

Control of Dissipated Energy

As a test for the code used, we propose a numerical algorithm to obtain the snapback curve. Our goal is to find the evolution of damage when the dissipated energy is controlled and compare it to the usual case, where we consider the displacement at one of the extremities as the loading parameter. This problem is hard to be solved numerically and, thus, will allow us to see the robustness of our code and the solvers used.

We consider a brittle 1D bar Ω = [0, L] under a quasi-static loading u(0, t) = 0 and u(L, t) = t. We are going to use the AT1 model, that is, E(α)=E 0 (1 -α) 2 and w(α)=w 1 α.

As we already know, this bar under traction deforms elastically for ε < ε c := w 1 /E 0 and breaks instantly after. As a result, σ and the energies are discontinuous in time.

We remember that the dissipated energy is

D(α) = Ω w(α) + 1 2 w 1 2 (α ) 2 dΩ, (2.119)
the elastic energy is

E el (u, σ, α) = Ω 1 2 E(α)(u ) 2 dΩ -σu(L) (2.120)
and the total energy is

E(u, σ, α) = E el (u, σ, α) + D(α). (2.121)
More precisely, in this chapter we want to find an admissible triplet (u, σ, α) such that we can control the value of D(α):

Chapter 2 -Numerical Implementation and Validation

Control of Dissipated Energy

For a fixed displacement u and scalar t, we want to solve the problem inf α≥0:D(α)=t E(u, α).

(2.128)

We define the Lagrangian for the constraint D(α) = t: For a given value of λ, we can easily calculate the above minimization in α (and thus find G(λ)) using the library TAO.

L(α, λ) = E(u, σ, α) + λ(D(α) -t). ( 2
We have now to find the maximum with respect to λ. We take ∆λ sufficiently small and we calculate the approximations to the derivatives of G:

G = G(λ+∆λ)-G(λ) ∆λ G = G(λ+2∆λ)-2G(λ+∆λ)-G(λ) (∆λ) 2 (2.132) 50

-Numerical Verification

The point of maximum satisfies G (λ) = 0. We could use Newton's algorithm to solve this problem:

λ n+1 = λ n - G G . (2.133)
This approach, however, does not work well numerically. As we can see in Figure 2.7, G is almost a straight line everywhere, except for a small region.

We notice that we can write the Lagrangian as

L(α, λ) = E el (u, σ) + (1 + λ)D(α) -λt. (2.134)
The minimization of the Lagrangian with respect to α is similar to our usual damage problem, where the dissipated energy is multiplied by 1+λ. As a result, for large values of λ, α = 0 (we do not reach the critical displacement to damage) and D(α) = 0. Therefore, the Lagrangian varies linearly with λ in this region. Hence

d 2 G dλ 2 = 0.
(2.135)

This shows that we cannot use Newton's algorithm for large values of λ. Numerical simulation have also shown convergence problems if λ is too small.

We chose to use a brute-force algorithm. It is clear that this approach is not optimal, but it works reasonably well and gives us good results.

To find the maximum of G, we first fix the non-empty interval [A 0 , B 0 ] and N ≥ 2. We use the iterative algorithm, starting at n = 0:

• We set ∆ = (B n -A n )/N .
• Set λ n as the argument tested which gives the largest value of G:

λ n = arg max 0≤i≤N G(A 0 + ∆ • i).
(2.136)

• Find α n such that L(α n , λ n ) = inf α≥0 L(α, λ n ).
(2.137)

• Stop if D(α n ) = t. • Otherwise      (A n+1 , B n+1 ) = (λ n -∆, λ n + ∆) if A n < λ n < B n decrease A n , B n = λ n + ∆ if λ n = A n A n = λ n -∆, increase B n if λ n = B n .
(2.138)

The above algorithm searches for the largest value of G in an interval that decreases with each iteration. It stops when it reaches the maximum, that we know is when D(α n ) = t. The concavity of G assures the convergence.

For a fixed displacement u, we have described how we can find α such that D(α) = t. We describe now the main algorithm.

We set σ 0 = σ c = E 0 ε c and α 0 = 0.

For n ≥ 1:

• Find the displacement u * u * = arg min u E(u, α n-1 , σ n-1 ) (2.139)
• Find a saddle point (α * , λ) of the Lagrangian (u * is the displacement used implicitly in the definition of L). Thus D(α * ) = t.

• Find the displacement u n :

u n := arg min u E(u, α * , σ n-1
).

(2.140)

• Find the damage profile α n :

α n := arg min α≥0 E(u n , α, σ n-1
).

(2.141)

• Find the new value of the imposed stress σ n as the mean stress on the bar:

σ n := Ω E(α n )ε(u n )dΩ Ω dΩ
.

(2.142)

• Stop if α n = α * .

Analytic Solution

We are now going to briefly remember construction of a localized solution in the interval [x 0 -D, x 0 + D], where x 0 is the center of the damage profile (assumed to be symmetric) and the value of D is unknown for now.

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

1 2 A (α)ε 2 + w (α) -w 1 2 α = 0. (2.143) 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. (2.144) 52

-Numerical Verification

We multiply this expression by α and we integrate to obtain

- 1 2 S(α)σ 2 + w(α) -w 1 2 α 2 2 = constant = - σ 2 2E(0) . ( 2

.145)

We define the function H by

H(α) := σ 2 1 w 1 E(0) - S(α) w 1 + 2w(α) w 1 .
(2.146)

Then 2 α 2 = H(α).
(2.147)

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) . ( 2

.148)

We can then find the damage profile using the relation

x(α) = x 0 αmax α H(β) dβ.
(2.149)

Results

We now compare the analytic and the numerical results. As we can see in Figure 2.8, our algorithm gives us the correct damage profile when comparing to the analytic solution (dashed red lines).

Figure 2.9 shows the different profiles we obtain. In the legend, we see the maximum value of damage. It is important to remark that the support of the damage profile decreases as damage increases. This, however, does not contradict the irreversibility condition as, in a quasi-static simulation, the cracks appear instantly without passing through the intermediate stages.

Finally, Figure 2.10 shows the stress in function of the displacement during the snapback. The analysis begins at σ=1, when damage begins and follows the decrease in the stress.

Overall, we obtained a purely numeric method to obtain the snapback phenomenon. We impose the dissipated energy at each instant and, from that, we obtain the damage profile, as well as the stress and the strain. We validated our results by a comparison with the analytic solution.

Qualitative Behaviour

In this section, we will briefly discuss the results obtained in our simulations. As for material behaviour, we considered E(α)=E 0 (1-α) 2 , w(α)=w 1 α and, when plastic strains are considered, σ Y (α)=σ 0 Y (1-α) 2 .

-Qualitative Behaviour

Multiple Cracks

When studying the rupture of a bar under a quasi-static loading, we obtain one single damage profile that appears instantly. When dynamics are considered, the damage profile increase gradually. If the internal length is small enough, we can obtain multiple cracks. These cracks do not appear necessarily at the same time, they depend on the waves propagating in the system.

We also see that the average damage in the bar changes when we use different strain rates: when we increase the strain rate, we also increase the average damage, as can be seen when comparing Figures 2.11 and 2.12.

For smaller strain rates, we obtain less cracks. This leads to the question of what would happen when ε→0. Unfortunately we cannot decrease ε enough in our simulations to know if there will be only one damage profile.

Convergence to Quasi-Static

When comparing the results obtained for a quasi-static loading and for a dynamic loading, the damage profiles are not the same. We consider a 1D bar, and we impose α=0 at the extremities. In quasistatic, the displacement at the borders is controlled. In dynamics, we impose uniform strain rate ε in the bar and we study what happens when ε→0. We see that the damage profile converges, but to a new profile that dissipates more energy (between 25% and 30% more) than the one in quasi-static, as shown in Figure 2.13. We have found two factors that contribute to this difference. The first one is the irreversibility condition. It is known that the intermediate damage profiles in the family of AT1 models have larger support than the final configuration, when α=1 (see, for instance, Figure 2.9). The quasi-static crack happens instantly and, therefore, it does not go through the intermediate profiles. Since the dynamic crack appears gradually, it goes through the intermediate stages and, due to the irreversibility condition, the final stage will be larger than the one in quasi-static. By removing the irreversiblity condition, we have verified that, in fact, its contribution in the difference of profiles is minor.

The main factor seems to be the snapback phenomenon. When it occurs, it is now well known that there is a difference between the energy dissipated by the damage process and the energy gathered in the elastic phase. As a result, there is a discontinuity in time of the total energy. In the dynamic model, almost all the energy is used to damage the structure, leading to a larger damage profile. To test this hypothesis, we consider a test case where there is no snapback. As shown in Pham [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], this can be obtained if the internal length is large enough. Therefore, we consider =0.5. In this case, we obtain in Figure 2.14 almost the same profiles for quasi-static and dynamic loadings.

Thickening of Cracks

If we continue the running a simulation after the value of a damage profile has reached the value of 1, the cracks continue to evolve. We notice that in a crack, α=1 in an interval that increases with time (see Figure 2.15) and the displacement profile has oscillations in this interval. We observed this phenomenon with all numerical schemes we implemented. This thickening of cracks was also noticed in other works in the damage community (Bleyer et al. [START_REF] Bleyer | Dynamic crack propagation with a variational phase-field model: limiting speed, crack branching and velocity-toughening mechanisms[END_REF], Li [START_REF] Li | Gradient Damage Modeling of Dynamic Brittle Fracture[END_REF]).

What we found out is that, in fact, α is just arbitrarily close to 1, as we are not allowed to have this value in an interval. In fact, by writing our energy as

E = Ω 1 2 E(α)|u | 2 + w 1 α + 1 2 w 1 2 |α | 2 dΩ (2.150)
and supposing that α=1 in an interval I, then for all β<0 sufficiently small and every point in the interval, we have

1 2 E (α)β + w 1 β -w 1 2 α β ≥ 0. (2.151)
We also have that α =0 and if we consider E(α)=(1 -α) 2 , then E (α)=0. As a result, we have w 1 β≥0 for all β≤0. This is a contradiction. We have thus proved that α cannot be equal to 1 in an interval. This absence of rigidity leads to another undesirable effect: a displacement field that is very irregular in the fractured regions (see Figure 2.16). This holds true even for the variational scheme, despite its tendency to regularize the numeric solution.

Direction of Cracks

When studying 2D and 3D ductile fracture, we can observe that the direction of cracks is not the same as in brittle fracture. In quasi-static, it has been shown in Tanne [START_REF] Tanne | Variational phase-field models from brittle to ductile fracture: nucleation and propagation[END_REF] that the direction of cracks depends on the yield stress σ Y , σ c := w 1 /E 0 and on the thickness of the specimen. When σ c >>σ Y , if we consider a rectangle under traction, the cracks that appear are inclined and close to one of the borders.

As in quasi-static, we obtain for dynamic loading that the angle also depends on the internal length and on height of the plaque.

One example of crack is shown in Figure 2.17. For this simulation, we considered E 0 =1, σ 0 Y =0.5 and =0.05. As for the boundary conditions, we considered α=0 at the extremities x=0 and x=1. At the instant t, the imposed displacements are u| x=0 = (0, 0) and u| x=1 = (t, 0). Figure 2.16: Displacement for a brittle material after rupture. We see that u is irregular in the region where the crack appeared.

Discontinuous Galerkin

So far, the discretization of all the variables was based on the classic Lagrange elements. We present now the Discontinuous Galerkin (DG) methods.

Discontinuous Galerkin methods are finite element methods that allow discontinuities in the discrete trial and test spaces. We are never able to create a real fracture if the displacement u is continuous: one (or more) element will suffer all the deformation. In a fragile quasi-static fracture, this caused errors in the profile of σ. From a practical point of view, it would be natural to allow discontinuities in the discrete space, making the DG approach very appealing.

Another reason is the study of a dynamic bar. We want to be able to continue the simulation after the first crack appeared and, when a crack appears, we want each part of the bar to act as an independent bar. This is only possible if the nodes of different elements are completely independent.

The objective of this section is to first present a heuristic derivation to the Poisson problem using the Symmetric Interior Penalty method (Arnold [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF]), as in Antonio Di Pietro and Ern [START_REF] Antonio | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF], and, from that, obtain a suitable form of energy that can be applied to gradient damage problems, both in quasi-static and dynamics.

There is, however, a price to pay: the number of nodes increases and the weak formulation is more complicated. In fact, the discontinuities appear on the variational formulation as integrals of jumps in the boundary between elements. Because of these new terms, the coercivity is no longer evident.

-Discontinuous Galerkin

One last important difference in terms of implementation is the application of boundary conditions: when using Lagrange elements, we imposed a displacement at the boundary. When using DG-methods, the boundary conditions are a part of the variational formulation.

Heuristic derivation

We consider the Poisson problem in a domain Ω:

-∆u = f in Ω, u = u 0 on ∂Ω, (2.152)
where f ∈ L 2 (Ω)the source term and the solution at the boundary is given by the function u 0 ∈ L 2 (Ω).

Chapter 2 -Numerical Implementation and Validation

We consider the Sobolev space H k (Ω) and define the broken Sobolev space on the mesh T :

H k broken (T ) := {v ∈ L 2 (Ω) : v| T ∈ H k (T )
, for all element T of the mesh T }.

(2.153) We search an approximate solution to (2.152) in the space H 1 broken (T ). Before we continue, we need some definitions: Definition 2.6.1. If F is a face in the mesh T , we say that (i) F is an interface if there are two different elements T 1 and T 2 in T such that F = ∂T 1 ∩ ∂T 2 ;

We denote the set of all interfaces by F i .

(ii) F is a boundary face if there are is one element T ∈ T such that F = ∂T ∩ ∂Ω;

We denote the set of all boundary faces (exterior faces) by F e .

Definition 2.6.2. Consider a function v ∈ H 1 broken (T ), two mesh elements T + and T -sharing a face, and the normal n from T -to T + . We define:

(i) Average of a function:

v = 1 2 v| T + + v| T -. (2.154) (ii) Jump of a function: v = v| T + -v| T -. (2.155)
It's easy to verify the jump identity:

uv = u v + u v . ( 2

.156)

We have now all the necessary definitions and we can proceed to the calculations.

Suppose that u ∈ H 1 broken (T ) ∩ H 2 (Ω) is the exact solution to (2.152). Then, for all v ∈ H 1 broken , we have

Ω -∆u • v = T ∈T T -∆u • v = Ω f • v. (2.157)
If we integrate by parts, we have where n T is the outward normal to T on ∂T . Each face F in F i is shared by two mesh elements T + and T -. Hence,

T ∈T T ∇u • ∇v - T ∈T ∂T ∇u • n T • v = Ω f • v, (2.158) 2.6 -Discontinuous Galerkin n T + T - v v v| T - v| T +
∂T + ∩F ∇u • n T + • v + ∂T -∩F ∇u • n T -• v = F ∇u • n T + • v - F ∇u • n T + • v = F ∇u • n • v , (2.159)
where

n = n T + = -n T -. Therefore, T ∈T T ∇u•∇v- F i ∈F i F i ∇u•n •v - Fe∈Fe Fe ∇u•n •v = Ω f •v. (2.160)
The solution u is continuous and its gradient is continuous. This means that u = 0 and ∇u = 0 in every face.

We then have

T ∈T T ∇u•∇v - F i ∈F i F i ∇u•n v - Fe∈Fe Fe ∇u•n •v = Ω f •v. (2.161)
For (w, v) ∈ H 1 broken (T ) × H 1 broken (T ), we define the forms

a 0 (w, v) = Ω ∇w • ∇v - F i ∈F i F i ∇w • n v - Fe∈Fe Fe ∇w • n • v, (2.162) and L 0 (v) = Ω f • v. (2.163)
By our calculations, the form a 0 is consistent. We also want the discrete bilinear form to preserve the original symmetry of the exact bilinear form. We define

a 1 (w, v) = Ω ∇w • ∇v+ - F i ∈F i F i ∇w • n v - F i ∈F i F i w ∇v • n + - Fe∈Fe Fe ∇w • n • v - Fe∈Fe Fe w • ∇v • n , (2.164)
and

L 1 (v) = Ω f • v - Fe∈Fe Fe ∇v • n • u 0 . (2.165)
The form a 1 is consistent, symmetric and a 1 (u, v) = L 1 (v), for all v ∈ H 1 broken (T ). We now fix two constants C 1 and C 2 . We define

a(w, v) = Ω ∇w • ∇v+ - F i ∈F i F i ∇w • n v - F i ∈F i F i w ∇v • n + C 1 h f i ∈F i F i w v + - Fe∈Fe Fe ∇w•n •v - Fe∈Fe Fe w• ∇v•n + C 2 h Fe∈Fe Fe w•v, (2.166)
and

L(v) = Ω f •v - Fe∈Fe Fe u 0 • ∇v•n + C 2 h Fe∈Fe Fe u 0 •v, (2.167)
where h represents the diameter of the smaller element of the mesh.

-Discontinuous Galerkin

If C 1 and C 2 are large enough, the bilinear form a is coercive (Antonio Di Pietro and Ern [START_REF] Antonio | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]). Using the Lax-Milgram Lemma, we know there exists one unique u ∈ H 1 broken (T ) such that a(u, v) = L(v), for all v ∈ H 1 broken (T ).

(2.168)

As a consequence, there exists a unique u ∈ H 1 broken (T ) that minimizes the functional

E(u) = 1 2 Ω ||∇u|| 2 + - F i ∈F i F i ∇u•n u + C 1 h F i ∈F i F i u 2 + - Fe∈Fe Fe ∇u•n •u + C 2 2h
fe∈Fe Fe

||u|| 2 + - Ω f •u + Fe∈Fe Fe u 0 • ∇u•n - C 2 h Fe∈Fe Fe u 0 •u.
(2.169)

Application to the gradient damage model

We now apply the previous results to the gradient damage model. Using the DG elements, the discrete problem allows discontinuities in u.

We write the quasi-static displacement problem as

-div E(α)(ε(u) -ε p ) = 0 in Ω, u = u 0 on ∂Ω.
(2.170)

We want to find u ∈ H 1 broken (T ) such that, for all v ∈ H 1 broken (T )

-Ω div E(α)ε(u) v = -Ω div E(α)ε p v in Ω, u = u 0 on ∂Ω.
(2.171)

We integrate by parts:

Ω div E(α)ε p v = T ∈T T div E(α)ε p v = - T ∈T T E(α)ε p ε(v)+ T ∈T ∂T E(α)ε p •n v (2.172)
From the calculations in the previous section, the bilinear form can be defined as

a(w, v) = Ω E(α)ε(w)•ε(v)+ - F i ∈F i F i E(α)ε(w)•n v - F i ∈F i F i E(α)w ε(v)•n + C 1 h f i ∈F i F i E(α)w v + - Fe∈Fe Fe E(α)ε(w)•n •v - Fe∈Fe Fe E(α)w• ε(v)•n + C 2 h Fe∈Fe Fe E(α)w•v.
(2.173)

The form L is now

L(v) = - Fe∈Fe Fe u 0 • ε(v)•n + C 2 h Fe∈Fe Fe u 0 •v - Ω E(α)ε p ε(v)+ T ∈T ∂T E(α)ε p •n v . ( 2

.174)

Thus we can find the displacement u by minimizing the functional

E DG (u) = 1 2 Ω E(α)(ε(u) -ε p )•(ε(u) -ε p )+ - F i ∈F i F i ε(u)•n E(α)u + C 1 h F i ∈F i F i E(α)u u + - Fe∈Fe Fe E(α) ε(u)•n •u + C 2 2h
Fe∈Fe Fe

E(α)||u|| 2 + Fe∈Fe Fe E(α)u 0 • ε(u)•n - C 2 h Fe∈Fe Fe E(α)u 0 •u+ F i ∈F i F i E(α)ε p u + Fe∈Fe Fe E(α)u 0 • ε p •n .
(2.175)

To solve the damage problem, we have to minimize this same energy with respect to α. By doing so, we are able to allow discontinuities in the displacement when α=1 and, reciprocally, enforce α=1 where there is a displacement in u.

It is important to emphasize that using the standard one-dimensional energy

E = Ω 1 2 E(α)|ε(u) -ε p | 2 + w 1 α + 1 2 w 1 2 |α | 2
for the damage problems yields invalid results. In fact, for a broken bar, we have a discontinuous 2.6 -Discontinuous Galerkin displacement and zero strain everywhere, leading to the absence of damage in the bar. This is not a problem with Lagrange elements because jumps are not allowed and all the deformation is located in one element. By refining the mesh, the value of u tends to +∞. The term 1 2 E(α)|u | 2 is large if E(α) = 0.

Quasi-static results

We now apply this same formulation for the quasi-static brittle damage evolution problem. We consider E(α)=α and w(α)=α, L=1 and =0.1.

In order to illustrate the advantages of the DG formulation, we show the results for a one-dimensional coarse mesh with only 20 elements.

We can see in Figure 2.19 the displacements obtained for each case. We can see that, for the standard formulation, since no discontinuities in the displacement are allowed, most of the strain is localized in the center of the bar, and the residual strain in the rest of the bar is much larger than the one obtained when using DG elements. When using the DG formulation, the residual strain is irrelevant and there is a clear discontinuity in the middle, as expected.

-Discontinuous Galerkin

For a quantitative result, we first compare the stress in the bar after fracture. The expect value is clearly zero. After numerically differentiating u and multiplying this value by E(α), we obtain σ=0.015 for the standard formulation and σ=0.0005 when using the DG formulation.

We then compare the energy dissipated in the damage process. As we already know, the dissipated energy for this problem is given by D = 4 √ 2 /3 = 0.189. We obtain D=0.242 for the standard formulation and D=0.179 for the DG formulation, an error five times smaller.

In summary, we can obtain a real discontinuity, a better damage profile and a lower residual stress when using DG methods.

Dynamic damage problem (Newmark scheme)

The natural progression of this work is the implementation of the Discontinuous Garlerkin method for solving the wave equation and the dynamics of the damage problem. We will follow the same approach as Grote et al. [START_REF] Grote | Discontinuous galerkin finite element method for the wave equation[END_REF].

The main calculations have already been performed for the quasi-static problem. In that case, the equation of displacement could be obtained by minimizing the energy given by equation (2.175). We obtained a problem defined by

a h (u, v) = L h (v), ∀v ∈ H 1 broken (2.176)
where a h (•, •) is continuous, symmetric and coercive, as long as α≤α max <1.

The exact expression for each form is written in equations (2.173) and (2.174). This problem is equivalent to solving divσ = 0 on the domain Ω. We recall that the boundary condition u 0 is already included in the energy (2.175).

The dynamic equation consists in solving ρü + divσ = 0 in Ω u = u 0 (t) on ∂Ω.

(2.177)

for given initial conditions u(t=0) and u(t=0). From these considerations, the displacement problem can be expressed as finding u such that

Ω ü•v + a h (u, v) = L h (v), ∀v ∈ H 1 broken . (2.178)
The integration in time will be performed using the explicit Newmark scheme. At the time instant t i , we calculate the acceleration Ü i+1 from the displacement U i , damage α and plastic strain ε p (these last two variables are present in the variational formulation of the problem and will be omitted here):

Ω Ü i+1 •v + a h (U i , v) = L h (v), ∀v ∈ H 1 broken . (2.179)
With this change in the calculation of the acceleration, the numerical solution is the same as in section 2.3.1.

With this approach, we do not gain precision and the numerical calculations are longer. The main advantage of this approach is that if α=1 in a certain point, a discontinuity in the displacement becomes possible. We consider a bar given by the interval [-L, L], where ε=1 and u=0 at t = 0, and we block the extremities u(±L)=±L. It is clear the this system does not move. We now impose a thin damage profile, with α=1 in the center. As a result, each side of the bar behaves as an independent bar, where one extremity is blocked and the other one is free. This does not happen when using the standard formulation. The comparison of damage profiles is show in This approach, unfortunately, does not solve the problem of thickening of cracks or irregular displacement near the crack tips in the general damage Chapter 3

Application: Explosion and Fragmentation of a Ring

The objective of this chapter is to model a cylinder in expansion, until it breaks in many parts. We then want to count the number of fragments obtained and determine how it depends on the parameters used. The main difficulty is that the cylinder is invariant by rotation and, therefore, we should obtain a damage configuration that is also axially symmetric. Surprisingly, this is not what happens, as we obtain radial fissures somewhat evenly spaced, as shown in Figure 3

Chapter 3 -Application: Explosion and Fragmentation of a Ring

We want to know what is causing damage profiles to localize in this axial symmetric problem. Our main focus is on the study of the role of each material parameter in the number of fragments in a one-dimensional periodic bar. The first step consists of establishing the dimensionless problem and knowing how many effective variables are in the problem. We then run a series of numerical simulations in order to obtain how many cracks appear in the bar for each set of parameter. We noticed that, in a homogeneous bar, small perturbations (even purely numerical perturbations) were sufficient to create localized damage profiles.

We then study the linearised system and show that some modes of perturbation are amplified faster than the others. Our main result is that, by finding these modes, we can predict the number of cracks that appear in the simulations. From that, we can obtain an analytic approximation to the number of the cracks in the bar, depending on the material and loading parameters.

We conclude by showing the simulations for a complete three dimensional ring to validate the results using the one dimensional simplification.

Fragmentation of a Brittle Ring

1D Periodic Bar

We are interested in obtaining the number of fragments of a ring under expansion. We consider first a brittle material, and the plasticity will be added in the next section. Instead of working with a ring in a 3D 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, (3.1)
for every t ∈ R.

We also suppose that we start our study with a completely sound 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]. (3.2)
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 expansion. At an instant t, the perimeter of the ring is L+ ε0 Lt. The displacement u at the 3.1 -Fragmentation of a Brittle Ring extremities and ε must satisfy

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

At each time instant t, we write

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

(3.5)

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. (3.6)
By differentiating the above equation with respect to x, we find that (u * ) = ε * . It is clear that that u * is periodic.

The stress can be written as σ(x, t) = E(α)((u * ) + ε0 t).

(3.7)

The dynamic equation is

ρü * (x, t) = ρü(x, t) = σ . (3.8)
Finally, we write the energy to be minimized in the evolution of damage as

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

Dimensionless Parameters

The objective of this section is to show that the problem of a bar made of a brittle material under traction depends only on two dimensionless parameters.

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

The study of the damage evolution consists in defining an energy and finding its minimum with respect to α:

E(u, α) = Ω 1 2 E 0 a(α) ε(u) 2 + w 1 α + 1 2 w 1 2 (α ) 2 dΩ. (3.10)
The dynamic equation in 1-D can be written as

ρü = σ = E 0 a(α)ε(u) (3.11)
and we impose a displacement on the extremities:

u(x=0, t) = 0 u(x=L, t) = ε0 Lt. (3.12)
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Ω. (3.13)
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=x/L and t=T t, for some constant T >0 that we will specify later.

If

-Fragmentation of a Brittle Ring

Suppose we have a constant U 0 >0. We define ũ(x, t) := 1 U 0 u(x, t) α(x, t) := α(x, t).

(3.17)

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Ω (3.18) and 1 T 2 d 2 ũ d t2 = E 0 ρL 2 d dx a(α) dũ dx . ( 3.19)
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Ω (3.20)
and the dynamics of the system is

d 2 ũ d t2 = d dx a(α) dũ dx . (3.21)
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 damage profile α minimizes the energy E taking into account the irreversibility condition, where E is given by Example 3.1.1. Suppose we want to study the fracture of a bar made of a brittle material. We suppose this bar is 10 cm long and is being stretched with a constant speed of 100m/s. For this material, we suppose that 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. We also consider E(α)=E 0 (1-α) 2 and w(α)=w 1 α.

E(u, α) = Ω 1 2 a(α) ε(u) 2 + α + 1 2 2 α 2 dΩ. ( 3
When using this damage gradient model in a quasi-static scenario, we know from Pham et al. [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF] that

w 1 = σ 2 C E 0 = (1000 • 10 6 ) 2 210 • 10 9 Pa = 4.7MPa. (3.25)
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. (3.26)
Again from [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF], we know that the energy dissipated by the fracturing process can be written as The deformation can be now be found:

G c = σ C 4 √ 2 3 . ( 3
ε = du dx = U 0 L dũ dx = 4.73 • 10 -3 dũ dx . (3.30)
The dimensionless deformation speed is

ε0 = ε0 T L U 0 = 4.3 • 10 -3 ε0 = 0.43. ( 3

.31)

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.

-Fragmentation of a Brittle Ring

To conclude, this section shows that we only need to parameters in our simulations: the internal length and the strain rate ε. For this reason, we set all the other parameters equal to one. The above example shows the values of and ε that we can obtain for a realistic material. Even though we are not interested in the behaviour of one specific material, the values we calculated in the example allow us to have an idea of the expected range of each parameter.

Influence of Each Parameter

The following step is understanding the influence of the parameters in the number of cracks in the bar. This problem is highly sensitive to numerical parameters, specially 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 set of parameters we want to study, we run several simulations with different values of ∆x and see if the 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 number of fragments seems to converge as ∆x → 0.

We illustrate the damage evolution in Figure 3.2, where the damage profiles are shown for different time instants. The system evolves uniformly until it reaches a critical value, causing stress and damage to localize in multiple places at the same time. Each point where α=1 represents a crack. We recall that the dimensionless problem has only two parameters ( ε and ).

For this bar of unitary length, the mesh size h is the inverse of the number of elements. In Figure 3.3, we can see 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 h ≤ /3. For other values of ε0 , this relation also seems to be valid.

In Figure 3.3, we can see 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 3.3, we can see that the number of fragments increases linearly with 1/ for ε0 = 0.5. This linear behavior holds true for every other We can see the influence of ε0 in Figure 3.4. We make ε0 vary between 10 -4 and 10 2 . The number of cracks change from less than 500 to over 1000.

We conclude with two remarks. Firstly, we can see that the behaviour 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 number of cracks in the end. In fact, we launched these simulations using an initial perturbation on the displacement or the damage profile at t=0. We tried different shapes but we emphasize that the number of cracks in the end is the same if the amplitude of the perturbation is sufficiently small. There is, however, some minor changes in how they appear and how many of these cracks develop until total failure.

Stability of the Homogeneous Solution (Brittle Material)

In this, section, we study the evolution of a small perturbation to the system in order to explain the number of cracks that appear in the bar. This

Study of a Perturbation

We consider a periodic bar [0, L] made of a brittle material such that E(α)=E 0 (1 -α) 2 and w(α)=w 1 α.

We recall that the homogeneous strain and damage can be written as .32) We are going to study the influence of a perturbation of the form (δε, δα). The dynamics of the system can be written as

ε 0 = εt and α 0 = 1 - w 1 E 0 ε 2 0 . ( 3
ρε = (1 -α) 2 E 0 ε . (3.33)
If we take (ε, α)=(ε 0 + δε, α 0 + δα) and consider only the terms of the first order, we obtain

ρδ ε = (1 -α 0 ) 2 E 0 δε -2E 0 (1 -α 0 )ε 0 δα (3.34)
or, equivalently, We now write the damage criterion and look at the terms of the first order:

δα = 1 2 w 1 E 0 ε 3 0 δε - ρε 0 2w 1 δ ε. ( 3
E 0 ε 2 0 δα + 2(α 0 -1)E 0 ε 0 δε -w 1 2 δα = 0. (3.36)
We now use the expression of δα found in the dynamic evolution equation and we find

E 0 ε 2 0 δα - 2w 1 ε 0 δε -w 1 2 1 2 w 1 E 0 ε 3 0 δε - ρε 0 2w 1 δ ε = 0. ( 3.37)
We now differentiate this expression twice in space

δα = 2w 1 E 0 ε 3 0 δε + w 2 1 2 2E 2 0 ε 5 0 δε - ρ 2 2E 0 ε 0 δ ε . (3.38)
and we replace the expression of δα in equation (3.34) to find

ρδ ε - ρw 1 2 E 0 ε 2 0 δ ε = - 3w 2 1 E 0 ε 4 0 δε - w 3 1 2 E 2 0 ε 6 0 δε . (3.39)
We now use that ε 0 = εt, to find the equation dictating the evolution of the perturbation:

ρδ ε - ρw 1 2 E 0 ε2 t 2 δ ε = - 3w 2 1 E 0 ε4 t 4 δε - w 3 1 2 E 2 0 ε6 t 6 δε (3.40)
We try to find a solution in the form δε= sin(kx)f (t). We thus obtain the ordinary differential equation to be solved: With that in mind, we propose a criterion to determine the number of cracks in the system by studying the amplitude of the each perturbation. We are going to consider that the first perturbation to reach the threshold S will be the one that dominates the system and dictates the number of cracks that are going to appear.

ρ 1 + w 1 2 k 2 E 0 ε2 t 2 f (t) = 3w 2 1 k 2 E 0 ε4 t 4 - w 3 1 2 k 4 E 2 0 ε6 t 6 f (t) (3.41)
More precisely, for each perturbation N , we denote by t N the time such that the amplitude reaches the threshold, that is, |f N (t N )|=S (we are going to ignore the modes that do not reach the threshold).

The mode N * is the first one to reach the threshold (t N * ≤t n , ∀n) and, therefore, N * cracks will appear. Table 3.1: Comparison between number of cracks found in the simulations and using the numeric calculations.

As we can see in Table 3.1, our method gives us a good approximation for the number of cracks. We have, however, to pay attention to some details: for instance, if we take =10 -3 , we find t 100 =1.0323 and t 1600 =1.0311. If we take =10 -4 , we find t 1000 =1.0032 and t 160 =1.0028.

-Stability of the Homogeneous Solution (Brittle Material)

Figure 3.8: Amplitude of f N (t) as a function of N and t=5 fixed. We consider =0.02.

If =10 -4 and ε=10 -5 , the values of t 100 and t 200 are very close (t 100 =100000.0357 and t 200 =100000.0328). This means that we cannot tell if the bar is more likely to have 100 or 200 cracks.

Another important aspect is that this method does not capture well the influence of the strain rate. If, for instance, we fix =10 -3 and change ε, we don't see a change in result when ε is very small. Even if the number of cracks, in some cases, isn't exactly the same found in the simulations, we find the correct instant for when the perturbations become visible.

If we know the time t N * when the threshold is reached for the mode N * , then we can calculate the damage value:

α * = 1 - w 1 E 0 ε2 t 2 N * . (3.42)
We compare the values of damage for which the perturbations become visible in Figure 3.15. We can see that we have very precise results (and it works well for every internal length and strain rate).

In order to validate our reasoning and calculations, we are going to make a test. Instead of supposing that all modes have the same initial amplitude, we are going to force one mode to have a larger amplitude. For modes that develop at roughly the same time, this should allow us to control the number of cracks that are going to appear. Figure 3.9: Value of α when the perturbations become visible ( =0.01). In blue, the values found in the simulations. In red, the results using the method described above.

We recall that:

• purely numerical perturbations are of the order of 10 -16 ;

• we are going to impose perturbations with amplitude between 10 -12 and 10 -3 .

• a perturbation sin(2πN x/L) creates N cracks (if N is one of the modes that increases fast enough);

• we can also find the value of α when cracks appear (the number of cracks is already known).

Since we only have one mode acting, we are going to impose a small initial perturbation A 0 sin(2πN x/L) and see when its amplitude reaches one.

In the table below, =0.02 and ε=1. Our method could not accurately predict how many cracks we should have since the instants when modes between 8 and 12 reach the threshold were almost the same. If we impose a perturbation sin(2πN x), with N between 8 et 12, we should find N cracks and the value of α when instabilities become visible should depend only on the initial amplitude.

This is exactly what we found with this simple test. We compared these values of α found using our method, and imposing a initial perturbation to ε and running the simulation of the whole system. The values can be seen in Table 3.2.

-Stability of the Homogeneous Solution (Brittle Material)

initial amplitude α in our method α in the complete simulation 10 -8 0.537 0.564 10 -7 0.392 0.419 10 -5 0.390 0.419 10 -5 0.390 0.419 10 -4 0.331 0.365 Table 3.2: Comparison between the value of damage for which instabilities become visible in the simulations and using the numeric calculations when a perturbation is imposed.

Analytic Approximation

The method developed in the previous section gives us good results when compared to the full simulation. In order to go further, we would like to be able to tell explicitly how each parameter influences the number of cracks and when they start to localize.

We recall that each perturbation is written as sin(kx)f (t), where

ρ 1 + w 1 2 k 2 E 0 ε2 t 2 f (t) = 3w 2 1 k 2 E 0 ε4 t 4 - w 3 1 2 k 4 E 2 0 ε6 t 6 f (t).
(3.43)

We now define the dimensionless ratios

τ = E 0 w 1 ε k t and λ = w 1 ρ 1 ε . (3.44)
We define y(τ )=f (t) and we have

y (τ ) = λ 2 g(τ )y(τ ) (3.45)
where

g(τ ) = 3τ 2 -1 τ 4 (τ 2 + 1) . ( 3

.46)

We are going to study the behaviour of g(τ ) in three regions:

τ << 1/ √ 3
In this region, g(τ )<0. As a result, the growth of y is not enough to influence the system. By using Airy's function, we can find a solution to this problem and see that it does not grow fast enough.

τ >> 1/ √ 3
We can find a numerical solution to this equation and see that the function y(τ ) increases really fast.

We consider the change of variable y = exp(λz). We have y = λ exp(z)z , y = (λz + λ 2 z 2 ) exp(z) and

1 λ z + (z ) 2 = g(τ ). (3.49)
Since we are working with large values of λ, we consider the Taylor expansion z = z 0 + z 1 /λ. By taking the terms that depend on λ 0 , we find z 2 0 = g(τ ). Therefore z 0 = ± g(τ ). Since we are interested by the functions that grow very fast, we are going to take z 0 =+ g(τ ).

Chapter 3 -Application: Explosion and Fragmentation of a Ring

The terms preceded by λ -1 , are z 0 + 2z 0 z 1 = 0. Therefore

z 1 = -z 0 /2z 0 = - 1 2 d dτ log z 0 = - 1 4 d dτ log g(τ ). (3.50)
If we want y(τ =τ 0 )=1, we must have z(τ =τ 0 )=0. We define

G(τ ) = τ 1 √ 3 g(s)ds (3.51) Thus z = G(τ ) -G(τ 0 ) -1 4λ (log g(τ ) -log g(τ 0
)). We find the approximation for y:

y(τ ) ≈ exp(λz) = exp(λ(G(τ ) -G(τ 0 ))) g(τ 0 ) 1/4 g(τ ) 1/4 (3.52)
With this approximation, if λ is large enough, we have a very good approximation and we cannot differentiate the exact from the approximate curve.

The adimensional time τ 0 is the correspondent to t c = w 1 /E 0 / ε. Damage and the perturbation start to evolve at τ 0 . We fix a threshold S and, for a perturbation that has unitary amplitude, we can find the adimensional time τ * for which y(τ * )=S.

By taking the logarithm on both sides of the equation, τ * satisfies

G(τ * ) -G(τ 0 ) - 1 4λ (log g(τ * ) -log g(τ 0 )) = log S λ . ( 3.53)
Our goal is to find the first perturbation to reach the threshold. The real instant where that happens will be denoted by t * . We remark that

t * = τ * t c τ 0 . ( 3.54)
Our problem consist in finding the value of k that minimises t * :

min k>0 t * . ( 3.55)
We recall that

τ = E 0 w 1 ε k t. (3.56)
The values of τ and t are connected by the material constants and k. By changing the value of k, we also change τ 0 and, if we know the value τ 0 , we can find k. For this reason, we are going to solve the problem min

τ 0 t * , (3.57)

-Stability of the Homogeneous Solution (Brittle Material)

which is equivalent to finding min

τ 0 τ * τ 0 . ( 3.58)
We define the auxiliary function

ψ(τ ) = G(τ ) - log g(τ ) 4λ . (3.59)
It is easy to see that

ψ(τ * ) = ψ(τ 0 ) + log S λ . ( 3.60)
We take the Taylor expansion of ψ to find

ψ(τ * ) ≈ ψ(τ 0 ) + ψ (τ 0 )(τ * -τ 0 ). (3.61)
These two equations allow us to find an approximation to τ * by using

ψ (τ 0 )(τ * -τ 0 ) = log S λ . (3.62)
Therefore, we have

τ * ≈ τ 0 + log S λψ (τ 0 ) . ( 3.63)
The minimisation problem (3.58) becomes min

τ 0 1 + log S λτ 0 ψ (τ 0 ) , (3.64)
which is equivalent to maximising the function τ ψ (τ ). Since we are working with large values of λ and the value of g(τ ) has only some small variations close to 1, we are going to consider that the function τ ψ (τ ) does not depend on λ.

We only have to calculate once the maximum value of τ ψ (τ ). This value is reached for τ 0 ≈1.0, and τ 0 ψ (τ 0 )≈1.05.

Since we now know τ 0 , we can find the number N of cracks :

N = kL 2π = L 2π E 0 w 1 εt c τ 0 . ( 3.65)
We simplify this expression to find We use equations (3.54) and (3.63) to find

N ≈ L 2π (3.66)
t * = t c 1 + log S λτ 0 ψ (τ 0 ) ≈ t c 1 + log S 1.05λ . ( 3.67)
The value of damage for which the localized solutions appear is

α * = 1 - w 1 E 0 ( εt * ) 2 ≈ 1 - 1 (1 + log S 1.05λ ) 2 = . (3.68)
By replacing λ, we find

α * ≈ 1 - 1 (1 + ρ ε 1.05 √ w 1 log S) 2 (3.69)

Results

We can see the results in Figures 3.14 and 3.15. To summarize:

• we have good approximations for the number of cracks and damage value;

• the number of cracks depends exclusively on , the small influence of the strain rate is not taken into account;

• the chosen threshold influences the instant where the perturbations are considered visible, but not the number of cracks;

• the approximation for α min is less precise for large strain rates.

Complete 3-D Simulation of a Ring

The last part of this work consists of comparing the results obtained for the 1-D periodic bar to those obtained for a 2-D or 3-D ring. More precisely, we want to know if the number of cracks we obtained in 1-D or 3-D are the same for a brittle material.

The multidimensional problem is clearly more complex. In the 1-D bar, we had two main parameters: the internal length and the strain rate. For a complete 3-D ring, we also have to study the influence of the thickness, the amplitude and the form of the applied loading, and there are more available options for the shape of the initial perturbation.

After running a set of simulations, we obtained answers for some of these questions: the thickness has little influence on the number of cracks, as long as it is small compared to the internal radius, but not too small as a very thinner rings pose problems for damage localization; the 2-D and the 3-D ring behave in the same way, so most of the tests were done in a 2-D scenario.

As for the loading, we considered an imposed displacement at the internal surface or an internal pressure. The form of this loading (that is, uniform or almost uniform with sufficiently small perturbations) has little influence on the number of cracks; it does, however, play an important role in facilitating damage localization and the instant when these damage profiles appears. We place ourselves in the case of small strain rates (in the dimensionless problem, when compared to one). For larger values of the strain rate (close or larger than one), the localization of damage is not guaranteed or is very close to α=1, making it difficult to distinguish the cracks.

When running our simulations, we used a structured mesh, as nodes are aligned with the direction of the movement and allow us to better represent the dynamics of the system. Different mesh sizes were used, giving us close (with a difference of less than 20%) numbers of cracks. The mesh is illustrated in Figure 3.17.

To be precise, the system starts at rest u(t=0)=v(t=0)=0. The external boundary is free of constraints (σ rr =0) and a displacement proportional to the time is imposed in the inner boundary (u r = ε0 t).

The radius of the internal hole is R i =1 and we tried different values of the external radius R e between 1.05 and 1.2, yielding similar results.

For the 1-D problem, we showed that the number of cracks was very close to L/(2π ) and that the strain rate has only a small influence in the number of damage profiles that appear. In a multidimensional ring, the length L of the bar represents the perimeter, so the expected number of cracks is 1/ .

For the simulations, we considered the dimensionless problem by setting E=1, ν=0.3, ρ=1, w 1 =1, R e =1.1 and placed ourselves in the case of plane stress.

Table 3.3 shows the results we obtained for a fixed strain rate ε=0.5 and different values of internal length. We expected the number of cracks to be proportional to 1/ , but that is not what happens here. For lower values of , we overestimated the number of cracks and, for higher values of , we underestimate the number of cracks. We used a mesh of size 10 We see in table 3.4 that the strain rate plays an important role in the number of cracks. We fixed =0.063 and we see a significant change in the number of cracks for ε=0.05 and ε=0.5. For values lower than ε=0.05, we didn't see a significant change in the number of cracks. For values higher than ε=0.5, the localization of cracks was less evident and the results were not considered. This means that the 2-D and 3-D cases are more complex and that the study of the one-dimensional bar is not enough to accurately predict the number of cracks. We have, however, the same tendency (increase for smaller internal lengths and higher strain rates) and the the number of cracks obtained is on the same order of the predicted value.

-Stability of the Homogeneous Solution (Ductile Material)

Other sources of difference might be the non-uniformity of the damage in the radial coordinate and the presence of more components in the stress tensor.

Stability of the Homogeneous Solution (Ductile Material)

We will now study the behaviour of the system when we also consider plastic strain. We consider the dimensionless model described by the energy

L 0 1 2 (1 -α) 2 (ε -p) 2 + α + 1 2 2 (α ) 2 + (1 -α) 2 σ p p dx, (3.70)
where the new parameters σ p is defined by the ratio of the original yield stress and the critical damage stress. Since the complete plasticity problem is difficult to study, we will consider p=p and search for p that minimises the energy.

We have the following dimensionless equations:

     dynamic : ε = (1 -α) 2 (ε -p) plasticity : ε -p = σ p damage : -(1 -α)(ε -p) 2 + 1 -2 α -2σ p (1 -α)p = 0. (3.71)
We first consider the homogeneous solution. It is easy to see that ε 0 = εt and p 0 = εt -σ p .

(3.72)

Damage begins to evolve at t=(1 + σ 2 p )/(2σ p ε) and can be written as

α 0 (t) = 1 - 1 (ε 0 -p 0 ) 2 + 2σ p p 0 = 1 - 1 2σ p εt -σ 2 p . ( 3.73)
We consider the perturbation (δε, δp, δα). We have

     δ ε = -2σ p (1 -α 0 )δα δε = δp (σ 2 p + 2σ p p 0 )δα -2 δα -2σ p (1 -α 0 )δp = 0. (3.74)
As in the brittle case, we replace δα in the last equation, differentiate it twice in space and find

(σ 2 p + 2σ p p 0 )δα + 2 2σ p (1 -α 0 ) δ ε -2σ p (1 -α 0 )δε = 0. (3.75)

Chapter 3 -Application: Explosion and Fragmentation of a Ring

We now replace δα in the movement equation:

δ ε = - 4σ 2 p (1 -α 0 ) 2 σ 2 p + 2σ p p 0 δε + 2 σ 2 p + 2σ p p 0 δ ε . (3.76)
We now inject a perturbation δε=sin(kx)f (t). We have

(1 + 2 k 2 σ 2 p + 2σ p p 0 )f (t) = 4σ 2 p (1 -α 0 ) 2 k 2 σ 2 p + 2σ p p 0 f (t). (3.77)
This gives us

f (t) = g k (t)f (t) with g k (t) = 4σ 2 p (1 -α 0 ) 2 k 2 σ 2 p + 2σ p p 0 + 2 k 2 .
(3.78)

We now replace p 0 and α 0 to find

f (t) = g k (t)f (t) with g k (t) = 4 (2 εt -σ p ) 2 2σp εt-σ 2 p k 2 + 2 (3.79)
The function g k (t) is always positive and decreases for sufficiently large times. Besides, the bigger the value of k, the slower this function decreases. For this reason, the value of f k (t) increases with k, as we can see in Figures 3. [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] We conclude that the number of cracks found will depend on the mesh used since finer meshes allow more modes to appear. This is in agreement with the simulations, both in 1-D and in 2-D. When refining the mesh, we the plastic strain is concentrated in very few elements and the simulations do not seem to converge when refining the mesh. This is illustrated in the twodimensional simulations of a ductile ring, where we impose the displacement of the interior surface. We see that if the ring is thick enough, the cracks are no longer in the radial direction, but inclined, as observed in rectangles by Tanne [START_REF] Tanne | Variational phase-field models from brittle to ductile fracture: nucleation and propagation[END_REF]. In the simulations, we use an external radius of 1.2 and an imposed displacement radial displacement u R =0.5•t. We use the same time step ∆t=10 3 and we vary the number of elements, always respecting the stability condition. We see that multiple cracks appear and they cross each other (Figure 3. [START_REF] Grote | Discontinuous galerkin finite element method for the wave equation[END_REF].

We see in Figure 3.22 that the cracks did not fully develop. To better analyse what happened, we look at damage profile at the external surface when angle varies between 0 and 2π. We show the results for 500 elements (Figure 3.23) and 2000 elements (Figure 3.24) in the perimeter. The number of elements in the radial direction has little influence in the number of cracks. We can see that the damage does not converge and we obtain more localized damage profiles when we refine the mesh.

Conclusion of the Chapter

We have studied the development of damage profiles in an one-dimensional periodic bar. The main objective was to understand the unexpected appearance of localized profiles in a homogeneous bar. We showed that these profiles appear because small perturbations (of amplitude equivalent to the numerical errors) were amplified to the point of becoming visible in the system. We then proceeded with an analytic study, where we found an approximation for the number of damage profiles that appeared. This analytic approximation confirms that the internal length is the main factor dictating the number of cracks. The last part was a comparison between the one-dimensional bar and a multidimensional ring. There were some differences in the results, but overall, the order of the number of cracks is preserved, as well as the tendency to obtain more cracks for smaller internal lengths or larger strain rates.

The main difference comes when comparing to energy models (Grady [START_REF] Grady | Local inertial effects in dynamic fragmentation[END_REF]), that predict that the number of fragments is proportional to ε2/3 . These models use the assumption that all kinetic energy is transformed in surface energy, and, therefore, increasing the number of cracks. In our model, this hypothesis no longer holds as an important part the of energy is dissipated as homogeneous damage.

Lastly, we see a difference of behaviour between the brittle and ductile cases. In the linearised system of equations obtained for a ductile material, the term corresponding to the perturbation of the elastic energy disappears. The equation that describes the amplitude of the perturbation changes drastically and there is no longer a value of the period k that maximises this amplitude. This means that, if the mesh is fine enough, the growth of the perturbation will be captured in the brittle case. In the ductile case, however, finer meshes allow us to capture more modes of perturbation. Since the growth of amplitude of these modes is not bounded, we do not obtain convergence.

Shockless Spalling

The identification of the parameters used in the gradient damage model is of extreme importance if such model is to be used in an industrial context. In this chapter, we detail how these parameters can be calibrated from a set of experiments.

The objective of this chapter is to study the shockless spalling test of a ceramic material and, from the results obtained in the experiments, we want to propose a model representing the material behaviour with as few parameters as possible, otherwise it would be impossible to identify them all. These tests also allow us to better understand how the strain rate influences the critical stress and the dissipated energy.

We remark that we are interested in the behavior of a brittle material since a ductile material would require more parameters to be calibrated and, therefore, observing the influence of the strain rate would be more difficult. Our model depends thus on only three parameters (instead of four, for the ductile case): the internal length , the function E(α) describing how damage affects de rigidity and the function w(α) describing the dissipation of energy.

The difficulties we faced come from the fact that we are working with shock-waves and the material in question does not have the same properties in quasi-static and in dynamics. We will begin this chapter by describing the experiments and showing that calibrating the parameters using only quasistatic tests is not enough to account for all the dynamic phenomena.

With that said, the emphasize that we want to propose material parameters that work well for both quasi-static and dynamic loadings.

Numerical Implementation

The cracks, the imposed pressure and the velocity at the free surface were uniform in the transverse direction, allowing us to model the specimen as a one-dimensional bar. Simulations were also performed in a 2D or 3D scenario, but there was no significant change in results, as we'll see in the next sections.

This approach has two main advantages:

1. working with a one-dimensional model allows us to perform many simulations, even if we have a very fine mesh;

2. we don't have to worry about different models of traction-compression asymetry, we will just consider the sign of the stress σ.

We model the specimen as bar [0, L]. At the extremity x=0, we impose a controlled compression wave. The elastic behaviour of the system allows the wave to travel through the whole bar and be reflected on the other side, returning as a traction wave.

Change in the critical stress

We consider the standard AT1 model. The first possible solution for the difference in the free surface velocity could be a simple change in the value of the constant w 1 , which is equivalent to changing the critical stress σ c , as

σ 2 c =E 0 w 1 .
By changing the value of w 1 , we can see its influence in the model. Since we don't want to change the value of the dissipated energy, changing w 1 forces us to also change .

We can see in Figure 4.5 that increasing the value of σ c causes a change to the free surface velocity. For higher values of σ c , we see that the velocity, after the initial peak, has lower values for the minimum, but keeps roughly the same values for the following peaks. In the experiments, we observe that the peaks are lower after rupture. This approach, however, is not recommended for two reasons: the first one is that by changing the critical stress, it will no longer coincide with the one of the quasi-static experiment. The second reason is that, to have a good approximate curve for each experiment, we would need a different value of σ c for each case, even though the same material was used.

Dependency of the deformation speed ε

Since for each experiment the velocity of the bar was different, we decided to use a law for w 1 that depended on the value of strain rate ε.

We proposed a law of the form

w 1 ( ε) = 1 + ε2 ε 2 ref w 1 ,
for some reference constant ε ref such that the results would be according to the experiments. This approach allows us to find the correct value for the critical stress when ε tends to zero and we obtain the quasi-static case.

The value of ε can be found in two ways: by finding the derivative of the velocity in space, or by dividing by ∆t the difference of the strain in two successive instants.

As we can see in Figure 4.6, this method gives us good results for the G692 test, and the critical stress is unchanged for a quasi-static simulation. This method, however, is not viable because for each experiment, we need a different value of ε ref . By taking ε ref =5000s -1 , we obtain a good velocity curve for the G692 test, but the same is not true for the G672 test, as shown in

Two-dimensional simulations

All the simulations discussed so far were of a one-dimensional bar. We tried the same simulations in 2D, imposing a uniform pressure at x=0.

We tried the plane-strain and the plane-stress cases, but there were no significant changes in results.

It is worth noting that in 2D and 3D, we have many more choices of the traction-compression asymmetry criterion and these choices were not explored in this work. We imposed the condition Trε>0 in order to have damage, so we could test the influence of moving from a one-dimensional to a multi-dimensional model.

When comparing the 1D to the 2D results, we see that we obtain the same damage region and that the free surface velocity is still far from that of the experiments.

Other laws for w(α)

We tried adding more terms to w(α). For instance, if we write w(α)=w 1 α+w 2 α 2 , we obtain some interesting results. Using this model, we can have a numer- ical curve that approaches the experimental one much better. Even though the results are still not perfect, we obtain better results for all the cases. It is important to notice that a change in w 2 does not change the critical stress (which depends on w (α=0), so we can change this value as much as we want and the model will still be compatible with the results obtained in the quasi-static test.

In Figures 4.10 We conclude by observing that the calibration is done only once and works for every case, as opposed to what happened in the previous sections.

Adding a dissipation term

We add, to the discrete energy, a term proportional to ∆t α i -α i-1 ∆t 2

. By doing so, we obtain the discrete version of a damage criterion that depends also on α.

So far, this is the most effective approach. We can see in Figure 4.12 that the velocity at the free surface decreases and, if we consider only its amplitude, is close to that of the experiments. This amplitude is related to the energy that remains between the cracks and the free surface of the specimen, meaning that we have conserved the correct amount of energy and, therefore, we have the appropriate dissipation.

In particular, once we obtain the value of w 1 from the quasi-static experiments, we define

w(α) = w 1 α + k 1 ∆t α i -α i-1 ∆t 2 . ( 4.1)
By doing so, the damage criterion depends also on α, causing the model to depend implicitly on the strain rate.

This new form of dissipation is compatible with the results obtained in the quasi-static test, as the last term disappears as the strain rate tends to zero.

It is important to notice that, with this approach, we obtain the correct dissipation for all the spalling tests for the same dissipation function as we can see in Figures 4.12

Conclusion of the Chapter

In this chapter, we proposed changes to the initial damage model in order to better approach the experimental results, in particular, the free surface velocity.

The results obtained are, qualitatively, close to the experiments. When considering only the region where damage appears, we are able to find a good approximation using our simulations. As for the free surface velocity, we are able to control the amplitude and the position of the the simulated velocity without changing the critical stress of the material. The next step is to find a solution for the difference in the period of the oscillations.

For the one-dimensional problem, the most effective change is in the function w(α), where adding a quadratic term or a term that depends on α yields good results.

For the 2D or 3D model, there is still room for future work if we consider the traction-compression asymmetry.

Regularization of the Plastic Strain

The main motivation behind the regularization of the damage field comes from the fact that brittle damage with softening behaviour are not viable, as there are problems related to the existence or unicity of solutions, or absence of stable configurations. The same is true for other phenomena in solid mechanics, such as plasticity. In particular, we are interested in finding a model where the plastic strain appears in a region of the structure, representing, for instance, the necking of bar, instead of being localized in a zero-thickness band.

In this chapter, we propose different ideas for the regularization the plastic strain:

1. temperature-plasticity coupling: when plasticity occurs, energy is dissipated as heat, increasing the temperature of the bar. The heat equation is then used to calculate the evolution of the temperature. The main objective is to discover if the regularization character of the heat equation is enough to regularize the plastic strain;

2. gradient plasticity: to solve the problem of localization of the plastic strain in softening materials, we follow the same approach as in the gradient damage models and we add a gradient term to the total energy.

These models are then applied in the setting of a quasi-static traction test to see which are viable and which are not.

If we consider that there are only small oscillations in temperature, and define C 0 := ρ 0 c 0 T 0 , we obtain a very simple model of temperature-plasticity coupling:

     σ = E(ε -ε p ) |σ| ≤ σ Y (T ) C 0 Ṫ = k∆T + σ Y (T ) εp .
(5.14)

To calculate the evolution of (u, ε p , T ), we use a variation of the standard alternate minimization algorithm, where we solve one problem at time, while keeping the other variables fixed, until all three variable converge. The FEniCS library was used.

Suppose that, for the discrete solution, we know u(t -∆t), ε ( t -∆t) and T (t -∆t).

At the instant t, we have:

• displacement problem: find u j such that σ = 0, where we control the displacement at the extremities;

• plasticity problem: find (ε p ) j such that σ = E((u j ) -(ε p ) j ) ≤ σ Y (T j-1 );

• temperature evolution:

C 0 T j -T (t -∆t) ∆t = k∆T j + σ Y (T j ) (ε p ) j -ε p (t -∆t) ∆t .
(5.15)

• Repeat this process until we obtain convergence for u j , (ε p ) j , and T j . We then take u(t) := u j , p(t) := (ε p ) j and T (t) := T j .

We fix one of the extremities of the bar, pull the other extremity with constant speed v. The boundary conditions used in the temperature problem will be detailed in the results section.

Dimensionless problem

In order to better understand the problem, the first step is to reduce the number of variables available. We consider a bar of length L and define x = x/L. We also consider a length scale U 0 > 0 and characteristic time τ > 0 (to be specified later), and define t = t/τ .

We define the dimensionless displacement and temperature by ũ(x, t) = u(x, t) U 0 and T (x, t) = T (x, t) T 0 .

(5.16)

-Temperature-Plasticity Coupling

We have

du dx = U 0 L dũ dx and dT dt = T 0 τ d T d t .
(5.17)

We begin by normalizing the stress. We define σ 0 := σ Y (T 0 ) and σ( T ) = σ Y (T )/σ 0 . This way, the plastic phase begins when the stress σ reaches one.

We define p(x, t) = L U 0 ε p (x, t).

(5.18)

We also have

σ = E σ 0 U 0 L ( dũ dx -p). (5.19)
Since we are free to specify U 0 , we define U 0 = Lσ 0 /E. Rewriting the evolution of the temperature, we obtain

C 0 T 0 τ σ 0 d T d t = kT 0 L 2 σ 0 d 2 T dx 2 + U 0 Lτ σ Y (T ) σ 0 dp d t (5.20)
or, equivalently,

C 0 T 0 L U 0 σ 0 d T d t = kT 0 τ U 0 Lσ 0 d 2 T dx 2 + dp d t . (5.21)
We can then chose the last free parameters τ as τ = LU 0 σ 0 /kT 0 . We can change notations without generating any confusion to consider the spatial and temporal derivatives to be in respect to the dimensionless position x and time t, and remove the tilde from the dimensionless variables to obtain

     σ = (u -p) |σ| ≤ σ Y (T ) C 1 Ṫ = ∆T + σ Y (T ) ṗ, (5.22)
where the initial temperature is T 0 = 1 and σ Y (1) = 1.

Example 5.1.1. In order to have a reference value for the dimensionless parameters of a realistic material, we consider a 10cm specimen made of steel.

Its density is ρ=8000kg/m 3 and the specific heat capacity c p =500J/kgK, so that the heat capacity is C 0 =ρc p =4•10 6 J/Km 3 . The Young's modulus is E=200GP a and the yield stress σ Y =250M P a. We consider the temperature T 0 =300K and thermal conductivity k=50W/mK.

By using the previous expressions, we obtain

U 0 = Lσ Y E = 0.1m•250M P a 200GP a = 0.000125m (5.23)
Chapter 5 -Regularization of the Plastic Strain

τ = LU 0 σ Y kT 0 = 0.1m•0.000125m•250M P a 50W/mK•300K = 2.08333s (5.24) C 1 = C 0 T 0 L U 0 σ Y = 4•10 6 J/Km 3 •300K•0.1m 0.000125m•250M P a = 3840.
(5.25)

Homogeneous results

We want to study the homogeneous response of this system. In the numerical simulations, this is achieved by controlling the displacement at the extremities (u(0) = 0 and u(L) = t), and imposing T (0) = T (L) = 0.

Analytically, we have the scalars σ, ε, p and T . The equations to be solved are given by (5.22).

The relation of the yield stress and the temperature still needs to be specified. We consider here the particular case

σ Y (T ) = 1 T . ( 5.26)
It is easy to see that we first have an elastic phase while ε ≤ 1.

After the beginning of the plastic phase, we have

ε -p = 1 T (5.27)
and

C 1 Ṫ = ṗ T . ( 5.28)
We multiply this expression by T and integrate from t = 0 to the instant t to obtain

T 2 = 1 + 2p C 1 . ( 5.29)
For a given strain ε in the plastic phase, we use equation (5.27) to find

p - 1 1 + 2p C 1 -ε = 0.
(5.30) Equation (5.30) can easily be solved to give us the answer of the system. We define

f (ε, p) = p - 1 1 + 2p C 1 -ε.
(5.31)

-Temperature-Plasticity Coupling

In order to study the snapback phenomenon when the plastic phase begin, we first observe that f (1, 0) = 0. In this same point (ε, p) = (1, 0), we have that ∂f ∂ε = -1 and ∂f ∂p = 1 -1 C 1 . Using the implicit function theorem, we have that the relation between a small change in the strain and in the plastic strain:

∆ε ≈ (1 - 1 C 1 )∆p.
(5.32)

Therefore, we have two situations:

1. if C 1 > 1, then the plastic strain will increase when the strain increases;

2. if C 1 < 1, we can see the snapback phenomenon.

We see in Figure 5.1 the evolution of the stress when we increase the strain, for the homogeneous problem. Figure 5.2 shows the temperature and the plastic strain. The imposed boundary conditions (T (0) = T (L) = 0) do not provoke any dissipation, and the temperature is always increasing. We see that the curves are smooth and that the analytic and the numeric curves are the same, validating our numerical implementation. In this case, we can see there is a snapback in the beginning of the plastic phase. We do not obtain the same curves analytically and numerically, even after refining the numerical parameters. This is not a surprise, since we obtain a discontinuity in p(ε) and T (ε) when ε = 1. To obtain convergence of the discretization that we used for the evolution of the temperature, given by equation (5.15), it is required that the temperature be continuous. Therefore, a more sophisticated time discretization is needed. One possible idea is to use a scheme that follows the snapback curve.

Non-homogeneous results

We then analyse a non-homogeneous case. We consider the same problem, that is, σ Y (T ) = 1/T , but with a different boundary condition for the temperature. We consider that the extremities are in contact with a surface with constant temperature. We impose T (0)=T (L)=T 0 =1. In this problem, the temperature is dissipated and, in the absence of a heat source, the temperature will become uniform in the bar and equal to 1.

We have no analytic results for this model. Numerically, however, this example works very poorly. For the examples below, we set C 1 =1 and use 100 elements. We remark that increasing the number of elements does not change the results.

These parameters, even though they are not realistic (when compared, for instance, to steel where C 1 ≈4000), are useful to show how the time step dt influences the results.

We see the evolution of the stress when the strain increases in Figure 5.5. We first observe that we do not obtain numeric convergence when changing the time-step. We can see that the stress decreases abruptly at certain values of strain, and then increases linearly until the next drop. The number and the instant where these drops happen changes with the time-step used.

Figure 5.6 shows the maximum temperature and plastic strain as the total strain increases. Again, we observe that there is no convergence. We also notice that temperature and plastic strain jump for some values of strain (same instants where the stress drops). The dissipation of heat is faster than the increase of strain. As a result, the stress does not reach the yield limit and the material is, once again, in elastic phase. This elastic phase lasts until the yield stress is reached again, and the cycle begins once again.

As a final example of why this model is not suitable, we consider a more realistic case by taking C 1 =4000. We show the plastic strain in Figure 5.7. We recall that the main objective of this section was to use the temperature as a way of controlling the localization of the plastic strain in softening materials. By taking the C 1 parameters to be that of a steel, the regularizing effect is obtained in excess, leading to an almost uniform profile, with some small variations near the boundaries.

In conclusion, the regularization using the dissipation of the temperature does not work well. The proposed model is very simple, but the numerical implementation needs some attention. For the homogeneous case, the model works well only when there is no snapback. For the non-homogeneous problem, the results are even worse, as we did not obtain convergence using our numerical scheme.

Gradient Plasticity

We follow the same idea of regularization used in gradient damage and we add, to the total energy, a term that depends on the gradient of the plastic strain. As we have seen in the introduction chapter, the quasi-static evolution of the plastic strain can be written as a problem of energy minimization. We present here the construction of a gradient-plasticity model as natural extension of the standard plasticity theory, but with a term that prevents localization in infinitely thin bands.

Even though this is not our main motivation, one important feature of this approach is the ability to better predict size effects, as no length scale appears in conventional plasticity theory. We will suppose that the bar is always in traction so that p = ε p ≥ 0, allowing us to remove the cumulated plasticity of the model.

We consider the general case and define the energy

E(ε, ε p ) = Ω 1 2 E(ε -ε p ) 2 + w(ε p ) + σ C 2 (ε p ) 2 .
(5.33)

As in gradient damage, we have a term of elastic energy, a term of energy dissipated by standard plasticity, and a last term consisting of the squared of the gradient of the plastic strain, weighed by a characteristic length and a normalization factor σ C , having the unit of a stress.

The standard Von-Mises criterion is equivalent to

w(ε p ) = σ Y ε p .
(5.34)

Localized Solution

We start the analysis ot this model by construction a localized solution p on the interval I = (-D, +D) when the displacement at the extremities is controlled. By differentiating the total energy with respect to the plastic strain, we can obtain (5.37)

-σ + w (p) -2σ C 2 p = 0. ( 5
The stress σ and the maximum plastic strain p max are related by w(p max ) = σp max .

(5.38)

We will suppose that the strain is symmetric and the maximum value is attained in x = 0. We can obtain x from the plastic strain p:

x = pmax p σ C -σp + w(p)
dp.

(5.39)

The value of the extremity D is thus given by

D = pmax 0 σ C -σp + w(p) dp. ( 5

Numerical Example

In order to test the model's ability to regularize, we consider a softening material.

We consider the particular case

w(ε p ) = σ Y ε p 1 + ε p /p C , ( 5.42)
for the constants σ Y and p C . In order to obtain a localized solution away from the extremities, we fix ε p (-L) = ε p (L) = 0 and run the numerical simulations using theFEniCS library.

We consider the bar [-L, L], with L=0.5, and with the parameters E=1, =0.05, σ C =0.5, p C = 1.0 and σ Y =1.0.

At the time instant t, we impose the displacements u(±L)= ± t 2 . The numerical problem consists of minimizing the total energy with respect to the strain and the plastic strain. For that, we use the well-known alternate minimization algorithm, where we minimize the energy with respect to one of the variables, while the other is fixed, until we obtain convergence. It is important to recall that we consider a bar under an increasing load, so that we can consider that εp ≥0 at all times. Numerically, this is done by considering that ε p (t i )≥ε p (t i-1 ).

As opposed to the damage problem, we noticed that we need a great number of iterations in order to obtain convergence. For a one-dimensional mesh with 400 elements with a tolerance of 10 -6 , a few hundred iterations are needed before convergence is reached.

In the previous section, we obtained the analytic expression for the profile of the plastic strain. This expressions could be easily solved by direct numerical integration. In Figure 5.8, we see that the profiles obtained using the finite-element simulation and the numerical integration coincide.

The support of the localized solution increases with time, until it reaches the borders and there are no longer localized solutions (Figure 5.9). We see the stress as a function of the displacement in Figure 5.10. As expected, we first have a linear relation, until σ reaches the critical value σ Y . There is then a sudden decrease in stress, which then continues to slowly decrease. This is true while a localized solution exists. Once the support of the solution is large enough, making it occupy the whole bar, the stress increases as the strain increases.

In conclusion, we have a model that allows the evolution of an elastic phase, followed then by the evolution of the plastic strain. The main feature of this model is its capacity to control the problem of localization in thin bands by using a gradient term.

As for the boundary conditions, since we wanted to study localization profiles in the center of the bar, they were taken to be zero and its importance could be ignored. As the system evolves, however, the support of the plastic strain becomes the whole bar, and the chosen boundary conditions are of extreme importance. A possible continuation to this work is understanding their influence in the profile obtained and what is their physical meaning.

Another possible continuation is the extension of this model to higher dimensions. As we have seen, we have considered the plastic strain to be equal to the cumulated strain. Since this hypothesis no longer holds in higher dimensions, the extension of this model to a plaque or a cube is not trivial.

-Gradient Plasticity Coupled with Damage

.9: Profile of the plastic strain at t=4.0.

Gradient Plasticity Coupled with Damage

Damage-plasticity coupling

A natural continuation of the gradient plasticity is its coupling with damage. The theory of gradient damage has one internal length, responsible for determining the thickness of a crack. In the same way, our model of gradient plasticity has also one internal length used to define the support of the damage zone.

We then propose a model of damage-plasticity coupling using two internal lengths α and p . The approach remains the same as in the other cases: we define a suitable form of energy and study the quasi-static evolution of the system by minimizing this energy, at each instant, with respect to the strain, the damage field and the plastic strain.

We consider the energy of a bar [0, L] under traction:

E(ε, α, p) = L 0 1 2 E(α)(ε -p) 2 + w(α) + w 1 2 α α 2 + σ Y (α)p + 1 2 γ(α)λ 2 p p 2 dx.
(5.43) The function γ(α) has the unit of a stress and is used for normalization. We are interested, in particular, in the behavior of the system when γ is constant or when γ(α)→0 as α→0.

We assume the loading is increasing. For this reason, the plastic strain is also always increasing ( ṗ≥0). Numerically, we solve this problem using an alternate minimization algorithm:

• For the instant t i , impose the displacements u(0)=0 and u(L)=t i :

minimize E with respect to ε, keeping α and p fixed;

minimize E with respect to p, keeping ε and α fixed;

minimize E with respect to α, keeping ε and p fixed;

check convergence on all variables.

• Go to the next time-step.

We remark, once again, that the irreversibility conditions of α and p are imposed as α i ≥α i-1 and p i ≥p i-1 .

The final question that needs to be addressed is that of the boundary conditions. In the following simulations, we allowed the material to have an elastic phase, followed by the uniform evolution of plasticity. The damage process is the main agent in localization. We therefore consider p =0 and α=0 at the extremities of the bar.

Results

We are interested in the behaviour of the dimensionless problem. More specifically, we want to know how the two characteristic length of the system

-Results

interact. We are going to consider a unitary bar and E(α)=(1-α) 2 , w(α)=α, σ Y (α)=(1-α) 2 σ P and, unless otherwise stated, σ P =0.5. For this set of parameters, it is easy to verify that the plastic phase begins at ε=0.5, and the damage phase begins when ε=1.25 and p=0.75.

The first steps consists of validating our numerical implementation. We consider two cases where we know what should happen:

1. when σ P >1, we obtain the standard gradient damage model; 2. when σ P <1, we obtain the uniform plastic evolution before the beginning of the damage phase.

After verifying that the code works for these two cases, we study the behaviour of the system for different values of α and p . We remark that, for each of the following cases discussed, different time-steps and mesh sizes were tested, always yielding the same results. For this reason, the numerical parameters used will be omited.

We consider the case α = p =0.1. This case seems to gives the expected result, that is, the damage profile and the plastic strain are spread in the center of the bar, as shown in Figures 5.11 and 5.12. In particular, we remark that the plastic strain evolves uniformly until the beginning of the damage phase, where it is equivalent to 0.75, and then evolves only in a interval near the center. We see in Figure5.13 the stress-strain relation for this model. We can clearly identify the elastic phase for ε<0.5, the homogeneous plastification for 0.5 < ε<1.25 and then the damage phase. During the damage phase, we see that the damage first increases slowly and then instantly jumps to one, representing a completely broken state with no stress left. p ≈0.09. We go from a model where there is no localization of the plastic strain to a model where the plastic strain evolves in an interval of approximate length 0.4. Figure 5.17 shows the influence of p when α =0.1. There is no longer two different types of behaviour, as in the previous case. When we decrease p , the region of plastic localization decreases continuously. By decreasing p and the mesh size, the support of the plastic zone converges to zero and the overall behaviour is the same of the model without the gradient of the plastic strain.

-Results

Conclusion

Due to the complexity of the differential equations involved, we were not able to obtain analytic expressions describing the evolution of the system. For α =0.1, the system behaves in the expected way, that is, there is a zone where there is damage localization and a zone where there is localization of the plastic strain. Furthermore, decreasing p causes a decrease in the zone of plastic localization that approaches zero as p tends to zero.

The behaviour of the system when α =0.05 is completely different. The dependency with respect to p is more complicated, as there is localization of the plastic strain for small values of p . The causes of this behaviour is not clearly and needs to be further investigated. By testing different mesh

-Conclusion

sizes and times-steps, the origin of this difference of behaviour does not seem to be numerical. Our guess is that this is a question of instability leading to a snapback, but so far we have no mathematical proof and further work is needed.

Overall, the model gives us promising results. For the functions taken as E(α), w(α) and σ Y (α), along with the boundary conditions, we always have an elastic phase and a homogeneous plastification phase. We can adjust the model's behaviour after the beginning of the damage phase by changing the characteristic lengths α and p , allowing us to obtain a a material that will break instantly, or a material that breaks slowly and allows the development of a localized plastic zone.

Conclusion and Future Work

The main objective of this work was the study of the fragmentation of a metallic shell. In order to achieve this goal, several intermediate steps were necessary: the construction of a damage model; its numerical implementation; calibration of the model parameters using experimental data and some analytical works.

Gradient damage models have been originally proposed for quasi-static brittle damage evolution, but have been extended to other scenarios, such as for ductile materials or dynamic loadings. In this work, we considered a model that couples these last two phenomena. This was done by writing a suitable form of energy that couples the elastic energy with the energy dissipated by the damage and plasticity process. Using this energy and the action of the system, we could obtain all the equations necessary to describe the dynamic ductile model: the equations of dynamics, the plasticity criterion and the damage criterion.

We then described, in Chapter 2, all the numerical aspects related to these models. The plasticity problem is solved in a straight-forward fashion, by locally comparing the stress to a threshold. The damage problem, as it is a global problem, requires us to assemble a large system that is then solved using the PETsC library. To solve the equation of dynamics and find the displacement field, we used two main schemes: the explicit Newmark scheme and a new version of the Generalized Midpoint Rule Scheme. In summary, these two schemes have a quadratic convergence rate for purely elastic or brittle damage problem, but a poorer convergence rate when plasticity is considered. Numeric experiments showed that both schemes can be considered to conserve the total energy, but we have no theoretical results supporting our claim. For the Genezalized Midpoint Rule scheme, since we have to solve each problem several times before converging, along with the fact that the dynamic problem is implicit, each iteration takes more time. This fact is, however, counterbalanced by the fact that there are no restrictions on the time-step. We were also able to prove stability for this scheme.

After the validation of our implementation, we were able to obtain some qualitative behaviours for our model. When considering the standard bar under traction, we see a very different behaviour depending whether the loading is quasi-static or dynamic. For the quasi-static scenario, only one crack appears. For the dynamic model, multiple cracks can appear, depending on the parameters used. We also recall that the system keeps evolving and other damage profiles can appear and develop, even after an initial crack appeared, due to the waves that keep propagating in the system.

A natural question to be posed in regards to this dynamic model is whether it converges to the quasi-static model when taking sufficiently small strain rates. When comparing the results obtained for a quasi-static loading and for a dynamic loading, the damage profiles are not the same. When decreasing the strain rates the damage profile seem to converge, but to a new profile that dissipates roughly 30% more energy than the one obtained under the assumption of a quasi-static test. The reason is not yet completely understood, but it might be related to the fact that, in the quasi-static experiment, a part of the energy is dissipated by the snapback process, leaving less energy to be dissipated by the crack.

The last aspect concerns the direction of cracks. The quasi-static traction test of a brittle specimen causes the crack to be perpendicular to the loading. When plasticity is added, the cracks forms an angle with the traction direction depending on the thickness of the specimen and on the material parameters used. The orientation of the crack can change of up to 45 • . In dynamics, this remains true. It was observed, especially in the expanding ring, that the cracks are in the radial direction for brittle materials and are inclined for ductile materials.

Since the study of the fragmentation of a complete shell is too costly in terms of computer power, we studied the case of a ring under expansion. This ring was modelled as a periodic one-dimensional bar and as a real three-dimensional ring. Due to the symmetry of the three-dimensional ring, we expected the response to depend only on the radius, and not on the angle. However, this is not what we obtained in the simulations, as some radial cracks appeared. We then focused on the questions of why and how many cracks appear in this ring. We have shown that the localized cracks result from the growth of small perturbations that are amplified of 10 10 times or even more. Using this idea, we were able to find the moment where these cracks develop and analytically show that, as a first approximation, the number of cracks depends only on the internal length.

After establishing a model, implementing it and testing it for the case of a ring, the last necessary step in order to run more realistic simulations 140

-Conclusion

in an industrial context is to calibrate the model parameters. Since this is an initial study, we focused on a brittle material, as a ductile material would require more parameters to be studied. Experiments have noticed that materials don't have the same properties for different strain rates. In Chapter 4, the spalling tests were used to propose material parameters that work well for both quasi-static and dynamic loadings by comparing the changes in our model to the experimental data obtained. We found that the best two approaches were to add, to the term of dissipated energy w(α), a term proportional to α 2 or a term proportional to α, so that there was no change to the critical quasi-static stress.

The final chapter consists of the study of other forms of regularization, but now applied to the plastic strain, so that it is no longer localized in a zero-thickness band. The first attempt was to use the dissipative properties of the temperature field to regularize the model. From a numerical point of view, this model was not optimal, as we had convergence problems.

We then talked about removing the temperature from the equations and adding a term depending on the gradient of the plastic strain. This approach allows us to control the growth of sharp profiles and work well numerically. The final test was coupling this gradient plasticity model to the gradient damage model. This new model has two characteristic lengths that dictate the size of the damage and the plastic strain profile. Increasing these characteristic lengths causes the damage and plastic strain profiles to be more spaced out in the bar, but overall, the behaviour or this model is not yet completely understood.

In terms of future work, there is clearly room for more simulations and comparison with experiments, as the main focus of this thesis was the theoretical understanding of damage models coupled with plasticity and dynamics. In this context, there are mainly two questions that remain open:

1. Why doesn't this model converge to the quasi-static model? Is it only the snapback process, or is the strain rate still too elevated? If we are able to run simulations with a sufficiently low strain rate, will we obtain only one damage profile?

2. After a damage profile appeared and evolved to the point of reaching one, there is a thickening of the region where α=1. Des résultats qualitatifs ont ensuite pu être obtenus, comme le nombre et la direction des fissures, ainsi que la convergence vers le modèle quasi-statique. Afin de mieux comprendre l'influence de chaque paramètre du modèle, nous avons étudié analytiquement le problème. A partir de l'observation de l'amplitude des perturbations, nous avons pu décrire comment obtenir une approximation analytique du nombre de fissures dans le cas d'un anneau en expansion. Cependant, pour être capable de simuler des problèmes réalistes, il est nécessaire de calibrer les paramètres du modèle. Nous nous sommes intéressés plus particulièrement au cas des matériaux fragiles. Les données expérimentales ont été obtenues par une série d'expériences réalisée par le CEA. Afin d'empêcher la localisation de la déformation plastique dans des bandes d'épaisseur nulle, d'autres formes de régularisation ont été étudiées, comme par exemple, l'utilisation des propriétés dissipatives du champ de température. Enfin, nous avons conclu ce travail en proposant des modèles de plasticité où l'énergie dépend aussi du gradient de la déformation plastique (modèles de plasticité à gradient).

Title: Variational approach to dynamic fracture and applications to the fragmentation of metals and ceramics Keywords: damage, dynamic fragmentation, plasticity Abstract: The main objective of this work was the study of the fragmentation of a metallic shell. This thesis is divided into four parts: construction of a damage model, numerical implementation, calibration of the model parameters using experimental data and analytical works. In this work, we considered a model that couples the standard gradient damage models with plasticity and dynamics. Using the energy and the action of the system, we could obtain all the equations necessary to describe the dynamic ductile model: the equations of dynamics, the plasticity criterion and the damage criterion. We then detail the numerical implementation of these models. Some qualitative behaviours are then obtained, such as the number and the direction of cracks, and the convergence to the quasi-static model. In order to better understand the influence of the parameters, we studied the problem analytically. By studying the amplitude of the perturbations, we describe how to obtain an analytic approximation for the number of cracks in a ring under expansion. In order to run realistic simulations, it is needed to calibrate the material parameters. We focus here on a simple case of brittle materials. The experimental data were obtained in a series of shockless spalling tests performed by the CEA. We also study other forms of regularization, now applied to the plastic strain, avoiding localization in zero-thickness bands. We considered using the dissipative properties of the temperature field to regularize the model. Finally, we conclude with plasticity models where we add a term depending on the gradient of the plastic strain (gradient plasticity models).

Figure 1 . 2 :

 12 Figure 1.2: Damage (dashed red), normalized stress (blue) and plastic strain (green dots), according to example 1.2.2.

Figure 1 . 3 :

 13 Figure 1.3: Brittle damage. E(α)=(1 -α) 2 and w(α)=α.

Figure 1 . 4 :

 14 Figure 1.4: Ductile damage. E(α)=σ Y (α)=(1 -α) 2 and w(α)=α.

Figure 1 . 5 :

 15 Figure 1.5: Brittle damage. E(α)=(1 -α) 2 and w(α)=α 2 .

Figure 1 . 6 :

 16 Figure 1.6: Ductile damage. E(α)=σ Y (α)=(1 -α) 2 and w(α)=α 2 .

( 2 )

 2 If |σ * |>σ Y , we put e=hσ * /|σ * |, with h>0. Then

( 3 . 1 )

 31 solve the plasticity problem; (3.2) solve the damage problem; (3.3) stop when the alternate minimisation converged for damage α i+1

Figure 2 . 1 :

 21 Figure 2.1: Error in the displacement for different values of time-step. The dashed red line represents the quadratic convergence rate.

  and w(α)=α. The errors in the displacement are shown in Figure 2.2 and the errors in damage are shown in Figure 2.3. By comparing the errors to the dashed red line, we verify the quadratic convergence of displacement and damage of both schemes.

Figure 2 . 2 :

 22 Figure 2.2: Error in the displacement for different time-steps at t=0.2.

Figure 2 . 3 :

 23 Figure 2.3: Error in the damage profile for different time-steps at t=0.2.

Figure 2 . 4 :

 24 Figure 2.4: Representation of the two regions of the bar.

Figures 2 .

 2 5 and 2.6 were obtained using 1000 elements.

Figure 2 . 5 :

 25 Figure 2.5: Example of stress and damage profiles on the bar. In these cases, t=0.22.

Figure 2 . 6 :

 26 Figure 2.6: Evolution of X(t).

Figure 2 . 7 :

 27 Figure 2.7: Function G(λ) for E 0 =w 1 =1

Figure 2 . 8 :

 28 Figure 2.8: Damage profile for D(α) = 0.1 for E 0 =w 1 =1.

Figure 2 . 9 :

 29 Figure 2.9: Damage profiles for E 0 =w 1 =1.

Figure 2 .

 2 Figure 2.10: Stress-displacement relation for E 0 =w 1 =1.

Figure 2 . 11 :

 211 Figure 2.11: Damage profile for a brittle material with E(α)=(1 -α) 2 and w(α)=α. We considered =0.1 and ε = 0.1.

Figure 2 . 12 :

 212 Figure 2.12: Damage profile for a brittle material with E(α)=(1 -α) 2 and w(α)=α. We considered =0.1 and ε = 1.0.

Figure 2 . 13 :

 213 Figure 2.13: Analytic damage profile for a brittle material under quasi-static loading (red) and damage profile under dynamic loading obtained numerically (blue).

Figure 2 . 14 :

 214 Figure 2.14: Analytic damage profile for a brittle material under quasi-static loading (red) and damage profile under dynamic loading obtained numerically (blue) for =0.5.

Figure 2 . 15 :

 215 Figure 2.15: Damage profile for a brittle material ( =0.1 and very low strain rate). We see that α=1 in an interval, indicating a thickening of the crack.

Figure 2 . 17 :

 217 Figure 2.17: Damage profile for a ductile 3D plaque under traction.

Figure 2 . 18 :

 218 Figure 2.18: One-dimensional example of the average and jump between two elements.

Figure 2 . 19 :

 219 Figure 2.19: Displacements obtained using the standard P1 elements and the DG formulation for the quasi-static brittle problem.

20

 20

  and 2.21 to compared the damage profile. It is clear the the profile obtained using the DG formulation provides a better approximation to the analytic one. We remind that only 20 elements were used.

Figure 2 . 20 :

 220 Figure 2.20: Damage profile obtained using the standard P1 elements for the quasi-static brittle problem.

Figure 2 . 21 :

 221 Figure 2.21: Damage profile obtained using the DG formulation for the quasistatic brittle problem.

  Figure 2.22.

Figure 2 . 22 :

 222 Figure 2.22: Displacements obtained using the standard P1 elements and the DG formulation for the dynamic problem of a bar with a sudden crack in the middle (t=0.25).

  .1.

Figure 3 . 1 :

 31 Figure 3.1: Example of fracture of a cylinder using the dynamic gradient damage approach.

. 22 )•

 22 The time evolution of the displacement u is given by ü = a(α)ε(u) , (3.23) under the imposed boundary conditions u(x=0, t) = 0 u(x=1, t) = ε0 t. (3.24)

. 27 )

 27 We obtain =1.2•10 -5 m and ˜ = /L=1.2•10 -4 . The values of T and U 0 are • 10 6 210 • 10 9 m = 4.7 • 10 -4 m. (3.29)

Figure 3 . 2 :

 32 Figure 3.2: Example of behaviour of evolution of the system. We see the damage profile α * for different time instants. We use the dimensionless parameters ε=0.5 and =0.02.

Figure 3 . 3 :

 33 Figure 3.3: We fix ε0 = 0.5 and see the number of cracks for different values of .

. 35 )Figure 3 . 4 :

 3534 Figure 3.4: Influence of ε0 for a fixed value of .

Figure 3 . 7 :

 37 Figure 3.7: N = 70, amplitude 8

Figure 3 . 10 :

 310 Figure 3.10: Value of y for τ << 1/ √ 3 for λ=100. τ close to 1/ √ 3 We write τ =1/ √ 3 + T . Thus g(τ ) = ( √ 3τ -1)( √ 3τ + 1) τ 4 (τ 2 + 1) ≈ 27 √ 3 2 T (3.47) and d 2 y dT 2 ≈ 27 √ 3 2T y(T ).(3.48)

Figure 3 . 11 :

 311 Figure 3.11: Value of y for τ close to 1/ √ 3 for λ=100.

Figure 3 . 12 :

 312 Figure 3.12: Value of y for τ >> 1/ √ 3 for λ=100.

Figure 3 . 13 :

 313 Figure 3.13: Function τ ψ (τ ) for large values of λ (λ=100 in the image).

Figure 3 . 14 :

 314 Figure 3.14: Number of cracks for ε=1 and different values of . The analytic approximation is obtained by the expression L/2π .

Figure 3 . 15 :

 315 Figure 3.15: Value of α 0 for different dimensionless strain rates ε.

Figure 3 . 16 :

 316 Figure 3.16: Illustration of a generic mesh used for the simulations.

Chapter 3 - 4 :

 34 Application: Explosion and Fragmentation of a Ring ε expected number of cracks cracks obtained in the simulation 0Influence of the strain rate in the number of cracks for =0.063.

Figure 3 . 17 :

 317 Figure 3.17: Damage profile for ε=0.5 and =0.0157. We obtained 58 cracks (64 expected).

  -3.21.

Figure 3 . 18 :

 318 Figure 3.18: Evolution of f (t) for k=1, σ p =0.5, =0.01 and ε=1.0. Amplitude 1.15.

Figure 3 . 19 :

 319 Figure 3.19: Evolution of f (t) for k=10, σ p =0.5, =0.01 and ε=1.0. Amplitude 100.

Figure 3 . 20 :

 320 Figure 3.20: Evolution of f (t) for k=100, σ p =0.5, =0.01 and ε=1.0. Amplitude 6•10 15 .

Figure 3 . 21 :

 321 Figure 3.21: Evolution of f (t) for k=1000, σ p =0.5, =0.01 and ε=1.0. Amplitude 10 24 .

Figure 3 . 22 :

 322 Figure 3.22: Damage on a ductile cylinder using 2000 elements on the perimeter.

Figure 3 . 23 :

 323 Figure 3.23: Damage on the external surface for a ductile cylinder using 500 elements on the perimeter. We obtain 48 profiles.

Figure 3 . 24 :

 324 Figure 3.24: Damage on the external surface for a ductile cylinder using 2000 elements on the perimeter. We obtain 68 profiles.

Figure 4 . 2 :

 42 Figure 4.2: Analysis of the specimen for experiment G692 (Zinszner et al. [53]).

4. 3 -

 3 Change in the critical stress

Figure 4 . 3 :

 43 Figure 4.3: Damage profile for the AT1 model in the G692 experiment.

Chapter 4 -

 4 Shockless Spalling

Figure 4 . 4 :

 44 Figure 4.4: Velocity of the free surface for E(α)=E 0 (1 -α) q and w(α)=w 1 α.

Figure 4 .

 4 7).

Figure 4 . 5 :

 45 Figure 4.5: Velocity of the free surface for the AT1 model and different values of the critical stress.

Figure 4 . 6 :

 46 Figure 4.6: Velocity of the free surface when w(•) depends on the strain rate (ε ref =5000s -1 ).

  and 4.11, we considered w(α)=w 1 (α + 12α 2 ).

4. 7 -

 7 Adding a dissipation term

Figure 4 . 7 :

 47 Figure 4.7: Velocity of the free surface when w(•) depends on the strain rate for different values of ε ref .

  and4.13.

Figure 4 . 8 :

 48 Figure 4.8: Damage profile in 2D.

Figure 4 . 11 :

 411 Figure 4.11: Velocity of the free surface after adding a quadratic term to w (experiment G692).

Figure 4 . 12 :

 412 Figure 4.12: Velocity of the free surface after adding a dissipative term to w (we consider k 1 =0.15MPa) -experiment G672.

Figure 5 . 1 :

 51 Figure 5.1: Stress-strain relation for C 1 =2 in the homogeneous temperatureplasticity coupling model.

3

 3

  and 5.4 the results for the problem when C 1 = 0.1.

Figure 5 . 2 :

 52 Figure 5.2: Evolution of the temperature and of the plastic strain for C 1 =2 in the homogeneous temperature-plasticity coupling model.

Figure 5 . 3 :

 53 Figure 5.3: Stress-strain relation for C 1 =0.1 in the homogeneous temperature-plasticity coupling model.

Figure 5 . 4 :

 54 Figure 5.4: Evolution of the temperature and of the plastic strain for C 1 =0.1 in the homogeneous temperature-plasticity coupling model.

Figure 5 . 5 :

 55 Figure 5.5: Evolution of the stress for C 1 =1 in the non-homogeneous temperature-plasticity coupling model for three different values of dt.

Figure 5 . 6 :

 56 Figure 5.6: Maximum temperature and plastic strain for C 1 =1 in the nonhomogeneous temperature-plasticity coupling model for three different values of dt.

. 35 )Figure 5 . 7 :

 3557 Figure 5.7: Maximum temperature and plastic strain for C 1 =1 in the nonhomogeneous temperature-plasticity coupling model for three different values of dt t=1.51.

. 40 )

 40 It is clear that the support of the region depends linearly on . Finally, we can write the displacement at the extremity L of the bar as a

Figure 5 . 8 :

 58 Figure 5.8: Profile of the plastic strain at t=1.5.

Figure 5 .

 5 Figure 5.10: Stress-strain curve for the gradient plasticity model.

Figure 5 . 11 :

 511 Figure 5.11: Damage profile for α = p =0.1 localized in the center of the bar for t=1.30.

Figure 5 . 14 :

 514 Figure 5.14: Damage profile for α = p =0.05 after fracture. The interval where α=1 has the size of one element of the mesh.

Figure 5 . 15 :

 515 Figure 5.15: Plastic strain for α = p =0.05 after fracture.

Figure 5 . 16 :

 516 Figure 5.16: Size of the zone where the plastic strain localizes α =0.05 and different values of p .

Figure 5 . 17 :

 517 Figure 5.17: Size of the zone where the plastic strain localizes α =0.1 and different values of p .

Table 3 .

 3 3: Influence of the internal length in the number of cracks for ε=0.5.

	-3 .

  Analytically, we have shown that this should not happen. Future work should be devoted to understand what is happening. Le principal objectif de ce travail est l'étude de la fragmentation d'enveloppes métalliques. Cette thèse est divisée en quatre parties : la construction d'un modèle d'endommagement, l'implémentation numérique, la calibration des paramètres du modèle en utilisant des données expérimentales, et des travaux analytiques. Tout d'abord, nous avons considéré des modèles qui couplent les modèles d'endommagement classiques avec la plasticité et la dynamique. En utilisant l'énergie et l'action du système, nous avons obtenu toutes les équations qui décrivent le modèle dynamique et ductile : l'équation de la dynamique,

	le	critère	de	plasticité	et	le	critère
	d'endommagement. Nous avons ensuite détaillé
	l'implémentation numérique de ces modèles.

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France Titre : Approche variationnelle de la rupture dynamique et applications à la fragmentation des métaux et céramiques Mots clés : endommagement, fragmentation dynamique, plasticité Résumé :

2.6 -Discontinuous Galerkinevolution problem. Discontinuities in the displacement are possible exclusively when α=1, but displacements become irregular before we reach this phase.

Remerciements

2. E(u, α) ≤ E(u, β), for all β ≥ α admissible; 3. D(u, α) = t fixed.

Mathematical Notions

We remember here (without much rigour) some results that were useful in our calculations. Details can be seen in Allaire [START_REF] Allaire | Analyse numérique et optimisation: Une introductionà la modélisation mathématique et à la simulation numérique[END_REF].

Let V be a Hilbert space and we define the subset K ⊂ V using the constraints F = {F i } 1≤i≤n :

(2.122)

Suppose we want to solve the constrained optimization problem inf v∈K J(v).

(2.123)

We first define the Lagrangian L : (V × R n ) → R associated to this problem as L(v, p) = J(v) + p • F (v).

(2.124)

We say that (u, p) ∈ V × R n is a saddle point or minimax point if

125)

Theorem : Suppose that (u, p) is a saddle point of L. Then u ∈ K and u is global minimum of J in K.

For v ∈ V and q ∈ R n , we set

(2.126)

Theorem : The pair (u, p) is a saddle point of L if and only if

(2.127)

Numerical Solution

As we have seen in section 3.1.2, the brittle problem depends only on the internal length and the strain rate. For this reason, as an initial study, we are going to fix all parameters equal to one (with the exception of , t and k). We take =0.02 and k=2πN and we look at the evolution of the mode N between the instant where damage begins to increase (t 0 =1) and t f =5.

For that, we have to solve the equation f (t)=g N (t)f (t). The solution will be denoted f N . We can solve these equations numerically, with initial values f (0)=1 and f (0)=0. The behaviour changes drastically depending on N :

• for small values of N , f N has an amplitude of order 10 5 ;

• for medium values of N (N ≈10), f N has an amplitude of order 10 20 ;

• for large values of N , f N has an amplitude of order 10 1 . For these same parameters, we plot the amplitude of f N (t) for different values of N , and t=5.0 fixed (Figure 3.8).

Suppose that we run the simulation of the whole system and find N cracks. For this same N , we find that the amplitude f N of the mode N is a lot bigger than for the other modes.

Since all the perturbation modes are present, we are going to consider that the one that dominates the system is the first one to be amplified of a factor S.

When using the library FEniCS, we find that max(ε)-min(ε) ≈ 10 -16 in the first iteration. We are going to consider that the perturbation becomes visible when they have the amplitude of 10 -3 . For this reason, we choose S=10 13 .

Description of the Spalling Tests

In this section, we will describe the experiments and their numerical modelling. The details of the experiments can be seen in Zinszner et al. [START_REF] Zinszner | Dynamic fragmentation of an alumina ceramic subjected to shockless spalling: An experimental and numerical study[END_REF].

Spalling Experiments

For each experiment, a specimen of height 10mm was placed in a pulse machine. The machine applied a compressive stress at one of the extremities of the specimen, while the other extremity was free of constraints. A compressive wave then propagated in the material and reached the free surface. The wave then returned, but now as a traction wave of same amplitude, causing the fracture of the specimen.

From each experiment, we obtained three informations:

1. the input stress as a function of time; This traction wave is then responsible for the fracture of the bar. We then compare the damaged zones obtained in the simulations to the regions where cracks appear in the experiments. The extremity x=L is free and vibrates. The velocity of this free surface is also measured and compared to numerical results.

For these simulations, we impose a damage criterion that allows fracture only in traction (σ>0).

The material used was the Al23 alumina. We have ρ = 3890 kg/m 3 , E=350 GPa, σ c = 202.8 MPa, L=10mm.

Using the standard AT1 model

We first run the simulations using the standard AT1 model. We recall that E(α)=E 0 (1 -α) 2 and w(α)=w 1 α.

The behaviour of a bar using this damage model is well-known. Using the analytic results for a quasi-static loading, we know that we have an homogeneous elastic phase and then sudden rupture when the stress reaches the critical value σ c . We can thus find

We can also find in the literature that K IC = 4 MPa • m 1/2 . Comparing it to the energy obtained in a quasi-static test, we have

We can then find the characteristic length = 0.000206m.

We can see in Figure 4.3 the damage profile. The region where α=1 is roughly the same as the regions where, in the experiments, we had a large crack density.

In Figure 4.4, we see the velocity of the free surface. We take E(α) = E 0 (1 -α) q for different values of q. Using the critical stress σ c , we can easily find w 1 = σ 2 c E q 2 . We can see that the value of q does not produce significant changes in the surface velocity. When comparing the experimental to the numerical results, we see that we have a good region for the cracks, but the free surface velocity is overestimated. In this work, we present possible modifications to the model in order to better approach the free surface velocity.

-Conclusion of the Chapter

Chapter 5 -Regularization of the Plastic Strain

The objective of this chapter is to propose models and study their onedimensional behaviour. The analysis performed here is far from extensive and a more rigorous work would be needed for each case.

Temperature-Plasticity Coupling

The objective of this section is to propose a one-dimensional model where the plastic strain does not localize in one single element for softening materials. The main idea is to use the dissipative properties of the heat equation in order to regularize the plastic strain.

We want to study a softening material such that the yield stress decreases when the temperature increases. We work with a bar subjected to a quasistatic loading. We will deduce the three equations that dictate the evolution problem:

• static equilibrium;

• evolution of the plastic strain by comparing the stress to the yield stress;

• evolution of the temperature, where the energy dissipated by the plastic process appears as a heat source.

We then conclude this section by discussing the difficulties found and why this model is not suitable.

Temperature-plasticity coupling

In this model, the evolution of the plastic strain is the only dissipative mechanism.

We consider a material with density density ρ 0 , internal energy e, speed v, temperature T , entropy s, heat flow q and volume heat source h.

We begin this section by stating the energy balance in its integral form:

for any volume Ω.

Locally, this is equivalent to

(5.2)

-Temperature-Plasticity Coupling

We use the dynamic equation (ρ 0 v = divσ + f ) to find ρ 0 ė = div(q) + ρ 0 h + σ ε.

(5.3)

We define the free energy ψ(ε, T, ε p ) = e -T s.

(5.4)

We recall that, using the second law of thermodynamics, we must have σ = ρ 0 ∂ψ ∂ε and s = -∂ψ ∂T . Therefore, the energy balance becomes ρ 0 ∂ψ ∂ε p εp + T ṡ = div(q) + ρ 0 h.

(5.5)

For this initial model, we suppose there is no dilation and that

(5.7)

We have the stress σ = E(ε -ε p ) (5.8)

and the equilibrium div σ = 0.

(5.9)

We can write the plasticity criterion as |σ| ≤ σ Y (T ).

(5.10)

By differentiating the free energy with respect to T and then differentiating it in time, we find ṡ = c 0 Ṫ . (5.11) Considering that there is no volume heat h and that the heat flow is q = k∇T , (5.12)

we find that ρ 0 c 0 T Ṫ = k∆T + σ Y (T ) εp .

(5.13)

Chapter 5 -Regularization of the Plastic Strain By dividing both characteristic lengths by a constant, one could expect that the damage and plastic strain profiles will also be divided by the same constant. This is not what we observe, however. We consider the case α = p =0.05. Figures 5.14 and 5.15 show that the behaviour of the system changes completely. The damage profile jumps instantly from zero to one and the plastic strain remains uniform, without any localization.

In order to better understand what happens we will fix the value of α and change the value of p . We will first study the support of the localized plastic zone. We know that the damage phase begins when p=0.75 and the plastic strain is always increasing. We are therefore going to calculate the size of the support of the plastic strain by measuring the interval where p>0.75.

We see in Figures 5.16 and 5.17 the size of the localized zone for different values of p , when α =0.05 and α =0.1, respectively.

We see in Figure 5.16 that there is a change in the behaviour when