Laura De Lorenzis

Claudia Comi

Jose Adachi

l'opportunité de rencontrer l'équipe R&D dont

Résumé

La défaillance d'une structure est souvent provoquée par la propagation de fissures dans le matériau initialement sain qui la compose. Une connaissance approfondie dans le domaine de la mécanique de la rupture est essentielle pour l'ingénieur. Il permet notamment de prévenir des mécanismes de fissurations garantissant l'intégrité des structures civiles, ou bien, de les développer comme par exemple dans l'industrie pétrolière. Du point de vue de la modélisation ces problèmes sont similaires, complexes et difficiles. Ainsi, il est fondamental de pouvoir prédire où et quand les fissures se propagent.

Ce travail de thèse se restreint à l'étude des fissures de type fragile et ductile dans les matériaux homogènes sous chargement quasi-statique. On adopte le point de vue macroscopique c'est-à-dire que la fissure est une réponse de la structure à une sollicitation excessive et est caractérisée par une surface de discontinuité du champ de déplacement. La théorie la plus communément admise pour modéliser les fissures est celle de Griffith. Elle prédit l'initiation de la fissure lorsque le taux de restitution d'énergie est égal à la ténacité du matériau le long d'un chemin préétabli. Ce type de critère requière d'évaluer la variation de l'énergie potentielle de la structure à l'équilibre pour un incrément de longueur de la fissure. Mais l'essence même de la théorie de Griffith est une compétition entre l'énergie de surface et l'énergie potentielle de la structure.

Cependant ce modèle n'est pas adapté pour des singularités d'entaille faible i.e. une entaille qui ne dégénère pas en pré-fissure. Pour pallier à ce défaut des critères de type contraintes critiques ont été développés pour des géométries régulières. Malheureusement ils ne peuvent prédire correctement l'initiation d'une fissure puisque la contrainte est infinie en fond d'entaille. Une seconde limitation de la théorie de Griffith est l'effet d'échelle. Pour illustrer ce propos, considérons une structure unitaire coupé par une fissure de longueur a. Le chargement critique de cette structure évolue en 1/ √ a, par conséquent le chargement admissible est infini lorsque la taille du défaut tend vers zéro.

Ceci n'a pas de sens physique et est en contradiction avec les expériences. Il est connu que cette limitation provient du manque de contrainte critique (ou longueur caractéristique) dans le modèle. Pour s'affranchir de ce défaut Dugdale et Barenblatt ont proposé dans leurs modèles de prendre en compte des contraintes cohésives sur les lèvres de la fissure afin d'éliminer la singularité de contraintes en fond d'entaille.

i Plus récemment, les modèles variationnels à champ de phase aussi connu sous le nom de modèles d'endommagements à gradient [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF] ont fait leurs apparition début des années 2000. Ces modèles permettent de s'affranchir des problèmes liés aux chemins de fissures et sont connus pour converger vers le modèle de Griffith lorsque le paramètre de régularisation tend vers 0. De plus les résultats numériques montrent qu'il est possible de faire nucléer une fissure sans singularité grâce à la présence d'une contrainte critique. Ces modèles à champ de phase pour la rupture sont-ils capables de surmonter les limitations du modèle de Griffith ?

Concernant les chemins de fissures, les modèles à champ de phase ont prouvé être redoutablement efficaces pour prédire les réseaux de fractures lors de chocs thermiques [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF]. Dans cette thèse, les résultats obtenus montrent que les modèles d'endommagement à gradient sont efficaces pour prédire la nucléation de fissure en mode I et de tenir compte de l'effet d'échelle. Naturellement ces modèles respectent le critère d'initiation de la théorie de Griffith et sont étendus à la fracturation hydraulique comme illustré dans le deuxième volet de cette thèse. Cependant ils ne peuvent rendre compte de la rupture de type ductile tel quel. Un couplage avec les modèles de plasticité parfaite est nécessaire afin d'obtenir des mécanismes de rupture ductile semblables à ceux observés pour les métaux.

Le manuscrit est organisé comme suit: Dans le premier chapitre, une large introduction est dédiée à l'approche variationnelle de la rupture en partant de Griffith vers une approche moderne des champs de phase en rappelant les principales propriétés. Le second chapitre étudie la nucléation de fissures dans des géométries pour lesquels il n'existe pas de solution exacte. Des entailles en U-et V-montrent que le chargement critique évolue continûment du critère en contrainte critique au critère de ténacité limite avec la singularité d'entaille. Le problème d'une cavité elliptique dans un domaine allongé ou infini est étudié. Le troisième chapitre se concentre autour de la fracturation hydraulique en prenant en compte l'influence d'un fluide parfait sur les lèvres de la fissure. Les résultats numériques montrent que la stimulation par injection de fluide dans d'un réseau de fissures parallèles et de même longueur conduit à la propagation d'une seule des fissures du réseau. Il s'avère que cette configuration respecte le principe de moindre énergie. Le quatrième chapitre se focalise uniquement sur le modèle de plasticité parfaite en partant de l'approche classique vers une l'approche variationnelle. Une implémentation numérique utilisant le principe de minimisation alternée de l'énergie est décrite et vérifiée dans un cas simple de Von Mises. Le dernier chapitre couple les modèles d'endommagement à gradient avec les modèles de plasticité parfaite. Les simulations numériques montrent qu'il est possible d'obtenir des fissures de type fragile ou ductile en variant un seul paramètre uniquement. De plus ces simulations capturent qualitativement le phénomène de nucléation et de propagation de fissures en suivant les bandes de cisaillement.

Introduction

Structural failure is commonly due to fractures propagation in a sound material. A better understanding of defect mechanics is fundamental for engineers to prevent cracks and preserve the integrity of civil buildings or to control them as desired in energy industry for instance. From the modeling point of view those problems are similar, complex and still facing many challenges. Common issues are determining when and where cracks will propagate.

In this work, the study is restricted to brittle and ductile fractures in homogeneous materials for rate-independent evolution problems in continuum mechanics. We adopt the macroscopic point of view, such that, the propagation of a macro fracture represents a response of the structure geometry subject to a loading. A fracture à la Griffith is a surface of discontinuity for the displacement field along which stress vanishes. In this widely used theory the fracture initiates along an a priori path when the energy release rate becomes critical, this limit is given by the material toughness. This criterion requires one to quantify the first derivative of potential energy with respect to the crack length for a structure at the equilibrium. Many years of investigations were focused on the notch tips to predict when the fracture initiates, resulting to a growing body of literature on computed stress intensity factors. Griffith is by essence a competition between the surface energy and the recoverable bulk energy. Indeed, a crack increment reduces the potential energy of the structure while it is compensated by the creation of a surface energy.

However such a fracture criterion is not appropriate to account for weak singularity i.e. a notch angle which does not degenerate into a crack. Conversely many criteria based on a critical stress are adapted for smooth domains, but fail near stress singularities. Indeed, a nucleation criterion based solely on pointwise maximum stress will be unable to handle with crack formation at the singularity point i.e. σ → ∞. A second limitation of Griffith's theory is the scale effects. To illustrate this, consider a unit structure size cut by a pre-fracture of length a. The critical loading evolves as ∼ 1/ √ a, consequently the maximum admissible loading is not bounded when the defect size decays. Again this is physically not possible and is inconsistent with experimental observations. It is well accepted that this discrepancy is due to the lack of a critical stress (or a critical length scale) in Griffith's theory. To overcome these aforementioned issues in Griffith's theory, Dugdale and Barenblatt pioneers of cohesive and ductile fractures theory proposed to kill the stress singularity at the tip by accounting of stresses on fracture lips.

iii Recently, many variational phase-field models [START_REF] Bourdin | The Variational Approach to Fracture[END_REF] are known to converge to a variational Griffith -like model in the vanishing limit of their regularization parameter. They were conceived to handle the issues of crack path. Furthermore, it has been observed that they can lead to numerical solution exhibiting crack nucleation without singularities. Naturally, these models raise some interesting questions: can Griffith limitations be overcome by those phase-field models?

Concening crack path, phase-field models have proved to be accurate to predict fracture propagation for thermal shocks [START_REF] Sicsic | Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF]. In this dissertation, numerical examples illustrate that Griffith limitations such as, nucleation and size effects can be overcome by the phase-field models referred to as gradient damage models in the Chapter 2. Naturally this models preserve Griffith's propagation criterion as shown in the extended models for hydraulic fracturing provided in Chapter 3. Of course Griffith's theory is unable to deal with ductile fractures, but in Chapter 5 we show that by coupling perfect plasticity with gradient damage models we are able to capture some of ductile fractures features, precisely the phenomenology of nucleation and propagation.

The dissertation is organized as follows: In Chapter 1, a large introduction of phasefield models to brittle fracture is exposed. We start from Griffith to the modern approach of phase-field models, and recall some of their properties. Chapter 2 studies crack nucleation in commonly encountered geometries for which closed-form solutions are not available. We use U-and V-notches to show that the nucleation load varies smoothly from that predicted by a strength criterion to that of a toughness criterion, when the strength of the stress concentration or singularity varies. We present validation and verifications of numerical simulations for both types of geometries. We consider the problem of an elliptic cavity in an infinite or elongated domain to show that variational phase-field models properly account for structural and material size effects. Chapter 3 focuses on fractures propagation in hydraulic fracturing, we extend the variational phase-field models to account for fluid pressure on the crack lips. We recover the closed form solution of a perfect fluid injected into a single fracture. For stability reason, in this example we control the total amount of injected fluid. Then we consider a network of parallel fractures stimulated. The numerical results show that only a single crack grows and this situation is always the best energy minimizer compared to a multi-fracking case where all fractures propagate. This loss of symmetry in the cracks patterns illustrates the variational structure and the global minimization principle of the phase-field model. A third example deals with fracture stability in a pressure driven laboratory test for rocks. The idea is to capture different stability regimes using linear elastic fracture mechanics to properly design the experiment. We test the phase-field models to capture fracture stability transition (from stable to unstable). Chapter 4 is concerned with the variational perfect plasticity models and its implementation and verification. We start by recalling main ingredients of the classic approach of perfect elasto-plasticity models and then recasting into the variational structure. Later the algorithm strategy is exposed with a verification example. The strength of the proposed algorithm is to solve perfect elasto-plasticity maiv terials by prescribing the yield surfaces without dealing with non differentiability issues. Chapter 5 studies ductile fractures, the proposed model couple gradient damage models with perfect plasticity independently exposed in Chapter 1 and 4. Numerical simulations show that transition from brittle to ductile fractures is recovered by changing only one parameter. Also the ductile fracture phenomenology, such as crack initiation at the center and propagation along shear bands are studied in plane strain specimens and round bars in three dimensions.

The main research contributions is in Chapter 2,3 and 5. My apologies to the reader perusing the whole dissertation which contains repetitive elements due to self consistency and independent construction of all chapters. v Chapter 1

Variational phase-field models of brittle fracture In Griffith's theory, a crack in brittle materials is a surface of discontinuity for the displacement field with vanishing stress along the fracture. Assuming an a priori known crack path, the fracture propagates when the first derivative of the potential energy with respect to the crack length at the equilibrium becomes critical. This limit called the fracture toughness is a material property. The genius of Griffith was to link the crack length to the surface energy, so the crack propagation condition becomes a competition between the surface energy and the recoverable bulk energy. By essence this criterion is variational and can be recast into a minimality principle. The idea of Francfort and Marigo in variational approach to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] is to keep Griffith's view and extend to any possible crack geometry and complex time evolutions. However cracks remain unknown and a special method needs to be crafted. The approach is to approximate the fracture by a damage field with a non zero thickness. In this region the material stiffness is deteriorated leading to decrease the sustainable stresses. This stress-softening material model is ill-posed mathematically [START_REF] Comi | On localisation in ductile-brittle materials under compressive loadings[END_REF] due to a missing term limiting the damage localization thickness size. Indeed, since the surface energy is proportional to the damage thickness size, we can construct a broken bar without paying any surface energy, i.e. by decaying the damaged area. To overcome this aforementioned issue, the idea is to regularize the surface energy. The adopted regularization takes its roots in Ambrosio and Tortorelli's [START_REF] Ambrosio | Existence theory for a new class of variational problems[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] functionals inspired by Mumford-Shah's work [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problem[END_REF] in image segmentation. Gradient damage models is closely related to Ambrosio and Tortorelli's functionals and have been adapted to brittle fracture. The introduction of a gradient damage term comes up with a regularized parameter. This parameter denoted is also called internal length and governs the damage thickness. Following Pham and Marigo [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Pham | From the onset of damage until the rupture: construction of the responses with damage localization for a general class of gradient damage models[END_REF], the damage evolution problem is built on three principles, the damage irreversibility, stability and balance of the total energy. The beauty of the model is that the unknown discrete crack evolution is approximated by a regularized functional evolution which is intimately related to Griffith by its variational structure and its asymptotic behavior.

This chapter is devoted to a large introduction of gradient damage models which 1

Chapter 1. Variational phase-field models of brittle fracture constitute a basis of numerical simulations performed in subsequent chapters. The presentation is largely inspired by previous works of Bourdin-Maurini-Marigo-Francfort and many others. In the sequel, section 1.1 starts with the Griffith point of view and recasts the fracture evolution into a variational problem. By relaxing the pre-supposed crack path constraint in Griffith's theory, the Francfort and Marigo's variational approach to fracture models is retrieved. We refer the reader to [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] for a complete exposition of the theory. Following the spirit of the variations principle, gradient damage models are introduced and constitute the basis of numerical simulations performed. Section 1.2 focuses on the application to a relevant one-dimensional problem which shows up multiple properties, such as, nucleation, critical admissible stress, size effects and optimal damage profile investigated previously by [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF]. To pass from a damage model to Griffith-like models, connections need to be highlighted, i.e letting the internal length to zero. Hence, section 1.3 is devoted to the Γ-convergence in one-dimensional setting, to show that gradient damage models behave asymptotically like Griffith. Finally, the implementation of such models is exposed in section 1.4 .

Gradient damage models

From Griffith model to its minimality principle

The Griffith model can be settled as follow, consider a perfectly brittle-elastic material with A the Hooke's law tensor and G c the critical energy release rate occupying a region Ω ⊂ R n in the reference configuration. The domain is partially cut by a fracture set Γ of length l, which grows along an a priori path Γ. Along the fracture, no cohesive effects or contact lips are considered here, thus, it stands for stress free on Γ(l). The sound region Ω \ Γ is subject to a time dependent boundary displacement ū(t) on a Dirichlet part of its boundary ∂ D Ω and time stress dependent g(t) = σ •ν on the remainder ∂ N Ω = ∂Ω\∂ D Ω, where ν denotes the appropriate normal vector. Also, for the sake of simplicity, body force is neglected. The infinitesimal total deformation e(u) is the symmetrical part of the spatial gradient of the displacement field u such that e(u) = ∇u + ∇ T u 2 .

In linear elasticity the free energy is a differentiable convex state function given by ψ e(u) = 1 2 Ae(u) : e(u) . Thereby, the stress-strain relation naturally follows σ = ∂ψ(e) ∂e = Ae(u).

By the quasi-static assumption made, the cracked solid is, at each time, in elastic equilibrium with the loads that it supports at that time. The problem is finding the unknown displacement u = u(t, l) for a given t and l = l(t) that satisfies the following constitutive equations,

1.1. Gradient damage models          div σ =0 in Ω \ Γ(l) u =ū(t) on ∂ D Ω \ Γ(l) σ • ν =g(t) on ∂ N Ω σ • ν =0 on Γ(l) (1.1)
At the time t and for l(t) let the kinematic field u(t, l) be at the equilibrium such that it solves (1.1). Hence, the potential energy can be computed and is composed of the elastic energy and the external work force, such that,

P(t, l) = Ω\Γ(l) 1 2
Ae(u) : e(u) dx -

∂ N Ω g(t) • u dH n-1
where dH n-1 denotes the Hausdorff n-1 -dimensional measure, i.e. its aggregate length in two dimensions or surface area in three dimensions. The evolution of the crack is given by Griffith's criterion:

Definition 1 (Crack evolution by Griffith's criterion)

i. Consider that the crack can only grow, this is the irreversibility condition, l(t) ≥ 0.

ii. The stability condition says that the energy release rate G is bounded from above by its critical value G c ,

G(t, l) = - ∂P(t, l) ∂l ≤ G c .
iii. The energy balance guarantee that the energy release rate is critical when the crack grows,

G(t, l) -G c l = 0
Griffith says in his paper [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], the "theorem of minimum potential energy" may be extended so as to of predicting the breaking loads of elastic solids, if account is taken of the increase of surface energy which occurs during the formation of cracks. Following Griffith, let us demonstrate that crack evolution criteria are optimality conditions of a total energy to minimize. Provided some regularity on P(t, l) and l(t), let formally the minimization problem be: for any loading time t such that the displacement u is at the Chapter 1. Variational phase-field models of brittle fracture equilibrium, find the crack length l which minimizes the total energy composed of the potential energy and the surface energy subject to irreversibility, min l≥l(t)

P(t, l) + G c l (1.2)
An optimal solution of the above constraint problem must satisfy the KKT1 conditions. A common methods consist in computing the Lagrangian, given by, L(t, l, λ) := P(t, l) + G c l + λ(l(t)l) (1.3) where λ denotes the Lagrange multiplier. Then, apply the necessary conditions, Substitute the Lagrange multiplier λ given by the stationarity into the dual feasibility and complementary slackness condition to recover the irreversibility, stability and energy balance of Griffith criterion.

Futhermore, let the crack length l and the displacement u be an internal variables of a variational problem. Note that the displacement does not depend on l anymore. Provided a smooth enough displacement field and evolution of t → l(t) to ensure that calculations make sense, the evolution problem can be written as a minimality principle, such as, Definition 2 (Fracture evolution by minimality principle) Find stable evolutions of l(t), u(t) satisfying at all t: i. Initial conditions l(t 0 ) = l 0 and u(t 0 , l 0 ) = u 0 ii. l(t), u(t) is a minimizer of the total energy,

E(t, l, u) = Ω\Γ(l) 1 2
Ae(u) : e(u)dx -

∂ N Ω g(t) • u dH n-1 + G c l (1.5) amongst all l ≥ l(t) and u ∈ C t := u ∈ H 1 (Ω \ Γ(l)) : u = ū(t) on ∂ D Ω \ Γ(l) .
1.1. Gradient damage models

iii. The energy balance,

E(t, l, u) = E(t 0 , l 0 , u 0 ) + t t 0 ∂ D Ω (σ • ν) • u dH n-1 - ∂ N Ω ġ(t) • u dH n-1 ds (1.6)
One observes that stability and irreversibility have been substituted by minimality, and the energy balance takes a variational form. To justify this choice, we show first that irreversibility, stability and kinematic equilibrium are equivalent to the first order optimality conditions of E(t, l, u) for u and l separately. Then, followed by the equivalence of the energy balance adopted in the evolution by minimality principle and within Griffith criterion.

Proof. For a fixed l, u is a local minimizer of E(t, l, u), if for all v ∈ H 1 0 (Ω \ Γ(l)), for some h > 0 small enough, such that u + hv ∈ C t ,

E(t, l, u + hv) = E(t, l, u) + hE (t, l, u) • v + o(h) ≥ E(t, l, u) (1.7)
thus,

E (t, l, u) • v ≥ 0 (1.8)
where E (t, l, u) denotes the first Gateaux derivative of E at u in the direction v. By standard arguments of calculus of variations, one obtains,

E (t, l, u) • v = Ω\Γ(l) 1 2 Ae(u) : e(v)dx - ∂ N Ω g(t) • v dH n-1 (1.9)
Integrating the term in e(v) by parts over Ω \ Γ(l), and considering both faces of Γ(l) with opposites normals, one gets,

E (t, l, u) • v = - Ω\Γ(l) div Ae(u) • v dx + ∂Ω Ae(u) • ν • v dH n-1 - Γ(l) Ae(u) • ν • v dH n-1 - ∂ N Ω g(t) • v dH n-1 (1.
E (t, l, u) • v = - Ω\Γ(l) div Ae(u) • v dx + ∂ N Ω Ae(u) • ν -g(t) • v dH n-1 - Γ(l) Ae(u) • ν • v dH n-1
(1.11)

Chapter 1. Variational phase-field models of brittle fracture

Taking v = -v ∈ H 1 0 (Ω\Γ(l))
, the optimality condition leads to E (t, l, u)•v = 0. Formally by a localization argument taking v such that it is concentrated around boundary and zero almost everywhere, we obtain that all integrals must vanish for any v. Since the stress-strain relation is given by σ = Ae(u), we recover the equilibrium constitutive equations,

         div Ae(u) = 0 in Ω \ Γ(l) u = ū(t) on ∂ D Ω \ Γ(l) Ae(u) • ν = g(t) on ∂ N Ω Ae(u) • ν = 0 on Γ(l) (1.12)
Now consider u is given. For any l > 0 for some h > 0 small enough, such that l + h l ≥ l(t), the derivative of E(t, l, u) at l in the direction l is,

E (t, l, u) • l ≥ 0 ∂P(t, l, u) ∂l + G c ≥ 0 (1.13)
this becomes an equality, G(t, l, u) = G c when the fracture propagates.

To complete the equivalence between minimality evolution principle and Griffith, let us verify the energy balance. Provided a smooth evolution of l, the time derivative of the right hand side equation (1.6) is,

dE(t, l, u) dt = ∂ D Ω (σ • ν) • u dH n-1 - ∂ N Ω ġ(t) • u dH n-1 (1.14)
and the explicit left hand side,

dE(t, l, u) dt = E (t, l, u) • u + E (t, l, u) • l - ∂ N Ω ġ(t) • u dH n-1 . (1.15)
The Gateaux derivative with respect to u have been calculated above, so E (t, l, u) • u stands for,

E (t, l, u) • u = - Ω\Γ div(Ae(u)) • u dx + ∂ D Ω Ae(u) • ν • u dH n-1 + ∂ N Ω Ae(u) • ν • u dH n-1 - ∂ N Ω g(t) • u dH n-1 .
(1.16) Since u respects the equilibrium and the admissibility u = u on ∂ D Ω, all kinematic contributions to the elastic body vanish and the energy balance condition becomes,

E (t, l, u) • l = 0 ⇔ ∂P ∂l + G c l = 0 (1.17)

Gradient damage models

At this stage minimality principle is equivalent to Griffith criterion for smooth evolution of l(t). Let's give a graphical interpretation of that. Consider a domain partially cut by a pre-fracture of length l 0 subject to a monotonic increasing displacement load, such that, ū(t) = tū on ∂ D Ω and stress free on the remainder boundary part. Hence, the elastic energy is ψ e(tu) = t 2 2 Ae(u) : e(u) and the irreversibility is l ≥ l 0 . The fracture stability is given by

t 2 ∂P (1, l) ∂l + G c ≥ 0
and for any loading t > 0, the energy release rate for a unit loading is bounded by

G(1, l) ≤ G c /t 2 .
Forbidden region The fracture evolution is smooth if G(1, l) is strictly decreasing in l, i.e. P(1, l) is strictly convex as illustrated on the Figure 1.1(left). Thus, stationarity and local minimality are equivalent. Let's imagine that material properties are not constant in the structure, simply consider the Young's modulus varying in the structure such that G(1, l) has a concave part, see Figure 1.1(right). Since G(1, l) is a deceasing function, the fracture grows smoothly by local minimality argument until being stuck in the local well for any loadings which is physically inconsistent. Conversely, considering global minimization allows up to a loading point, the nucleation of a crack in the material, leading to a jump of the fracture evolution.

Extension to Francfort-Marigo's model

In the previous analysis, the minimality principle adopted was a local minimization argument because it considers small perturbations of the energy. This requires a topology, which includes a concept of distance defining small transformations, whereas for global minimization principle it is topology-independent. Without going too deeply into details, arguments in favor of global minimizers are described below. Griffith's theory does not Chapter 1. Variational phase-field models of brittle fracture hold for a domain with a weak singularity. By weak singularity, we consider any free stress acute angle that does not degenerate into a crack (as opposed to strong singularity). For this problem, by using local minimization, stationary points lead to the elastic solution. The reason for this are that the concept of energy release rate is not defined for a weak singularity and there is no sustainable stress limit over which the crack initiates. Hence, to overcome the discrepancy due to the lack of a critical stress in Griffith's theory, double criterion have been developed to predict fracture initiation in notched specimen, more details are provided in Chapter 2. Conversely, global minimization principle has a finite admissible stress allowing cracks nucleation, thus cracks can jump from a state to another, passing through energy barriers. For physical reasons, one can blame global minimizers to not enforce continuity of the displacement and damage field with respect to time. Nevertheless, it provides a framework in order to derive the fracture model as a limit of the variational damage evolution presented in section 1.3. This is quite technical but global minimizers from the damage model converge in the sens of Γ-convergence to global minimizers of the fracture model. Finally, under the assumptions of a pre-existing fracture and strict convexity of the potential energy, global or local minimization are equivalent and follow Griffith. In order to obtain the extended model of Francfort-Marigo variational approach to fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Marigo | Initiation of cracks in griffith's theory: An argument of continuity in favor of global minimization[END_REF][START_REF] Bourdin | Variational Models and Methods in Solid and Fluid Mechanics, chapter Fracture[END_REF] one has to keep the rate independent variational principle and the Griffith fracture energy, relax the constrain on the pre-supposed crack path by extending to all possible crack geometries Γ and consider the global minimization of the following total energy

E(u, Γ) := Ω\Γ 1 2
Ae(u) : e(u) dx -

∂ N Ω g(t) • u dH n-1 + G c H n-1 (Γ) (1.18)
associated to cracks evolution problem given by, Definition 3 (Crack evolution by global minimizers) u(t), Γ t satisfies the variational evolution associated to the energy E(u, Γ) if the following three conditions hold:

i. t → Γ t is increasing in time, i.e Γ t ⊇ Γ s for all t 0 ≤ s ≤ t ≤ T . ii. for any configuration (v, Γ) such that v = g(t) on ∂ D Ω \ Γ t and Γ ⊇ Γ t , E v, Γ ≥ E u(t), Γ t (1.19)
iii. for all t,

E u(t), Γ t = E u(t 0 ), Γ t 0 + t t 0 ∂ D Ω (σ • ν) • u(t) dH n-1 - ∂ N Ω ġ(t) • u dH n-1 ds (1.20)

Gradient damage models

It is convenient, to define the weak energy by extending the set of admissibility function to an appropriate space allowing discontinuous displacement field, but preserving "good" properties.

SBD(Ω) = u ∈ SBV (Ω); Du = ∇u + (u + -u -) • ν dH n-1 J(u) (1.21)
where, Du denotes the distributional derivative, J(u) is the jump set of u. Following De Giorgi in [START_REF] De Giorgi | Existence theorem for a minimum problem with free discontinuity set[END_REF], the minimization problem is reformulated in a weak energy form functional of SBV , such as,

min u∈SBV (Ω) Ω 1 2 Ae(u) : e(u) dx - ∂ N Ω g(t) • u dH n-1 + G c H n-1 J(u) (1.22)
For existence of solution in the discrete time evolution and time continuous refer to [START_REF] Francfort | Existence and convergence for quasi-static evolution in brittle fracture[END_REF][START_REF] Babadjian | Existence of strong solutions for quasi-static evolution in brittle fracture[END_REF]. The weak energy formulation will be recalled in section 1.3 for the Γ-convergence in one dimension.

Gradient damage models to brittle fracture

Because the crack path remains unknown a special method needs to be crafted. The approach is to consider damage as an approximation of the fracture with a finite thickness where material properties are modulated continuously. Hence, let the damage α being an internal variable which evolves between two extreme states, up to a rescaling α can be bounded between 0 and 1, where α = 0 is the sound state material and α = 1 refers to the broken part. Intermediate values of the damage can be seen as "micro cracking", a partial disaggregation of the Young's modulus. A possible choice is to let the damage variable α making an isotropic deterioration of the Hooke's tensor, i.e. a(α)A where a(α) is a stiffness function. Naturally the recoverable energy density becomes, ψ(α, e) = 1 2 a(α)Ae(u) : e(u), with the elementary property that ψ(α, e) is monotonically decreasing in α for any fixed u. The difficulty lies in the choice of a correct energy dissipation functional. At this stage of the presentation a choice would be to continue by following Marigo-Pham [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Pham | From the onset of damage until the rupture: construction of the responses with damage localization for a general class of gradient damage models[END_REF][START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF][START_REF] Pham | Stability of homogeneous states with gradient damage models: Size effects and shape effects in the three-dimensional setting[END_REF] for a full and self consistent construction of the model. Their main steps are, assume a dissipation potential k(α), apply the Drucker-Ilushin postulate, then, introduce a gradient damage term to get a potential dissipation of the form k(α, ∇α). Instead, we will continue by following the historical ideas which arose from the image processing field with Mumford-Shah [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problem[END_REF] where continuous functional was proposed to find the contour of the image in the picture by taking into account strong variations of pixels intensity across boundaries. Later, Ambrosio-Tortorelli [START_REF] Ambrosio | Existence theory for a new class of variational problems[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] proposed the following functional which constitute main ingredients of the regularized damage models,

Ω 1 2 (1 -α) 2 |∇u| 2 dx + Ω α 2 + 2 |∇α| 2 dx
where > 0 is a regularized parameter called internal length. One can recognize the second term as the dissipation potential composed of two parts, a local term depending only on the damage state and a gradient damage term which penalizes sharp localization of the damage. The regularized parameter came up with the presence of the gradient damage term which has a dimension of the length. Following [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF][START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF], we define the regularized total energy of the gradient damage model for a variety of local dissipations and stiffness functions denoted w(α) and a(α), not only w(α) = α 2 and a(α

) = (1 -α) 2 by E (u, α) = Ω 1 2 a(α)Ae(u) : e(u) dx - ∂ N Ω g(t) • u dH n-1 + G c 4c w Ω w(α) + |∇α| 2 dx (1.23)
where G c is the critical energy release rate, c w = 1 0 w(α)dα with w(α) and a(α) following some elementary properties.

1. The local dissipation potential w(α) is strictly monotonically increasing in α. For a sound material no dissipation occurs hence w(0) = 0, for a broken material the dissipation must be finite, and up to a rescaling we have w(1) = 1.

2. The elastic energy is monotonically decreasing in α for any fixed u. An undamaged material should conserve its elasticity property and no elastic energy can be stored in a fully damaged material such that, the stiffness function a(α) is a decreasing function with a(0) = 1 and a(1) = 0.

3.

For numerical optimization reasons one can assume that a(α) and w(α) are continuous and convex.

A large variety of models with different material responses can be constructed just by choosing different functions for a(α) and w(α). A non exhaustive list of functions used in the literature is provided in Table 1.1. Despite many models used, we will mainly focus on AT 1 and sometimes refers to AT 2 for numerical simulations. Now let us focus on the damage evolution of E (u, α) defined in (1.23). First, remark that to get a finite energy, the gradient damage is in L 2 (Ω) space. Consequently, the trace can be defined at the boundary, so, damage values can be prescribed. Accordingly let the set of admissible displacements and admissible damage fields C t and D, equipped with their natural H 1 norm,

C t = u ∈ H 1 (Ω) : u = ū(t) on ∂ D Ω , D = α ∈ H 1 (Ω) : 0 ≤ α ≤ 1, ∀x ∈ Ω .
The evolution problem is formally similar to one defined in Definition 2 and reads as,

1.1. Gradient damage models Name a(α) w(α) AT 2 (1 -α) 2 α 2 AT 1 (1 -α) 2 α LS k 1 -w(α) 1 + (c 1 -1)w(α) 1 -(1 -α) 2 KKL 4(1 -α) 3 -3(1 -α) 4 α 2 (1 -α) 2 /4 Bor c 1 (1 -α) 3 -(1 -α) 2 + 3(1 -α) 2 -2(1 -α) 3 α 2 SKBN (1 -c 1 ) 1 -exp (-c 2 (1 -α) c 3 ) 1 -exp (-c 2 ) α
Table 1.1: Variety of possible damage models, where c 1 , c 2 , c 3 are constants. AT 2 introduced by Ambrosio Tortorelli and used by Bourdin [START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF], AT 1 model initially introduced by Pham-Amor [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF], LS k in Alessi-Marigo [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF], KKL for Karma-Kessler-Levine used in dynamics [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF], Bor for Borden in [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], SKBN for Sargadoa-Keilegavlena-Berrea-Nordbottena in [START_REF] Sargado | High-accuracy phase-field models for brittle fracture based on a new family of degradation functions[END_REF].

Definition 4 (Damage evolution by minimality principle) For all t find (u, α) ∈ (C t , D) that satisfies the damage variational evolution:

i. Initial condition α t 0 = α 0 and u t 0 = u 0 ii. (u, α) is a minimizer of the total energy, E (u, α)

E (u, α) = Ω 1 2 a(α)Ae(u) : e(u) dx - ∂ N Ω g(t) • u dH n-1 + G c 4c w Ω w(α) + |∇α| 2 dx (1.24)
amongst all α ≥ α(t)

iii. Energy balance,

E (u t , α t ) = E (u 0 , α 0 ) + t t 0 ∂ D Ω (σ • ν) • u dH n-1 - ∂ N Ω ġ(t) • u dH n-1 ds. (1.25)
This damage evolution is written in a weak form in order to obtain the damage criterion in a strong formulation, we have to explicit the first order necessary optimality conditions of the constraint minimization of E for (u, α) given by,

E (u, α)(v, β) ≥ 0 ∀(v, β) ∈ H 1 0 (Ω) × D (1.26)
Chapter 1. Variational phase-field models of brittle fracture

Using calculus of variation argument, one gets,

E (u, α)(v, β) = Ω a(α)Ae(u) : e(v) dx - ∂ N Ω (g(t) • ν) • v dH n-1 + Ω 1 2 a (α)Ae(u) : e(u)β dx + G c 4c w Ω w (α) β + 2 ∇α • ∇β dx. (1.27)
Integrating by parts the first term in e(v) and the last term in ∇α • ∇β, the expression leads to,

E (u, α)(v, β) = - Ω div a(α)Ae(u) • v dx + ∂ N Ω [a(α)Ae(u) -g(t)] • ν • v dH n-1 + Ω 1 2 a (α)Ae(u) : e(u) + G c 4c w w (α) -2 ∆α β dx + G c 4c w ∂Ω 2 ∇α • ν β dH n-1 .
(1.28)

This holds for all β ≥ 0 and for all v ∈ H 1 0 (Ω), thus, one can take β = 0 and v = -v. Necessary, the first two integrals are equal to zero. Again, we recover the kinematic equilibrium with the provided boundary condition since σ = a(α)Ae(u),

     div a(α)Ae(u) = 0 in Ω a(α)Ae(u) = g(t) on ∂ N Ω u = ū(t) on ∂ D Ω (1.29)
The damage criteria and its associated boundary conditions arise for any β ≥ 0 and by taking v = 0 in (1.28), we obtain that the third and fourth integrals are non negative.

     1 2 a (α)Ae(u) : e(u) + G c 4c w w (α) -2 ∆α ≥ 0 in Ω ∇α • ν ≥ 0 on ∂Ω (1.30)
The damage satisfies criticality when (1.30) becomes an equality.

Before continuing with the energy balance expression, let us focus a moment on the damage criterion. Notice that it is composed of an homogeneous part depending in w (α) and a localized contribution in ∆α. Assume the structure being at an homogeneous damage state, such that α is constant everywhere, hence the laplacian damage term vanishes. In that case, the elastic domain in a strain space is given by,

1.1. Gradient damage models Ae(u) : e(u) ≤ G c 2c w w (α) -a (α) (1.31)
and in stress space, by, .32) this last expression requires to be bounded such that the structure has a maximum admissible stress,

A -1 σ : σ ≤ G c 2c w w (α)a(α) 2 -a (α) (1
max α w (α) c (α) < C (1.33)
where c(α) = 1/a(α) is the compliance function.

If α → w (α)/c (α) is increasing the material response will be strain-hardening. For a decreasing function it is a stress-softening behavior. This leads to,

w (α)a (α) > w (α)a (α) (Strain-hardening) w (α)c (α) < w (α)c (α) (Stress-softening) (1.34)
Those conditions restrict proper choice for w(α) and a(α).

Let us turn our attention back to find the strong formulation of the problem using the energy balance. Assuming a smooth evolution of damage in time and space, the time derivative of the energy is given by,

dE (u, α) dt = E (u, α)( u, α) - ∂ N Ω ( ġ(t) • ν) dH n-1 (1.35)
The first term has already been calculated by replacing (v, β) with ( u, α) in (1.27), so that,

dE (u, α) dt = - Ω div a(α)Ae(u) • u dx + ∂ N Ω [a(α)Ae(u) -g(t)] • ν • u dH n-1 + ∂ D Ω a(α)Ae(u) • ν • u dH n-1 - ∂ N Ω ġ(t) • ν • u dH n-1 + Ω 1 2 a (α)Ae(u) : e(u) + G c 4c w w (α) -2 ∆α α dx + G c 4c w ∂Ω 2 ∇α • ν α dH n-1
(1.36)

Chapter 1. Variational phase-field models of brittle fracture

The first line vanishes with the equilibrium and boundary conditions, the second line is equal to the right hand side of the energy balance definition (1.25). Since the irreversibility α ≥ 0 and the damage criterion (1.30) hold, the integral is non negative, therefore the energy balance condition gives,

     1 2 a (α)Ae(u) : e(u) + G c 4c w w (α) -2 ∆α • α = 0 in Ω (∇α • ν) • α = 0 on ∂Ω (1.37)
Notice that the first condition in (1.37) is similar to the energy balance of Griffith, in the sense that the damage criterion is satisfied when damage evolves. Finally, the evolution problem is given by the damage criterion (1.30), the energy balance (1.37) and the kinematic admissibility (1.29).

The next section is devoted to the construction of the optimal damage profile by applying the damage criterion to a one-dimensional traction bar problem for a given . Then, defined the critical energy release rate as the energy required to break a bar and to create an optimal damage profile.

Application to a bar in traction

The one-dimension problem

The aim of this section is to apply the gradient damage model to a one-dimensional bar in traction. Relevant results are obtained with this example such as, the role of critical admissible stress, the process of damage nucleation due to stress-softening, the creation of an optimal damage profile for a given and the role of gradient damage terms which ban spacial jumps of the damage.

In the sequel, we follow Pham-Marigo [START_REF] Pham | Construction et analyse de modèles d'endommagement à gradient[END_REF][START_REF] Pham | From the onset of damage until the rupture: construction of the responses with damage localization for a general class of gradient damage models[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF]] by considering a one-dimensional evolution problem of a homogeneous bar of length 2L stretched by a time controlled displacement at boundaries and no damage value is prescribed at the extremities, such that, the admissible displacement and damage sets are respectively,

C t := {u : u(-L) = -tL, u(L) = tL}, D := {α : 0 ≤ α ≤ 1 in [0, L]} (1.38)
with the initial condition u 0 (x) = 0 and α 0 (x) = 0. Since no external force is applied, the total energy of the bar is given by,

E (u, α) = L -L 1 2 a(α)Eu 2 dx + G c 4c w L -L w(α) + |α | 2 dx (1.39)
where E is the Young's modulus, > 0 and (•) = ∂(•)/∂x. For convenience, let the compliance being the inverse of the stiffness such that c(α) = a -1 (α). Assume that α is at least continuously differentiable, but a special treatment would be required for α = 1

1.2. Application to a bar in traction which is out of the scope in this example. The pair (u t , α t ) ∈ C t × D is a solution of the evolution problem if the following conditions holds:

1. The equilibrium,

σ t (x) = 0, σ t (x) = a(α t (x))Eu t (x), u t (-L) = -tL and u t (L) = tL
The stress is constant along the bar. Hence it is only a function of time, such that,

2tLE = σ t L -L c α t (x) dx (1.40)
Once the damage field is known. The equation (1.40) gives the stress-displacement response.

2. The irreversibility, αt (x) ≥ 0 (1.41)

3. The damage criterion in the bulk,

- c (α t (x)) 2E σ 2 t + G c 4c w w (α t (x)) -2 α t (x) ≥ 0 (1.42)
4. The energy balance in the bulk,

- c (α t (x)) 2E σ 2 t + G c 4c w w (α t (x)) -2 α t (x) αt (x) = 0 (1.43)
5. The damage criterion at the boundary, α t (-L) ≥ 0 and α t (L) ≤ 0 (1.44)

6. The energy balance at the boundary,

α t (±L) αt (±L) = 0 (1.45)
For smooth or brutal damage evolutions the first order stability enforce α t (±L) = 0 to respect E (u, α) = 0. Thus the damage boundary condition is replaced by α t (±L) = 0 when damage evolves. All equations are settled to solve the evolution problem. Subsequently, we study a uniform damage in the bar and then focus on the localized damage solution.

The homogeneous damage profile

Consider a case of a uniform damage in the bar α t (x) = α t , which is called the homogeneous solution. We will see that the damage response depends on the evolution of α → w (α)/c (α), i.e. for stress-hardening (increasing function) the damage evolves uniformly in the bar, and localizes for stress-softening configuration. Now suppose that the damage does not evolve and remains equal to its initial value, α t = α 0 = 0. Then using the damage criterion in the bulk (1.42) the admissible stress must satisfy,

σ 2 t ≤ 2EG c 4c w w (0) c (0) (1.46)
and the response remains elastic until the loading time t e , such that,

t 2 ≤ - G c 2E c w w (0) a (0) = t 2 e (1.47)
Suppose the damage evolves uniformly belongs to the bar, using the energy balance (1.43) and the damage criterion (1.42) we have,

σ 2 t ≤ 2EG c 4c w w (α t ) c (α t ) , σ 2 t - 2EG c 4c w w (α t ) c (α t ) αt = 0 (1.48)
The homogeneous damage evolution is possible only if α t → w (α)/c (α) is growing, this is the stress-hardening condition. Since αt > 0, the evolution of the stress is given by,

σ 2 t = 2EG c 4c w w (α t ) c (α t ) ≤ max 0<α<1 2EG c 4c w w (α t ) c (α t ) = σ 2 c (1.49)
where σ c is the maximum admissible stress for the homogeneous solution. One can define the maximum damage state α c obtained when σ t = σ c . This stage is stable until the loading time t c ,

t 2 ≤ - G c 2E c w w (α c ) a (α c ) = t 2 c (1.50)
Since w (α)/c (α) is bounded and > 0, a fundamental property of gradient damage model is there exists a maximum value of the stress called critical stress, which allows crack to nucleate using the minimality principle.

The localized damage profile

The homogeneous solution is no longer stable if the damage α t → w (α)/c (α) is decreasing after α c . To prove that, consider any damage state such that, α t (x) > α c and the stress-softening property, leading to, 1.2. Application to a bar in traction

0 ≤ 2EG c 4c w w (α t (x)) c (α t (x)) ≤ 2EG c 4c w w (α c ) c (α c ) = σ 2 c (1.51)
By integrating the damage criterion (1.42) over (-L, L) and using (1.44), we have,

σ 2 t 2E L -L c (α t (x)) dx ≤ G c 4c w L -L w (α t (x)) dx + 2 α t (L) -α t (-L) ≤ G c 4c w L -L w (α t (x)) dx (1.52)
then, put (1.51) into (1.52) to conclude that σ t ≤ σ c and use the equilibrium (1.40) to obtain σ t ≥ 0. Therefore using (1.52) we get that α t (x) ≥ 0, consequently the damage is no longer uniform when stress decreases 0 ≤ σ t ≤ σ c . Assume α t (x) is monotonic over (-L, x 0 ) with α t (-L) = α c and the damage is maximum at x 0 , such that, α t (x 0 ) = max x α t (x) > α c . Multiplying the equation (1.42) by α t (x), an integrating over [-L, x) for x < x 0 we get,

α 2 t (x) = - 2c w σ 2 t EG c c(α t (x)) -c(α c ) + w(α t (x)) -w(α c ) (1.53)
Plugging this above equation into the total energy restricted to the (-L, x 0 ) part,

E (u t (x), α t (x)) (-L,x 0 ) = x 0 -L σ 2 t 2a(α t (x))E dx + G c 4c w x 0 -L w(α t (x)) + α 2 t (x) dx = x 0 -L σ 2 t 2a(α c )E dx + G c 4c w x 0 -L 2w(α) -w(α c ) dx (1.54)
Note that the energy does not depend on α anymore, we just have two terms: the elastic energy and the surface energy which depends on state variation of w(α).

The structure is broken when the damage is fully localized α(x 0 ) = 1. From the equilibrium (1.40), the ratio stress over stiffness function is bounded such that |σ t c(α)| < C, thus, |σ 2 t c(1)| → 0 and (1.53) becomes,

α 2 t (x) = w(α t (x)) -w(α c ) , ∀x ∈ (-L, x 0 )
Remark that, the derivative of the damage and u across the point x 0 where α = 1 is finite. By letting the variable β = α t (x), the total energy of the partial bar (-L, x 0 ) is Chapter 1. Variational phase-field models of brittle fracture

E (u t (x), α t (x)) (-L,x 0 ) = lim x→x 0 G c 4c w x -L 2w(α t (x)) -w(α c ) dx = lim β→1 G c 4c w β αc 2w(α) -w(α c ) β dβ = lim β→1 G c 4c w β αc 2w(β) -w(α c ) w(β) -w(α c ) dβ = lim β→1 G c 4c w β αc 2 w(β) -w(α c ) + w(α c ) w(β) -w(α c ) dβ = G c 2c w k(α c ) (1.55) with, k(α c ) := 1 αc w(β) -w(α c ) dβ + w(α c ) D 4 ,
where D is the damage profile size between the homogeneous and fully localized state, given by,

D = L -L dx α (x) = 1 αc 2 w(β) -w(α c ) dβ. (1.56)
Note that the right side of the bar (x 0 , L) contribute to the exact same total energy than the left one (-L, x 0 ). Different damage response is observed depending on the choice of w(α) and a(α). The model AT 1 for instance has an elastic part, thus α c = 0 and the energy release during the breaking process of a 1d bar is equal to G c . Models with an homogeneous response before localization, AT 2 for example, overshoot G c due to the homogeneous damage profile. A way to overcome this issue, is to consider that partial damage do not contribute to the dissipation energy, it can be relaxed after localization by removing the irreversibility. Another way is to reevaluate c w such as, c w = k(α c ).

Limit of the damage energy

From inception to completion gradient damage models follows the variational structure of Francfort-Marigo's [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] approach seen as an extension of Griffith, but connections between both need to be highlighted. Passing from damage to fracture, i.e. letting → 0 requires ingredients adapted from Ambrosio Tortorelli [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] on convergence of global minimizers of the total energy. A framework to study connections between damage and fracture variational models is that of Γ-convergence which we briefly introduce below. We refer the reader to [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF][START_REF] Braides | Gamma-convergence for Beginners[END_REF][START_REF] Dal | An introduction to Γ-convergence[END_REF] for a complete exposition of the underlying theory.

Limit of the damage energy

In the sequel, we restrict the study to a 1d case structure of interval Ω ⊂ R whose size is large compare to the internal length and with a unit Young's modulus. We prescribe a boundary displacement ū on a part ∂ D Ω and stress free on the remaining part ∂ N Ω := ∂Ω \ ∂ D Ω. We set aside the issue of damage boundary conditions for now and we define the weak fracture energy,

E(u, α, Ω) = F(u, Ω) if u ∈ SBV (Ω) +∞ otherwise (1.57)
and

F(u, Ω) := 1 2 Ω (u ) 2 dx + G c #(J(u)) (1.58)
where #(J(u)) denotes the cardinality of jumps in the set of u. Derived from E(u, α, Ω) its associated regularized fracture energy is,

E (u, α, Ω) = F (u, α, Ω) if u ∈ W 1,2 (Ω), α ∈ W 1,2 (Ω; [0, 1]) +∞ otherwise (1.59)
and

F (u, α, Ω) := 1 2 Ω a(α)(u ) 2 dx + G c 4c w Ω w(α) + α 2 dx (1.60)
To prove that up to a subsequence minimizers for E converge to global minimizers of E we need the fundamental theorem of the Γ-convergence given in the Appendix A.

We first show the compactness of the sequence of minimizers of E , then the Γconvergence of E to E. Before we begin, let the truncation and optimal damage profile lemma be, Lemma 1 Let u (resp. (u, α)) be a kinematically admissible global minimizer of F (resp.

F ). Then u L ∞ (Ω) ≤ ū L ∞ (Ω) Proof. Let M = ū L ∞ , and u * = inf {sup{-M, u}, M }. Then F(u * ) ≤ F(u) with equality if u = u * .
Lemma 2 Let α be the optimal profile of

S (α ) := I w(α ) + (α ) 2 dx
where I ⊂ R, then S (α ) = 4c w .

Proof. In order to construct α we solve the optimal profile problem: Let γ be the solution of the following problem: find γ ∈ C 1 [-δ, x 0 ) such that γ(-δ) = 0 and lim x→x 0 γ(x) = ϑ, and which is a minimum for the function,

F (γ) = x 0 -δ f (γ(x), γ (x), x)dx (1.61)
where

f (γ(x), γ (x), x) := w(γ(x)) + γ 2 (x) (1.62)
Note that the first derivative of f is continuous. We will apply the first necessary optimality condition to solve the optimization problem described above, if γ is an extremum of F , then it satisfies the Euler-Lagrange equation,

2γ = w (γ) 2 and γ (-δ) = 0 (1.63)
Note that w (γ) ≥ 0 implies γ convex, thus γ is monotonic in [-δ, x 0 ). Multiplying by γ and integrating form -δ to x, we obtain,

γ 2 (x) -γ 2 (-δ) = w(γ(x)) -w(γ(-δ)) 2 (1.64)
Since γ (-δ) = 0 and w(γ(-δ)) = 0, one gets,

γ (x) = w(γ(x)) 2 (1.65) Let us define, α (x) = γ (|x -x 0 |) then, α (x) := γ (|x -x 0 |) if |x -x 0 | ≤ δ 0 otherwise (1.66)
Note that α is continuous at x 0 and values ϑ, we have that,

S (α ) = I w(α ) + (α ) 2 dx = 2 x 0 -δ w(γ ) + (γ ) 2 dx (1.67)
Plug (1.65) into the last integral term, and change the variables β = γ (x), it turns into

S (α ) = 2 x 0 -δ w(γ ) + (γ ) 2 dx = 2 γ(x 0 ) γ(-δ) w(β) β dβ = 4 ϑ 0 w(β) dβ (1.68)
The fully damage profile is obtained once ϑ → 1, we get, This will be usefull for the recovery sequence in higer dimensions.

S (α ) = lim

Compactness

Theorem 1 Let (x) :=

x 0 w(s)ds, and assume that there exists C > 0 such that 1 -(s) ≤ C a(s) for any 0 ≤ s ≤ 1. Let (u , α ) be a kinematic admissible global minimizer of E . Then, there exists a subsequence (still denoted by (u , α ) ), and a function u ∈ SBV (Ω) such that u → u in L 2 (Ω) and α → 0 a.e. in Ω as → 0 Proof. Note that the technical hypothesis is probably not optimal but sufficient to account for the AT 1 and AT 2 functionals. Testing α = 0 and an arbitrary kinematically admissible displacement field ũ, we get that,

E (u , α ) ≤ E (ũ, 0) ≤ 1 2 Ω |ũ | 2 dx ≤ C (1.69)
So that E (u , α ) is uniformly bounded by some C > 0. Also, this implies that w(α ) → 0 almost everywhere in Ω, and from properties of w, that α → 0 almost everywhere in Ω.Using the inequality a 2 + b 2 ≥ 2|ab| on the surface energy part, we have that,

Ω 2 w(α )|α | dx ≤ Ω w(α ) + (α ) 2 dx ≤ C (1.70)
In order to obtain the compactness of the sequence u , let v := (1 -(α )) u and using the truncation Lemma 1, v is uniformly bounded in L ∞ (Ω). Then,

v = (1 -(α ))u -(α )α u ≤ (1 -(α ))|u | + w(α )|α ||u | ≤ a(α )|u | + w(α )|α ||u | (1.71)
From the uniform bound on E (u , α ), we get that the first term is bounded in L 2 (Ω), while (1.70) and the truncation Lemma 1 show that the second term is bounded in

L 1 (Ω) thus in L 2 (Ω). Finally, i. v is uniformly bounded in L ∞ (Ω)
ii. v is uniformly bounded in L 2 (Ω)

Chapter 1. Variational phase-field models of brittle fracture iii. J(v ) = ∅ invoking the Ambrosio's compactness theorem in SBV (in the Appendix A), we get that there exists v ∈ SBV (Ω) such that v → v strongly in L 2 (Ω). To conclude, since u = v (1-(α )) and α → 0 almost everywhere, we have,

u → u in L 2 (Ω)
Remark the proof above applies unchanged to the higher dimension case.

Gamma-convergence in 1d

The second part of the fundamental theorem of Γ-convergence requires that E Γ-converges to E. The definition of the Γ-convergence is in the Appendix A. The first condition means that E provides an asymptotic common lower bound for the E . The second condition means that this lower bound is optimal. The Γ-convergence is performed in 1d setting and is decomposed in two steps as follow: first prove the lower inequality, then construct the recovery sequence.

Lower semi-continuity inequality in 1d

We want to show that for any u ∈ SBV (Ω), and any (u , α ) such that u → u and α → 0 almost everywhere in Ω, we have,

lim inf →0 E (u , α , Ω) ≥ 1 2 Ω (u ) 2 dx + G c #(J(u)) (1.72)
Proof. Consider any interval I ⊂ Ω ⊂ R, such that,

lim inf →0 E (u , α , I) ≥ 1 2 I (u ) 2 dx if u ∈ W 1,2 (I) (1.73)
and,

lim inf →0 E (u , α , I) ≥ G c otherwise (1.74)
If lim inf →0 E (u , α , I) = ∞, both statements are trivial, so we can assume that there exist 0 ≤ C < ∞ such that,

lim inf →0 E (u , α , I) ≤ C (1.75)
We focus on (1.73) first, and assume that u ∈ W 1,2 (I). From (1.75) we deduce that w(α ) → 0 almost everywhere in I. Consequently, α → 0 almost everywhere in I. By Egoroff's theorem, for any > 0 there exists I ⊂ I such that |I | < and such that α → 0 uniformly on I \ I . For any δ > 0, thus we have,

1.3. Limit of the damage energy 1 -δ ≤ a(α ) on I \ I ,
for all and small enough, so that,

I\I (1 -δ) (u ) 2 dx ≤ I\I a(α ) (u ) 2 dx ≤ I a(α ) (u ) 2 dx (1.76)
Since u → u in W 1,2 (I) , and taking the lim inf on both sides, one gets,

(1 -δ) 2 I\I (u ) 2 dx ≤ lim inf →0 1 2 I a(α ) (u ) 2 dx (1.77)
we obtain the desired inequality (1.73) by letting → 0 and δ → 0.

To prove the second assertion (1.74), we first show that lim →0 sup x∈I α = 1, proceeding by contradiction. Suppose there exists δ > 0 such that α < 1δ on I. Then,

I a(1 -δ) (u ) 2 dx ≤ I a(α ) (u ) 2 dx
Taking the lim inf on both sides and using (1.75) , we get that,

lim inf →0 I (u ) 2 dx ≤ C a(1 -δ)
So u is uniformly bounded in W 1,2 (I), and therefore u ∈ W 1,2 (I), which contradicts our hypothesis. Reasoning as before, we have that α → 0 almost everywhere in I. Proceeding the same way on the interval (b , c ), one gets that,

lim inf →0 G c 4c w I w(α ) + (α ) 2 dx ≥ G c
which is (1.74). In order to obtain (1.72), we apply (1.74) on arbitrary small intervals centered around each points in the jump set of u and (1.73) on each remaining intervals in I.

Recovery sequence for the Γ-limit in 1d

The construction of the recovery sequence is more instructive. Given (u, α) we need to buid a sequence (u , α ) such that lim sup F (u , α ) ≤ F(u, α).

Proof. If F (u, α) = ∞, we can simply take u = u and α = α, so that we can safely assume that F(u, α) < ∞. As in the lower inequality, we consider the area near discontinuity points of u and away from them separately. Let (u, α) be given, consider an open interval I ⊂ R and a point x 0 ∈ J(u) ∩ I. Without loss of generality, we can assume that x 0 = 0 and I = (-δ, δ) for some δ > 0 . The construction of the recovery sequence is composed of two parts, first the recovery sequence for the damage, then one for the displacement.

The optimal damage profile obtained in the Lemma 2, directly gives,

lim sup →0 G c 4c w δ -δ w(α ) + (α ) 2 dx ≤ G c , (1.81)
this is the recovery sequence for the damage. Now, let's focus on the recovery sequence for the bulk term. We define b and

u (x) :=    x b u(x) if -b ≤ x ≤ b u(x) otherwise (1.82) Since a(α ) ≤ 1, we get that, -b -δ a(α ) (u ) 2 dx ≤ -b -δ (u ) 2 dx (1.
a(α ) (u ) 2 dx ≤ b -b (u ) 2 dx ≤ b -b u b + xu b 2 dx ≤ 2 b -b u b 2 dx + 2 b -b xu b 2 dx ≤ 2 b 2 b -b |u| 2 dx + 2 b -b (u ) 2 dx (1.85)
Since |u| ≤ M , the first term vanish when b → 0. Combining (1.83),(1.85) and (1.84).

Then, taking the lim sup on both sides and using I |u | 2 dx < ∞, we get that,

lim sup →0 1 2 δ -δ a(α ) (u ) 2 dx ≤ 1 2 δ -δ (u ) 2 dx (1.86)
Finally combining (1.81) and (1.86), one obtains

lim sup →0 δ -δ 1 2 a(α ) (u ) 2 + G c 4c w δ -δ w(α ) + (α ) 2 dx ≤ 1 2 δ -δ (u ) 2 dx + G c (1.87)
For the final construction of the recovery sequence, notice that we are free to assume that #(J(u)) is finite and chose δ ≤ inf{|x ix j |/2 s.t. x i , x j ∈ J(u), x i = x j }. For each x i ∈ J(u), we define I i = (x iδ, x i + δ) and use the construction above on each I i whereas on I \ I i we chose u = u and α linear and continuous at the end points of the I i . With this construction, is easy to see that α → 1 uniformly in I \ I i and that, lim sup

→0 I\ I i 1 2 a(α )(u ) 2 dx ≤ I (u ) 2 dx, (1.88)
and, lim sup

→0 I\ I i w(α ) + (α ) 2 dx = 0 (1.89)
Altogether, we obtain the upper estimate for the Γ-limit for pairs (u, 1) of finite energy, i.e.

lim sup →0 F (u , α ) ≤ F (u , 1) (1.90)

Extension to higher dimensions

To extend the Γ-limit to higher dimensions the lower inequality part is technical and is not developed here. But, the idea is to use Fubini's theorem, to build higher dimension by taking 1d slices of the domain, and use the lower continuity on each section see [START_REF] Ambrosio | Existence theory for a new class of variational problems[END_REF][START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]. The recovery sequence is more intuitive, a possible construction is to consider a smooth Γ ⊂ Ω and compute the distance to the crack J(u), such that,

d(x) = dist(x, J(u)) (1.91)
and let the volume of the region bounded by p-level set of d, such that,

s(y) = |{x ∈ R n ; d(x) ≤ y}| (1.92)
Figure 1.2: Iso distance to the crack J(u) for the level set b and δ

Following [START_REF] Evans | Measure theory and fine properties of functions[END_REF][START_REF] Evans | On the partial regularity of energy-minimizing, areapreserving maps[END_REF], the co-area formula from Federer [START_REF] Federer | Geometric measure theory[END_REF] is,

Ω f (x) ∇g(x) dx = +∞ -∞ g -1 (y) f (x)dH n-1 (x) dy (1.93)
In particular, taking g(x) = d(x) which is 1-Lipschitz, i.e. ∇d(x) = 1 almost everywhere. We get surface s(y),

s(y) = s(y) ∇d(x) dx = y 0 H n-1 ({x; d(x) = t})dt (1.94)
and

s (y) = H n-1 ({x; d(x) = y}) (1.95)
In particular,

s (0) = lim y→0 s(y) y = 2H n-1 (J(u)) (1.96)

Limit of the damage energy

Consider the damage,

α (d(x)) :=      1 if d(x) ≤ b γ (d(x)) if b ≤ d(x) ≤ δ 0 otherwise (1.97)
The surface energy term is,

Ω w(α ) + |∇α | 2 dx = 1 d(x)≤b dx + b ≤d(x)≤δ w(α (d(x))) + |∇α (d(x))| 2 dx (1.98)
The first integral term, is the surface bounded by the iso-contour distant b from the crack, i.e

s(b ) = d(x)≤b dx = b 0 H n-1 ({x; d(x) = y}) dy (1.
≤ δ/ 0 w(α (x )) + α (x ) 2 s (x ) dx (1.102)
Passing the limit → 0 and using the Remark 1 on the optimal profile invariance, we get,

lim sup →0 G c 4c w Ω w(α (x)) + |∇α (x)| 2 dx ≤ G c H n-1 (J(u)) (1.103)
For the bulk term, consider the displacement,

u (x) :=    d(x) b u(x) if d(x) ≤ b u(x) otherwise (1.104)
Similarly to the 1d, one gets,

lim sup →0 Ω 1 2 a(α )(∇u ) 2 dx ≤ Ω 1 2 (∇u ) 2 dx (1.105)
Therefore,

lim sup →0 Ω 1 2 a(α )(∇u ) 2 dx + G c 4c w Ω w(α + |∇α | 2 dx ≤ Ω 1 2 (∇u ) 2 dx + G c H n-1 (J(u))
(1.106)

Numerical implementation

In a view to numerically implement gradient damage models, it is common to consider time and space discretization. Let's first focus on the time-discrete evolution, by considering a time interval [0, T ] subdivided into (N + 1) steps, such that,

0 = t 0 < t 1 < • • • < t i-1 < t i < • • • < t N = T .
At any step i, the sets of admissible displacement and damage fields C i and D i are, For any i find (u i , α i ) ∈ (C i , D i ) that satisfies the discrete evolution by local minimizer if the following hold: i. Initial condition α t 0 = α 0 and u t 0 = u 0 ii. For some

C i := u ∈ H 1 (Ω) : u = ūi on ∂ D Ω D i := β ∈ H 1 (Ω) : α i-1 (x) ≤ β ≤ 1, ∀x ∈ Ω , (1.107
h i > 0, find (u i , α i ) ∈ (C i , D i ), such that, (v, β) -(u i , α i ) ≤ h i , E (u i , α i ) ≤ E (v, β) (1.108)
where,

E (u, α) = Ω 1 2 a(α)Ae(u) : e(u) dx - ∂ N Ω g(t) • u dH n-1 + G c 4c w Ω w(α) + |∇α| 2 dx (1.109)
One observes that our time-discretization evolution do not enforce energy balance. Since a(α) and w(α) are convex, the total energy E (u, α) is separately convex with respect to u and α, but that is not convex. Hence, a proposed alternate minimization algorithm guarantees to converge to a critical point of the energy satisfying the irreversibility condition [START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF][START_REF] Burke | An adaptive finite element approximation of a variational model of brittle fracture[END_REF]. The idea is for each time-step t i , we minimize the problem with respect to any kinematic admissible u for a given α, then, fixed u and minimize E (u, α) with respect to α subject to the irreversibility α i ≥ α i-1 , repeat the procedure until the variation of the damage is small. This gives the following algorithm see Algorithm 1, where δ α is a fixed tolerance parameter.

For the space discretization of E (u, α), we use the finite element methods considering linear Lagrange elements for u and α. To solve the elastic problem preconditioned conjugate gradient solvers is employed, and the constraint minimization with respect to the damage is implemented using the variational inequality solvers provided by PETSc [START_REF] Balay | PETSc Web page[END_REF][START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF]. All computations were performed using the open source mef902 .

Due to the non-convexity of E , solution satisfying irreversibility and stationarity might not be unique.

For remainder solutions, a study selection can be performed. For instance looking at solutions which satisfy the energy balance, or selecting displacement and damage fields which are continuous in time. Another way is to compare results with all previous one in order to avoid local minimizers solution (see [START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Bourdin | The variational formulation of brittle fracture: numerical implementation and extensions[END_REF] for more details on the backtracking idea). This method will select global minimizers from the set of solutions.

1: Let j = 0 and α 0 := α i-1 2: repeat 3:

Compute the equilibrium,

u j+1 := argmin u∈C i E (u, α j ) 4:
Compute the damage,

α j+1 := argmin α∈D i α≥α i-1 E (u j+1 , α) 5: j := j + 1 6: until α j -α j-1 L ∞ ≤ δ α 7: Set, u i := u j and α i := α j

Conclusion

The strength of the phase fields models to brittle fracture is the variational structure of the model conceived as an approximation of Griffith and its evolution based on three principles: irreversibility of the damage, stability and energy balance of the total energy. A fundamental property of the model is the maximum admissible stress illustrated in the one dimensional example. This also constrained the damage thickness size, since it governs . Numerically the fracture path is obtained by alternate searching of the damage trajectory which decreases the total energy and the elastic solution of the problem.

Appendix A

Theorem 2 (Ambrosio's compactness and lower semicontinuity on SBV)

Let (f n ) n be a sequence of functions in SBV (Ω) such that there exists non-negative constants

C 1 , C 2 and C 3 such that, i. f n is uniformly bounded in L ∞ (Ω) ii. ∇f n is uniformly bounded in L q (Ω, R n ) with q > 1 iii. H n-1 (J(f n )) is uniformly bounded Then, there exists f ∈ SBV (Ω) and a subsequence f k(n) such that, i. f k(n) → f strongly in L p (Ω), for all p < ∞ ii. ∇f k(n) → ∇f weakly in L q (Ω; R n ) iii. H n-1 (J(f )) ≤ lim inf n H n-1 (J(f n )) Theorem 3 (Fundamental theorem of Γ-convergence)
If E Γ -converges to E, u is a minimizer of E , and (u ) is compact in X, then there exists u ∈ X such that u is a minimizer of E, u → u, and E (u ) → E(u).

Definition 6 (Γ-convergence)

Let E : X → R and E : X → R, where X is a topological space. Then E Γ converges to E if the following two conditions hold for any u ∈ X i. Lower semi continuity inequality: for every equence (u ) ∈ Xsuch that u → u

E(u) ≤ lim inf →0 E (u ),
ii. Existence of a recovery sequence: there exists a sequence (u )

∈ X with u → u such that lim sup →0 E (u ) ≤ E(u).
Chapter 2

Crack nucleation in variational phase-field models of brittle fracture

Despite its many successes, Griffith's theory of brittle fracture [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] and its heir, Linear Elastic Fracture Mechanics (LEFM), still faces many challenges. In order to identify a crack path, additional branching criteria whose choice is still unsettled have to be considered. Accounting for scale effects in LEFM is also challenging, as illustrated by the following example: consider a reference structure of unit size rescaled by a factor L.

The critical loading at the onset of fracture scales then as 1/ √ L, leading to a infinite nucleation load as the structure size approaches 0, which is inconsistent with experimental observation for small structures [START_REF] Bažant | Scaling of quasibrittle fracture: asymptotic analysis[END_REF][START_REF] Issa | Size effects in concrete fracture: Part I, experimental setup and observations[END_REF][START_REF] Chudnovsky | Slow crack growth, its modeling and crack-layer approach: A review[END_REF].

It is well accepted that this discrepancy is due to the lack of a critical stress (or a critical lengthscale) in Griffith's theory. Yet, augmenting LEFM to account for a critical stress is very challenging. In essence, the idea of material strength is incompatible with the concept of elastic energy release rate near stress singularity, the pillar of Griffith-like theories, as it would imply crack nucleation under an infinitesimal loading. Furthermore, a nucleation criterion based solely on pointwise maximum stress will be unable to handle crack formation in a body subject to a uniform stress distribution.

Many approaches have been proposed to provide models capable of addressing the aforementioned issues. Some propose to stray from Griffith fundamental hypotheses by incorporating cohesive fracture energies [START_REF] Ortiz | Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis[END_REF][START_REF] Del Piero | A diffuse cohesive energy approach to fracture and plasticity: the one-dimensional case[END_REF][START_REF] De Borst | Cohesive-zone models, higher-order continuum theories and reliability methods for computational failure analysis[END_REF][START_REF] Charlotte | Initiation of cracks with cohesive force models: a variational approach[END_REF] or material non-linearities [START_REF] Gou | Modeling fracture in the context of a strain-limiting theory of elasticity: A single plane-strain crack[END_REF]. Others have proposed dual-criteria involving both elastic energy release rate and material strength such as [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF], for instance. Models based on the peridynamics theory [START_REF] Silling | Reformulation of elasticity theory for discontinuities and long-range forces[END_REF] may present an alternative way to handle these issues, but to our knowledge, they are still falling short of providing robust quantitative predictions at the structural scale.

Francfort and Marigo [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] set to devise a formulation of brittle fracture based solely on Griffith's idea of competition between elastic and fracture energy, yet capable of handling the issues of crack path and crack nucleation. However, as already pointed-out in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], their model inherits a fundamental limitation of the Griffith theory and LEFM: the lack of an internal length scale and of maximum allowable stresses.

Amongst many numerical methods originally devised for the numerical implemen-tation of the Francfort-Marigo model [START_REF] Bourdin | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Negri | Numerical minimization of the Mumford-Shah functional[END_REF][START_REF] Fraternali | Free discontinuity finite element models in two-dimensions for inplane crack problems[END_REF][START_REF] Schmidt | Eigenfracture: An eigendeformation approach to variational fracture[END_REF], Ambrosio-Tortorelli regularizations [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], originally introduced in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], have become ubiquitous. They are known nowadays as phase-field models of fracture, and share several common points with the approaches coming from Ginzburg-Landau models for phase-transition [START_REF] Karma | Phase-field model of mode III dynamic fracture[END_REF]. They have been applied to a wide variety of fracture problems including fracture of ferro-magnetic and piezo-electric materials [START_REF] Abdollahi | Phase-field modeling of crack propagation in piezoelectric and ferroelectric materials with different electromechanical crack conditions[END_REF][START_REF] Wilson | A phase-field model for fracture in piezoelectric ceramics[END_REF], thermal and drying cracks [START_REF] Maurini | Crack patterns obtained by unidirectional drying of a colloidal suspension in a capillary tube: experiments and numerical simulations using a two-dimensional variational approach[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF], or hydraulic fracturing [START_REF] Bourdin | A variational approach to the numerical simulation of hydraulic fracturing[END_REF][START_REF] Wheeler | An augmented-lagrangian method for the phase-field approach for pressurized fractures[END_REF][START_REF] Chukwudozie | Application of the Variational Fracture Model to Hydraulic Fracturing in Poroelastic Media[END_REF][START_REF] Wilson | Phase-field modeling of hydraulic fracture[END_REF] to name a few. They have been expended to account for dynamic effects [START_REF] Larsen | Existence of solutions to a regularized model of dynamic fracture[END_REF][START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF][START_REF] Hofacker | A phase field model of dynamic fracture: Robust field updates for the analysis of complex crack patterns[END_REF], ductile behavior [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. coupled brittle-to-ductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF], cohesive effects [START_REF] Crismale | Viscous approximation of quasistatic evolutions for a coupled elastoplastic-damage model[END_REF][START_REF] Conti | Phase field approximation of cohesive fracture models[END_REF][START_REF] Freddi | Numerical insight of a variational smeared approach to cohesive fracture[END_REF], large deformations [START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Miehe | Phase field modeling of ductile fracture at finite strains. a variational gradient-extended plasticity-damage theory[END_REF][START_REF] Borden | A phasefield formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF], or anisotropy [START_REF] Li | Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy[END_REF], for instance.

Although phase-field models were originally conceived as approximations of Francfort and Marigo's variational approach to fracture in the vanishing limit of their regularization parameter, a growing body of literature is concerned with their links with gradient damage models [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF]. In this setting, the regularization parameter is kept fixed and interpreted as a material's internal length [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF][START_REF] Freddi | Regularized variational theories of fracture: A unified approach[END_REF][START_REF] Del | A variational approach to fracture and other inelastic phenomena[END_REF]. In particular, [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF] proposed an evolution principle for an Ambrosio-Tortorelli like energy based on irreversibility, stability and energy balance, where the regularization parameter is kept fixed and interpreted as a material's internal length. This approach, which we refer to as variational phase-field models, introduces a critical stress proportional to 1/ √ . As observed in [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF], it can potentially reconcile stress and toughness criteria for crack nucleation, recover pertinent size effect at small and large length-scales, and provide a robust and relatively simple approach to model crack propagation in complex two-and three-dimensional settings. However, the few studies providing experimental verifications [START_REF] Pham | Experimental validation of a phase-field model for fracture[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF] are still insufficient to fully support this conjecture.

The goal of this chapter is precisely to provide such evidences, focusing on nucleation and size-effects for mode-I cracks. We provide quantitative comparisons of nucleation loads near stress concentrations and singularities with published experimental results for a range of materials. We show that variational phase-field models can reconcile strength and toughness thresholds and account for scale effect at the structural and the material length-scale. In passing, we leverage the predictive power of our approach to propose a new way to measure a material's tensile strength from the nucleation load of a crack near a stress concentration or a weak singularity. In this study, we focus solely on the identification of the critical stress at the first crack nucleation event and are not concerned by the post-critical fracture behavior.

The chapter is organized as follows: in Section 2.1, we introduce variational phasefield models and recall some of their properties. Section 2.2 focuses on the links between stress singularities or concentrations and crack nucleation in these models. We provide validation and verification results for nucleation induced by stress singularities using Vshaped notches, and concentrations using U-notches. Section 2.3 is concerned with shape and size effects. We investigate the role of the internal length on nucleation near a defect, focusing on an elliptical cavity and a mode-I crack, and discussing scale effects at the material and structural length scales.

Chapter 2. Crack nucleation in variational phase-field models of brittle fracture

Variational phase-field models

We start by recalling some important properties of variational phase-field models, focussing first on their construction as approximation method of Francfort and Marigo's variational approach to fracture, then on their alternative formulation and interpretation as gradient-damage models.

Regularization of the Francfort-Marigo fracture energy

Consider a perfectly brittle material with Hooke's law A and critical elastic energy release rate G c occupying a region Ω ⊂ R n , subject to a time dependent boundary displacement ū(t) on a part ∂ D Ω of its boundary and stress-free on the remainder ∂ N Ω. In the variational approach to fracture, the quasi-static equilibrium displacement u i and crack set Γ i at a given discrete time step t i are given by the minimization problem (see also [START_REF] Bourdin | The variational approach to fracture[END_REF])

(u i , Γ i ) = argmin u=ū i on ∂ D Ω Γ⊃Γ i-1 E(u, Γ) := Ω\Γ 1 2 Ae(u) • e(u) dx + G c H n-1 (Γ ∩ Ω \ ∂ N Ω), (2.1)
where H n-1 (Γ) denotes the Hausdorff n -1-dimensional measure of the unknown crack Γ, i.e. its aggregate length in two dimensions or surface area in three dimensions, and e(u) := 1 2 (∇u + ∇ T u) denotes the symmetrized gradient of u. Because in (2.1) the crack geometry Γ is unknown, special numerical methods had to be crafted. Various approaches based for instance on adaptive or discontinuous finite elements were introduced [START_REF] Bourdin | Implementation of an adaptive finite-element approximation of the Mumford-Shah functional[END_REF][START_REF] Giacomini | A discontinuous finite element approximation of quasi-static growth of brittle fractures[END_REF][START_REF] Fraternali | Free discontinuity finite element models in two-dimensions for inplane crack problems[END_REF]. Variational phase-field methods, take their roots in Ambrosio and Tortorelli's regularization of the Mumford-Shah problem in image processing [START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF][START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF], adapted to brittle fracture in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. In this framework, a regularized energy E depending on a regularization length > 0 and a "phase-field" variable α taking its values in [0, 1] is introduced. A broad class of such functionals was introduced in [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF]. They are

E (u, α) = Ω a(α) + η 2 Ae(u) • e(u) dx + G c 4c w Ω w(α) + |∇α| 2 dx, (2.2)
where a and w are continuous monotonic functions such that a(0) = 1, a(1) = 0, w(0) = 0, and w(1) = 1, η = o( ), and c w := 1 0 w(s) ds is a normalization parameter. The approximation of E by E takes place with the framework of Γ-convergence (see [START_REF] Maso | An introduction to Γ-convergence[END_REF][START_REF] Braides | Gamma-convergence for Beginners[END_REF] for instance). More precisely, if E Γ-converges to E, then the global minimizers of E converge to that of E. The Γ-convergence of a broad class of energies, including the ones above was achieved with various degrees of refinement going from static scalar elasticity to time discrete and time continuous quasi-static evolution linearized elasticity, and their finite element discretization [START_REF] Bellettini | Discrete approximation of a free discontinuity problem[END_REF][START_REF] Bourdin | Image segmentation with a finite element method[END_REF][START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF][START_REF] Giacomini | A discontinuous finite element approximation of quasi-static growth of brittle fractures[END_REF][START_REF] Chambolle | An approximation result for special functions with bounded variations[END_REF][START_REF] Chambolle | Addendum to "An Approximation Result for Special Functions with Bounded Deformation[END_REF][START_REF] Giacomini | Ambrosio-Tortorelli approximation of quasi-static evolution of brittle fractures[END_REF][START_REF] Burke | An adaptive finite element approximation of a variational model of brittle fracture[END_REF][START_REF] Burke | An adaptive finite element approximation of a generalized Ambrosio-Tortorelli functional[END_REF][START_REF] Iurlano | A density result for gsbd and its application to the approximation of brittle fracture energies[END_REF].

Throughout this chapter, we focus on two specific models:

E (u, α) = Ω (1 -α) 2 + η 2 Ae(u) • e(u) dx + G c 2 Ω α 2 + |∇α| 2 dx, (AT 2 )
2.1. Variational phase-field models introduced in [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF] for the Mumford-Shah problem and in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] for brittle fracture, and

E (u, α) = Ω (1 -α) 2 + η 2 Ae(u) • e(u) dx + 3G c 8 Ω α + |∇α| 2 dx (AT 1 )
used in [START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF].

The "surfing" problem introduced in [START_REF] Hossain | Effective toughness of heterogeneous media[END_REF] consists in applying a translating boundary displacement on ∂Ω given by ū(x, y) = ūI (x-V t, y), where ūI denotes the asymptotic farfield displacement field associated with a mode-I crack along the x-axis with tip at (0, 0), V is a prescribed loading "velocity", and t a loading parameter ("time"). ] with an initial crack Γ 0 = [0, l 0 ] × {0} for several values of . The AT 1 model is used, assuming plane stress conditions, and the mesh size h is adjusted so that /h = 5, keeping the "effective" numerical toughness G eff := G c 1 + h 4cw fixed (see [START_REF] Bourdin | The variational approach to fracture[END_REF]). The Poisson ratio is ν = 0.3, the Young's modulus is E = 1, the fracture toughness is G c = 1.5, and the loading rate V = 4. As expected, after a transition stage, the crack length depends linearly on the loading parameter with slope 3.99, 4.00 and 4.01 for =0.1, 0.05 and 0.025 respectively. The elastic energy release rate G, computed using the G θ method [START_REF] Destuynder | Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile[END_REF][START_REF] Sicsic | From gradient damage laws to Griffith's theory of crack propagation[END_REF][START_REF] Li | Gradient damage modeling of brittle fracture in an explicit dynamics context[END_REF] is very close to G eff . Even though Γ-convergence only mandates that the elastic energy release rate in the regularized energy converges to that of Griffith as → 0, we observe that as long as is "compatible" with the discretization size and domain geometry, its influence on crack propagation is insignificant. Similar observations were reported in [START_REF] Klinsmann | An assessment of the phase field formulation for crack growth[END_REF][START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF][START_REF] Pham | Experimental validation of a phase-field model for fracture[END_REF]. Figure 2.1(right) repeats the same experiment for a curve propagating along a circular path. Here, the boundary displacement is given by Muskhelishvili's exact solution for a crack propagating in mode-I along a circular path [START_REF] Muskhelishvili | Some Basic Problems of the Mathematical Theory of Elasticity: Fundamental Equations, Plane Theory of Elasticity, Torsion, and Bending (translated from Russian)[END_REF]. The Young's modulus, fracture toughness, and loading rate are set to 1. Again, we see that even for a fixed regularization length, the crack obeys Griffith's criterion.

Chapter 2. Crack nucleation in variational phase-field models of brittle fracture When crack nucleation is involved, the picture is considerably different. Consider a one-dimensional domain of length L, fixed at one end and submitted to an applied displacement ū = e L at the other end. For the lack of an elastic singularity, LEFM is incapable of predicting crack nucleation here, and predicts a structure capable of supporting arbitrarily large loads without failing. A quick calculation shows that the global minimizer of (2.1) corresponds to an uncracked elastic solution if e < e c := 2Gc EL , while at e = e c , a single crack nucleates at an arbitrary location (see [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]). The failure stress is σ c = 2G c E/L, which is consistent with the scaling law σ c = O 1/ √ L mentioned in the introduction. The uncracked configuration is always a stable local minimizer of (2.1), so that if local minimization of (2.1) is considered, nucleation never takes place. Just as before, one can argue that for the lack of a critical stress, an evolution governed by the generalized Griffith energy (2.1) does not properly account for nucleation and scaling laws.

When performing global minimization of (2.2) using the backtracking algorithm of [START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF] for instance, a single crack nucleates at an -dependent load. As predicted by the Γ-convergence of E to E, the critical stress at nucleation converges to 2G c E/L as → 0. Local minimization of (2.2) using the alternate minimizations algorithm of [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], or presumably any gradient-based monotonically decreasing scheme, leads to the nucleation of a single crack at a critical load e c , associated with a critical stress σ c = O G c E/ , as described in [START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF] for example. In the limit of vanishing , local and global minimization of (2.2) inherit therefore the weaknesses of Griffith-like theories when dealing with scaling properties and crack nucleation.

Variational phase-field models as gradient damage models

More recent works have sought to leverage the link between σ c and . Ambrosio-Tortorelli functionals are then seen as the free energy of a gradient damage model [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Lorentz | Analysis of non-local models through energetic formulations[END_REF][START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF] where α plays the role of a scalar damage field. In [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], a thorough investigation of a one-dimensional tension problem led to interpreting as a material's internal or characteristic length linked to a material's tensile strength. An overview of this latter approach, which is the one adopted in the rest of this work, is given below.

In all that follows, we focus on a time-discrete evolution but refer the reader to [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF] for a time-continuous formulation which can be justified within the framework of generalized standard materials [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF] and rate-independent processes [START_REF] Mielke | Evolution of rate-independent systems[END_REF]. At any time step i > 1, the sets of admissible displacement and damage fields C i and D i , equipped with their natural H 1 norm, are

C i = u ∈ H 1 (Ω) : u = ūi on ∂ D Ω , D i = β ∈ H 1 (Ω) : α i-1 (x) ≤ β(x) ≤ 1, ∀x ∈ Ω ,
where the constraint α i-1 (x) ≤ β(x) ≤ 1 in the definition of D i mandates that the damage be an increasing function of time, accounting for the irreversible nature of the 2.1. Variational phase-field models damage process. The damage and displacement fields (u i , α i ) are then local minimizers of the energy E , i.e. there exists h i > 0 such that

∀(v, β) ∈ C i × D i such that (v, β) -(u i , α i ) ≤ h i , E (u i , α i ) ≤ E (v, β), (2.3)
where • denotes the natural H 1 norm of C i × D i . We briefly summarize the solution of the uniaxial tension of a homogeneous bar [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], referring the reader to the recent review [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF] for further details: As one increases the applied strain, the damage field remains 0 and the stress field constant until it reaches the elastic limit

σ e = G c E c w w (0) 2s (0) . (2.4)
where E is the Young modulus of the undamaged material, and s(α) = 1/a(α). If the applied displacement is increased further, the damage field increases but remains spatially constant. Stress hardening is observed until peak stress σ c , followed by stress softening. A stability analysis shows that for long enough domains (i.e. when L ), the homogeneous solution is never stable in the stress softening phase, and that a snapback to a fully localized solution such that max x∈(0,L) α(x) = 1 is observed. The profile of the localized solution and the width D of the localization can be derived explicitly from the functions a and w. With the choice of normalization of (2.2), the surface energy associated to the fully localized solution is exactly G c and its elastic energy is 0, so that the overall response of the bar is that of a brittle material with toughness G c and strength σ c .

Knowing the material's toughness G c and the Young's modulus E, one can then adjust in such a way that the peak stress σ c matches the nominal material's strength.

Let us denote by

ch = G c E σ 2 c = K 2 Ic σ 2 c (2.5)
the classical material's characteristic length (see [START_REF] Rice | The mechanics of earthquake rupture[END_REF][START_REF] Falk | A critical evaluation of cohesive zone models of dynamic fracture[END_REF], for instance), where E = E in three dimensions and in plane stress, or E = E/(1ν 2 ) in plane strain, and K Ic = √ G c E is the mode-I critical stress intensity factor. The identification above gives

1 := 3 8 ch ; 2 := 27 256 ch , (2.6)
for the AT 1 and AT 2 models, respectively. Table 2.1 summarizes the specific properties of the AT 1 and AT 2 models. The AT 1 model has some key conceptual and practical advantages over the AT 2 model used in previous works, which were leveraged in [START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF] for instance: It has a non-zero elastic limit, preventing diffuse damage at small loading. The length localization band D is finite so that equivalence with Griffith energy is obtained even for a finite value of , and not only in the limit of → 0, as predicted by Γ-convergence [START_REF] Sicsic | From gradient damage laws to Griffith's theory of crack propagation[END_REF]. By remaining quadratic in the α and u variables, its numerical implementation using alternate minimizations originally introduced in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] is very efficient.

Chapter 2. Crack nucleation in variational phase-field models of brittle fracture

Model w(α) a(α) c w σ e σ c D ch AT 1 α (1 -α) 2 2 3 3GcE 8 3GcE 8 4 8 3 AT 2 α 2 (1 -α) 2 1 2 0 3 16 3GcE ∞ 256 27
Table 2.1: Properties of the gradient damage models considered in this work: the elastic limit σ e , the material strength σ c , the width of the damage band D, and the conventional material length ch defined in (2.5). We use the classical convention E = E in three dimension and in plane stress, and

E = E/(1 -ν 2 ) in plane strain.
In all the numerical simulations presented below, the energy (2.2) is discretized using linear Lagrange finite elements, and minimization performed by alternating minimization with respect to u and α. Minimization with respect to u is a simple linear problem solved using preconditioned gradient conjugated while constrained minimization with respect to α is reformulated as a variational inequality and implemented using the variational inequality solvers provided by PETSc [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | PETSc Web page[END_REF]. All computations were performed using the open source implementations mef901 and gradient-damage2 .

Effect of stress concentrations

The discussion above suggests that variational phase-field models, as presented in Section 2.1.2, can account for strength and toughness criteria simultaneously, on an idealized geometry. We propose to investigate this claim further by focusing on more general geometries, a V-shaped notch to illustrate nucleation near stress singularities and a Ushaped notch for stress concentrations. There is a wealth of experimental literature on crack initiation in such geometries using three-point bending (TPB), four-point bending (FPB), single or double edge notch tension (SENT and DENT) allowing us to provide qualitative validation and verification simulations of the critical load at nucleation.

Initiation near a weak stress singularity: the V-notch

Consider a V-shaped notch in a linear elastic isotropic homogeneous material. Let (r, θ) be the polar coordinate system emanating from the notch tip with θ = 0 corresponding to the notch symmetry axis, shown on Figure 2.2(left). Assuming that the notch lips Γ + ∪ Γ -are stress-free, the mode-I component of the singular part of the stress field in 2.2. Effect of stress concentrations plane strain is given in [START_REF] Leguillon | Computation of Singular Solutions in Elliptic Problems and Elasticity[END_REF]:

σ θθ = kr λ-1 F (θ), σ rr = kr λ-1 F (θ) + (λ + 1)F (θ) λ(λ + 1) , σ rθ = -kr λ-1 F (θ) (λ + 1) , (2.7)
where

F (θ) = (2π) λ-1 cos((1 + λ)θ) -f (λ, ω) cos((1 -λ)θ) 1 -f (λ, ω) , (2.8)
and

f (λ, ω) = (1 + λ) sin((1 + λ)(π -ω)) (1 -λ) sin((1 -λ)(π -ω)) , (2.9)
and the exponent of the singularity λ ∈ [1/2, 1], see (2.11)

Note that this definition differs from the one often encountered in the literature by a factor (2π) λ-1 , so that when ω = 0 (i.e. when the notch degenerates into a crack), k corresponds to the mode-I stress intensity factor whereas when ω = π/2, k is the tangential stress, and that the physical dimension of [k] ≡ N/m -λ -1 is not a constant but depends on the singularity power λ.

If ω < π/2 (i.e. ω > π/2), the stress field is singular at the notch tip so that a nucleation criterion based on maximum pointwise stress will predict crack nucleation for any arbitrary small loading. Yet, as long as ω > 0 (ω < π), the exponent of the singularity is sub-critical in the sense of Griffith, so that LEFM forbids crack nucleation, regardless of the magnitude of the loading.

ūr = r λ E (1 -ν 2 )F (θ) + (λ + 1)[1 -νλ -ν 2 (λ + 1)]F (θ) λ 2 (λ + 1) ūθ = r λ E (1 -ν 2 )F (θ) + [2(1 + ν)λ 2 + (λ + 1)(1 -νλ -ν 2 (λ + 1)]F (θ) λ 2 (1 -λ 2 ) .
(2.12)

In the mode-I Pac-Man test, we apply a boundary displacement on the outer edge of the domain ∂ D Ω of the form tū on both components of u, t being a monotonically increasing loading parameter. We performed series of numerical simulations varying the notch angle ω and regularization parameter for the AT 1 and AT 2 models. Up to a rescaling and without loss of generality, it is always possible to assume that E = 1 and G c = 1. The Poisson ratio was set to ν = 0.3. We either prescribed the value of the damage field on Γ + ∪ Γ -to 1 (we refer this to as "damaged notch conditions") or let it free ("undamaged notch conditions"). The mesh size was kept at a fixed ratio of the internal length h = /5.

For "small" enough loadings, we observe an elastic or nearly elastic phase during which the damage field remains 0 or near 0 away from an area of radius o( ) near the notch tip. Then, for some loading t = t c , we observed the initiation of a "large" add-crack associated with a sudden jump of the elastic and surface energy. Figure 2.3 shows a typical mesh, the damage field immediately before and after nucleation of a macroscopic crack and the energetic signature of the nucleation event.

Figure 2.4 shows that up to the critical loading, the generalized stress intensity factor can be accurately recovered by averaging σ θθ (r, 0)/(2π r) λ-1 along the symmetry axis of the domain, provided that the region r ≤ 2 be excluded.

Figure 2.5(left) shows the influence of the internal length on the critical generalized stress intensity factor for a sharp notch (ω = 0.18°) for the AT 1 and AT 2 models, using damaged and undamaged notch boundary conditions on the damage field. In this case, with the normalization (2.11), the generalized stress intensity factor coincides with the standard mode-I stress intensity factor K Ic . As suggested by the surfing experiment in

t := k AT c K Ic = √ G c E .
As reported previously in [START_REF] Klinsmann | An assessment of the phase field formulation for crack growth[END_REF] for instance, undamaged notch conditions lead to overestimating the critical load. We speculate that this is because with undamaged notch condition, the energy barrier associated with bifurcation from an undamaged (or partially damaged) state to a fully localized state needs to be overcome. As expected, this energy barrier is larger for the AT 1 model than for the AT 2 model for which large damaged areas ahead of the notch tip are observed.

For flat notches (2ω = 179.64°) as shown in Figure 2.5(right), the generalized stress intensity factor k takes the dimension of a stress, and crack nucleation is observed when k c reaches the -dependent value σ c given in Table 2.1, i.e. when σ θθ | θ=0 = σ c , as in the uniaxial tension problem. In this case the type of damage boundary condition on the notch seems to have little influence. For intermediate values of ω, we observe in Figure 2.6 that the critical generalized stress intensity factor varies smoothly and monotonically between its extreme values and remains very close to K Ic for opening angles as high as 30°, which justifies the common numerical practice of replacing initial cracks with slightly open sharp notches and damaged notch boundary conditions. See Table 2.3 for numerical data.

k c /(K Ic ) ef f ω = 0.18 • AT 1 -U AT 1 -D AT 2 -U AT 2 -D 10 -1 10 0 ch 1 2 4 6 8 k c 1 -0.5 ω = 89, 82 • AT 1 -U AT 1 -D AT 2 -U AT 2 -D

Validation

For intermediate values 0 < 2ω < π, we focus on validation against experiments from the literature based on measurements of the generalized stress intensity factor at a V-shaped notch. Data from single edge notch tension (SENT) test of soft annealed tool steel, (AISI O1 at -50 • C) [START_REF] Strandberg | Fracture at V-notches with contained plasticity[END_REF], four point bending (FPB) experiments of Divinycell® H80, H100, H130, and H200 PVC foams) [START_REF] Grenestedt | On cracks emanating from wedges in expanded PVC foam[END_REF], and double edge notch tension (DENT) experiments of poly methyl methacrylate (PMMA) and Duraluminium [START_REF] Seweryn | Brittle fracture criterion for structures with sharp notches[END_REF], were compiled in [START_REF] Gómez | A fracture criterion for sharp V-notched samples[END_REF]. We performed a series of numerical simulations of Pac-Man tests using the material properties reported in [START_REF] Gómez | A fracture criterion for sharp V-notched samples[END_REF] and listed in Table 2.2. In all cases, the internal length was computed using (2.6). Plexiglass DENT (Seweryn) numerical simulations with experimental values reported in the literature for V-notch with varying aperture. The definition (2.11) for k is used. For the AT 1 model, we observe a good agreement for the entire range of notch openings, as long as damaged notch conditions are used for small notch angles and undamaged notch conditions for large notch angles. For the AT 2 model, the same is true, but the agreement is not as good for large notch angles, due to the presence of large areas of distributed damage prior to crack nucleation.

Effect of stress concentrations

Material E ν K Ic σ c source [MPa] [MPa √ m] [MPa] Al 2 O 3 -
k c [MPa.m 1-λ ] Steel SENT (Strandberg) AT 1 -U AT 1 -D AT 2 -U AT 2 -D 0
AT 1 -U AT 1 -D AT 2 -U AT 2 -D
The numerical values of the critical generalized stress intensity factors for the AT 1 models and the experiments from the literature are included in Tables 2.4, 2.5, 2.6, and 2.7 using the convention of (2.11) for k. As suggested by Figure 2.5 and reported in the literature see [START_REF] Klinsmann | An assessment of the phase field formulation for crack growth[END_REF], nucleation is best captured if damaged notch boundary conditions are used for sharp notches and undamaged notch conditions for flat ones.

These examples strongly suggest that variational phase-field models of fracture are capable of predicting mode-I nucleation in stress and toughness dominated situations, as seen above, but also in the intermediate cases. Conceptually, toughness and strength (or equivalently internal length) could be measured by matching generalized stress intensity factors in experiments and simulations. When doing so, however, extreme care has to be exerted in order to ensure that the structural geometry has no impact on the measured generalized stress. Similar experiments were performed in [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF][START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] for three and four point bending experiments on PMMA and Aluminum oxide-Zirconia ceramics samples. While the authors kept the notch angle fixed, they performed three and four point bending experiments or varied the relative depth of the notch as a fraction of the sample height (see Figure 2.9). Figure 2.9: Schematic of the geometry and loading in the four point bending experiments of [START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] (left) and three point bending experiments of [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] (right). The geometry of the three point bending experiment of [START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] is identical to that of their four point bending, up to the location of the loading devices.

Figure 2.10 compares numerical values of the generalized stress intensity factor using the AT 1 model with experimental measurements, and the actual numerical values are included in Table 2.8 and 2.9.

For the Aluminum oxide-Zirconia ceramic, we observe that the absolute error between measurement and numerical prediction is typically well within the standard deviation of the experimental measurement. As expected, damaged notch boundary conditions lead Chapter 2. Crack nucleation in variational phase-field models of brittle fracture to better approximation of k c for small angles, and undamaged notches are better for larger values of ω.

k c [MPa.m 1-λ ] Al 2 O 3 -7%ZrO 2 FPB (Yoshibash) Al 2 O 3 -7%ZrO 2 TPB (Yoshibash) AT 1 -U AT 1 -D 20
For the three point bending experiments in PMMA of [START_REF] Dunn | Fracture initiation at sharp notches: correlation using critical stress intensities[END_REF] later reported in [START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF], the experimental results suggest that the relative depth a/h of the notch has a significant impact on k c . We therefore performed full-domain numerical simulation using the geometry and loading from the literature, and compared the critical force upon which a crack nucleates in experiments and simulations. All computations were performed using the AT 1 model in plane strain with undamaged notch boundary conditions. Figure 2.11 compares the experimental and simulated value of the critical load at failure, listed in Table 2.10 and 2.11.

These simulations show that a robust quantitative prediction of the failure load in geometries involving a broad range of stress singularity power can be achieved numerically with the AT 1 model, provided that the internal length be computed using (2.6), which involves only material properties. In other words, our approach is capable of predicting crack nucleation near a weak stress singularity using only elastic properties, fracture toughness G c , the tensile strength σ c , and the local energy minimization principle (2.3).

In light of Figure 2.11, we suggest that both toughness and tensile strength (or equivalently toughness and internal length) can be measured by matching full domain or Pac-Man computations and experiments involving weak elastic singularity of various power (TPB, FPB, SENT, DENT with varying notch depth or angle) instead of measuring σ c directly. We expect that this approach will be much less sensitive to imperfections than the direct measurement of tensile strength, which is virtually impossible. Furthermore, since our criterion is not based on crack tip asymptotics, using full domain computations do not require that the experiments be specially designed to isolate the notch tip singularity from structural scale deformations. PMMA TPB (Dunn)

Al 2 O 3 -7%ZrO 2 FPB (Yoshibash) Al 2 O 3 -7%ZrO 2 TPB (Yoshibash) AT 1 -U FPB AT 1 -U TPB
a h = .1 PMMA TPB (Dunn) a h = .2 PMMA TPB (Dunn) a h = .3 PMMA TPB (Dunn) a h = .4 AT 1 -U a h = .1 AT 1 -U a h = .2 AT 1 -U a h = .3 AT 1 -U a h = .4 Figure

Initiation near a stress concentration: the U-notch

Crack nucleation in a U-shaped notch is another classical problem that has attracted a wealth of experimental and theoretical work. Consider a U-shaped notch of width ρ and length a ρ subject to a mode-I local loading (see Figure 2.12 for a description of notch geometry in the context of a double edge notch tension sample). Assuming "smooth" loadings and applied boundary displacements, elliptic regularity mandates that the stress field be non-singular near the notch tip, provided that ρ > 0. Within the realm of Griffith fracture, this of course makes crack nucleation impossible. As it is the case for the Vnotch, introducing a nucleation principle based on a critical stress is also not satisfying as it will lead to a nucleation load going to 0 as ρ → 0, instead of converging to that of an infinitely thin crack given by Griffith's criterion. There is a significant body of literature on "notch mechanics", seeking to address this problem introducing stress based criteria, generalized stress intensity factors, or intrinsic material length and cohesive zones. A survey of such models, compared with experiments on a wide range of brittle materials is given [START_REF] Gómez | Failure criteria for linear elastic materials with U-notches[END_REF].

In what follows, we study crack nucleation near stress concentrations in the AT 1 and AT 2 models and compare with the experiments gathered in [START_REF] Gómez | Failure criteria for linear elastic materials with U-notches[END_REF]. The core of their analysis consist in defining a generalized stress intensity factor

K U = K t σ ∞ c πρ 4 , (2.13)
where K t , the notch stress concentration factor, is a parameter depending on the local (a and ρ), as well as global sample geometry and loading. Through a dimensional analysis, they studied the dependence of the critical generalized stress intensity factor at the onset Chapter 2. Crack nucleation in variational phase-field models of brittle fracture 25, and 0.5 for which the value K t , computed in [START_REF] Lazzarin | A generalized stress intensity factor to be applied to rounded v-shaped notches[END_REF] is respectively 5.33, 7.26, and 11.12. In each case, we leveraged the symmetries of the problem by performing computations with the AT 1 and AT 2 models on a quarter of the domain for a number of values of the internal length corresponding to ρ/ ch between 0.05 and 20. In all cases, undamaged notch boundary conditions were used.

In Figure 2.13, we overlay the outcome of our simulations over the experimental results gathered in [START_REF] Gómez | Failure criteria for linear elastic materials with U-notches[END_REF]. As for the V-notch, we observe that the AT 2 model performs poorly for weak stress concentrations (large values of ρ/ ch ), as the lack of an elastic phase leads to the creation of large partially damaged areas. For sharp notches (ρ 0), our simulations concur with the experiments in predicting crack nucleation when K U = K Ic . As seen earlier, the AT 1 slightly overestimates the critical load in this regime when undamaged notch boundary conditions are used. In light of Figure 2.13, we claim that numerical simulations based on the variational phase-field model AT 1 provides a simple way to predict crack nucleation that does not require the computation of a notch stress concentration factors K t or the introduction of an ad-hoc criterion.

Size effects in variational phase-field models

Variational phase-field models are characterized by the intrinsic length , or ch . In this section, we show that this length-scale introduces physically pertinent scale effects, corroborating its interpretation as a material length. To this end, we study the nucleation of a crack in the uniaxial traction of a plate (-W, W ) × (-L, L) with a centered elliptical hole with semi-axes a and ρa (0 ≤ ρ ≤ 1) along the x-and y-axes respectively, see Figure 2.14. In Section 2.3.1, we study the effect of the size and shape of the cavity, assumed to be small with respect to the dimension of the plate (a W, L). In Section 2.3.2, we investigate material and structural size effects for a plate of finite width in the limit case of a perfect crack (ρ = 0). For a small hole (a W, L), up to a change of scale, the problem can be fully characterized by two dimensionless parameters: a/ , and ρ. For a linear elastic and isotropic material occupying an infinite domain, a close form expression of the stress field as a function of the hole size and aspect ratio is given in [START_REF] Inglis | Stresses in plates due to the presence of cracks and sharp corners[END_REF]. The stress is maximum at the points A = (a, 0) and A = (-a, 0), where the radial stress is zero and the hoop stress is given by:

σ max = t 1 + 2 ρ , (2.14)
t denoting the applied tensile stress along the upper and lower edges of the domain, i.e. the applied macroscopic stress at infinity. We denote by ū the corresponding displacement field for t = 1, which is given in [START_REF] Gao | A general solution of an infinite elastic plate with an elliptic hole under biaxial loading[END_REF].

As for the case of a perfect bar, (2.14) exposes a fundamental issue: if ρ > 0, the stress remains finite, so that Griffith-based theories will only predict crack nucleation if ρ = 0. In that case the limit load given by the Griffith's criterion for crack nucleation is

t = σ G := G c E aπ . (2.15)
However, as ρ → 0, the stress becomes singular so that the critical tensile stress σ c is exceeded for an infinitesimally small macroscopic stress t.

Following the findings of the previous sections, we focus our attention on the AT 1 model only, and present numerical simulations assuming a Poisson ratio ν = 0.3 and plane-stress conditions. We perform our simulations in domain of finite size, here a disk of radius R centered around the defect. Along the outer perimeter of the domain, we apply a boundary displacement u = tū, where ū is as in [START_REF] Inglis | Stresses in plates due to the presence of cracks and sharp corners[END_REF], and we use the macroscopic stress t a loading parameter. Assuming a symmetric solution, we perform our computations on a quarter domain. For the circular case ρ = 1, we use a reference mesh size h = min /10, where min is the smallest value of the internal length of the set of simulations. For ρ < 1, we selectively refine the element size near the expected nucleation site (see Figure 2.14-right). In order to minimize the effect of the finite size of the domain, we set R = 100a.

We performed numerical simulations varying the aspect ratio a/ from 0.1 to 50 and the ellipticity ρ from 0.1 to 1.0. In each case, we started from an undamaged state an monotonically increased the loading. In all numerical simulations, we observe two critical loading t e and t c , the elastic limit and structural strength, respectively. For 0 ≤ t < t e the solution is purely elastic, i.e. the damage field α remains identically 0 (see Figure 2.15left). For t e ≤ t < t c , partial distributed damage is observed. The damage field takes its maximum value α max < 1 near point A (see Figure 2.15-center). At t = t c , a fully developed crack nucleates, then propagates for t > t c (see Figure 2.15-right). As for the Pac-Man problem, we identify the crack nucleation with a jump in surface energy, and focus on loading at the onset of damage. From the one-dimensional problem of Section 2.1.2 and [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], we expect damage nucleation to take place when the maximum stress σ max reaches the nominal material strength σ c = 3G c E /8 (see Table 2.1), i.e. for a critical load

t e = ρ 2 + ρ σ c = ρ 2 + ρ 3G c E 8 . (2.16)
Figure 2.16-left confirms this expectation by comparing the ratio t e /σ c to its expected value ρ/(2 + ρ) for ρ ranging from 0.1 to 1. Figure 2.16-right highlights the absence of size effect on the damage nucleation load, by comparing t e /σ c for multiple values of a/ while keeping ρ fixed at 0.1 and 1.

Figure 2.17 focuses on the crack nucleation load t c , showing its dependence on the defect shape (left) and size (right). Figure 2.17-right shows the case of circular hole (ρ = 1) and an elongated ellipse, which can be identified to a crack (ρ = 0.1). It clearly highlights a scale effect including three regimes: i. For "small" holes (a ), crack nucleation takes place when t = σ c , as in the uniaxial traction of a perfect bar without the hole: the hole has virtually no effect on crack nucleation. In this regime the strength of a structure is completely determined by that of the constitutive material. Defects of this size do not reduce the structural strength and can be ignored at the macroscopic level.

ii. Holes with length of the order of the internal length (a = O( )), have a strong impact on the structural strength. In this regime the structural strength can be approximated by

log(t c /σ c ) = D log(a/ ) + c, (2.17)
where D is an dimensionless coefficient depending on the defect shape. For a circular hole ρ = 1, we have D ≈ -1/3.

iii. When a , the structural failure is completely determined by the stress distribution surrounding the defect. We observe that for weak stress singularities (ρ ≡ 1), nucleation takes place when the maximum stress reaches the elastic limit σ e , whereas the behavior as ρ ≡ 0 is consistent with Griffith criterion, i.e. the nucleation load scales as 1/ √ a.

Figure 2.17-right shows that the shape of the cavity has a significant influence on the critical load only in the latter regime, a

. Indeed, for a/ of the order of the unity or smaller, the critical loads t c for circular and highly elongated cavities are almost indistinguishable. This small sensitivity of the critical load on the shape is the result of the stress-smoothing effect of the damage field, which is characterized by a cut-off length of the order of . Figure 2.17-left shows the critical stress t c at nucleation when varying the aspect ratio ρ for a/ = 48, for which σ G /σ c = 2/15. As expected, the critical stress varies smoothly from the value σ G (2.15) predicted by the Griffith theory for a highly elongated cavity identified to a perfect crack, to t e (2.16) for circular cracks, where the crack nucleates as soon as the maximum stress σ max attains the elastic limit.

This series of experiments is consistent with the results of Section 2.2.2 showing that variational phase-field models are capable of simultaneously accounting for critical elastic energy release rate and critical stress. Furthermore, they illustrate how the internal length can be linked to critical defect size as the nucleation load for a vanishing defect of size less than approaches that of a flawless structure.

Competition between material and structural size effects

We can finally conclude the study of size effects in variational phase-field models by focusing on the competition between material and structural size effects. For that matter, we study the limit case ρ = 0 of a perfect crack of finite length 2a in a plate of finite width 2W (see Figure 2.18-left). Under the hypotheses of LEFM, the critical load upon which the crack propagates is

σ G (a/ ch , a/W ) = G c E cos( aπ 2W ) aπ = σ c 1 π ch a cos aπ 2W , (2.18)
which reduces to (2.15) for large plate (W/a → ∞). As before, we note that σ G /σ c → ∞ as a/ ch → 0, so that for any given load, the material's tensile strength is exceeded for short enough initial crack. We performed series of numerical simulations using the AT 1 model on a quarter of the domain with W = 1, L = 4, ν = 0.3, = W/25, h = /20, and the initial crack's halflength a ranging from from 0.025 to 12.5 (i.e. 0.001W to 0.5W ). The pre-existing crack was modeled as a geometric feature and undamaged crack lip boundary conditions were prescribed. The loading was applied by imposing a uniform normal stress of amplitude t to its upper and lower edge. theories linking size-effect on the strength of the material [START_REF] Bažant | Scaling of Structural Strength[END_REF]. When a , i.e. when the defect is large compared to the material's length, crack initiation is governed by Griffith's criterion (2.18). As noted earlier, the choice of undamaged notch boundary conditions on the damage fields leads to slightly overestimating the nucleation load. Our numerical simulations reproduce the structural size effect predicted by LEFM when the crack length is comparable to the plate width W .

When a , we observe that the macroscopic structural strength is very close to the material's tensile strength. Again, below the material's internal length, defects have virtually no impact on the structural response. LEFM and Griffith-based models cannot account for this material size-effect. These effects are introduced in variational phase-field model by the additional material parameter .

In the intermediate regime a = O( ), we observe a smooth transition between strength and toughness criteria, where the tensile strength is never exceeded. When a , our numerical simulations are consistent with predictions from Linear Elastic Fracture Mechanics shown as a dashed line in Figure 2.18, whereas when a , the structural effect of the small crack disappear, and nucleation takes place at or near the material's tensile strength, i.e. t c /σ c 1.

Conclusion

In contrast with most of the literature on phase-field models of fracture focusing validation and verification in the context of propagation "macroscopic" cracks [START_REF] Mesgarnejad | Validation simulations for the variational approach to fracture[END_REF][START_REF] Pham | Experimental validation of a phase-field model for fracture[END_REF], we have studied crack nucleation and initiation in multiple geometries. We confirmed observations reported elsewhere in the literature that although they are mathematically equivalent in the limit of → 0, damaged notch boundary conditions lead to a more accurate computation near strong stress singularities whereas away from singularities, undamaged notch boundary conditions are to be used. Our numerical simulations also highlight the superiority of phase-field models such as AT 1 which exhibit an elastic phase in the one-dimensional tension problem over those who don't (such as AT 2 ), when nucleation away from strong singularity is involved. Our numerical simulations suggest that it is not possible to accurately account for crack nucleation near "weak" singularities using the AT 2 model. We infer that a strictly positive elastic limit σ e is a required feature of a phase-field model that properly account for crack nucleation.

We have shown that as suggested by the one-dimensional tension problem, the regularization parameter must be understood (up to a model-dependent multiplicative constant) as the material's characteristic or internal length ch = G c E/σ 2 c , and linked to the material strength σ c . With this adjustment, we show that variational phasefield models are capable of quantitative prediction of crack nucleation in a wide range of geometries including three-and four-point bending with various type of notches, single and double edge notch tests, and a range of brittle materials, including steel and Duraluminium at low temperatures, PVC foams, PMMA, and several ceramics.

We recognize that measuring a material's tensile strength is difficult and sensitive to the presence of defect, so that formulas (2.6) may not be a practical way of computing a material's internal length. Instead, we propose to perform a series of experiments such as three point bending with varying notch depth, radius or angle, as we have demonstrated in Figure 2.11 that with a properly adjusted internal length, variational phase-field models are capable of predicting the nucleation load for any notch depth or aperture. Furthermore, since variational phase-field models do not rely on any crack-tip asymptotic, this identification can be made even in a situation where generalized stress or notch intensity factors are not known or are affected by the sample's structural geometry.

We have also shown that variational phase-field models properly account for size effects that cannot be recovered from Griffith-based theories. By introducing the material's internal length, they can account for the vanishing effect of small defects on the structural response of a material, or reconcile the existence of a critical material strength with the existence of stress singularity. Most importantly, they do not require introducing ad-hoc criteria based on local geometry and loading. On the contrary, we see that in most situation, criteria derived from the asymptotic analysis of a micro-geometry can be recovered a posteriori. Furthermore, variational phase-field models are capable of quantitative prediction of crack path after nucleation. Again, they do so without resolving to introduce additional ad-hoc criteria, but only rely on a general energy minimization principle.

In short, we have demonstrated that variational phase-field models address some of the most vexing issues associated with brittle fracture: scale effects, nucleation, existence of a critical stress, and path prediction.

Of course, there are still remaining issues that need to be addressed. Whereas the models are derived from irreversibility, stability and energy balance, our numerical simulations do not enforce energy balance as indicated by a drop of the total energy upon crack nucleation without string singularities. Note that to this day, devising an evolution principle combining the strength of (2.3) while ensuring energy balance is still an open

Appendix B

Tables of experimental an numerical data for V-notch experiments [START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] using three point bending experiments of a PMMA sample compared to full domain numerical simulations using the AT 1 model with undamaged notch boundary conditions. The value a/h refers to the ratio depth of the notch over sample thickness. See Figure 2.9 for geometry and loading.

ω λ k c k c k c k c (AT 1 -U) (AT 1 -D) (AT 2 -U) (AT 2 -D) 0 01°0.
Chapter 3

A phase-field model for hydraulic fracturing in low permeability reservoirs: propagation of stable fractures

Hydraulic fracturing is a process to initiate and to extend fractures by injecting fluid into subsurface. Mathematical modeling of hydraulic fracturing requires coupled solution of models for fluid flows and reservoir-fracture deformation. The governing equations for these processes are fairly well understood and includes, for example, the Reynold's equation, cubic law, diffusivity equation and Darcy's law for fluid flow modeling, linear poro-elasticity equation for reservoir-fracture deformation and Griffith's criterion for fracture propagation. Considering that fracture propagation is a moving boundary problem, the numerical and computational challenges of solving these governing equations on the fracture domain limit the ability to comprehensively model hydraulic fracturing. These challenges include but are not limited to, finding efficient ways of representing numerically the fracture and reservoir domains in the same computational framework while still ensuring hydraulic and mechanical coupling between both subdomains. To address these issues, several authors have assumed a known propagation path that is limited to a coordinate direction of the computational grid [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF][START_REF] Boone | A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media[END_REF] while some others simply treated fractures as external boundaries of the reservoir computational domain [START_REF] Ji | A novel hydraulic fracturing model fully coupled with geomechanics and reservoir simulation[END_REF][START_REF] Dean | Hydraulic-fracture predictions with a fully coupled geomechanical reservoir simulator[END_REF]. Special interface elements called zero-thickness elements have also been used to handle fluid flow in fractures embedded in continuum media [START_REF] Carrier | Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model[END_REF][START_REF] Segura | On zero-thickness interface elements for diffusion problems[END_REF][START_REF] Segura | Coupled hm analysis using zero-thickness interface elements with double nodes. part I: Theoretical model[END_REF][START_REF] Segura | Coupled hm analysis using zero-thickness interface elements with double nodes part II: Verification and application[END_REF][START_REF] Boone | A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media[END_REF][START_REF] Lobão | Modelling of hydrofracture flow in porous media[END_REF]. Despite the simplicity of these approaches and contrary to field evidence of complex fracture geometry and propagation paths, they have limited ability to reproduce realistic fracture behaviors.

Where attempts have been made to represent fractures and reservoir in the same computational domain, for instance using the extended finite element method (XFEM) [START_REF] Mohammadnejad | An extended finite element method for hydraulic fracture propagation in deformable porous media with the cohesive crack model[END_REF][START_REF] Dahi | Analysis of hydraulic fracture propagation in fractured reservoirs: an improved model for the interaction between induced and natural fractures[END_REF] and the generalized finite element method (GFEM) [START_REF] Gupta | Simulation of non-planar three-dimensional hydraulic fracture propagation[END_REF], the computational cost is high and the numerics cumbersome, characterized by continuous remeshing to provide grids 3.1. A phase fields model for hydraulic fracturing that explicitly match the evolving fracture surface. Some of these challenges can be overcome using a phase field representation for fractures as evident in the work of [START_REF] Bourdin | A variational approach to the numerical simulation of hydraulic fracturing[END_REF] and [START_REF] Bourdin | A variational approach to the modeling and numerical simulation of hydraulic fracturing under in-situ stresses[END_REF]. This chapter extends the works of [START_REF] Bourdin | A variational approach to the numerical simulation of hydraulic fracturing[END_REF] by applying the variational phase field model to a network of fractures. The hydraulic fracture model is developed by incorporating fracturing fluid pressure in Francfort and Marigo's variational approach to fracture [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Specifically, the fracture model recast Griffith's propagation criteria into a total energy minimization problem, where the global energy is the sum of the elastic and fracture surface energies, the fracturing fluid pressure force and the work done by in-situ stresses. We assume quasi static fracture propagation and in this setting, the fractured state of the reservoir is the solution of a series of minimizations of this total energy with respect to any kinematically admissible crack sets and displacement field. Numerical implementation of the model is based on a phase field representation of the fracture and subsequent regularization of the total energy function. The phase field technique avoids the need for explicit knowledge of fracture location, it permits the use of a single computational domain for fracture and reservoir representation. The strength of this method is to provide a unified setting for handling path determination, nucleation and growth of arbitrary number of stable cracks in any dimensions based on the energy minimization principle. This work focuses on the fracture propagation stability through various examples such as, a pressurized single fracture stimulated by a controlled injected volume in a large domain, a network of multiple parallel fractures and a pressure driven laboratory experiment to measure rocks toughness.

The Chapter is organized as follows: Section 3.1 is devoted to recall phase field models for hydraulic fracturing in the toughness dominated regime with no fluid loss to the impermeable elastic reservoir [START_REF] Detournay | The near tip region of a fluid driven fracture propagating in a permeable elastic solid[END_REF]. Then, our numerical implementation scheme and algorithm for volume driven hydraulic fracturing simulations is exposed in section 3.1.3. Tough the toughness dominated regime may not cover the whole spectrum of fracture propagation but provides an appropriate framework for verifications since it does not require the solution of a flow model. Therein, section 3.2 is concerned with comparisons between our numerical results and the closed form solutions provided by Sneddon [START_REF] Sneddon | The opening of a griffith crack under internal pressure[END_REF][START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] for the fluid pressure, fracture length/radius and fracture volume in a single crack case. Section 3.3 focuses on the propagation of infinite pressurized parallel fractures and it is compared with the derived solution. Section 3.4 is devoted to study the pre-fracture stability in the burst experiment at a controlled pressure. This test proposed by Abou-Sayed [START_REF] Abou-Sayed | An Experimental Technique for Measuring the Fracture Toughness of Rocks under Downhole Stress Conditions[END_REF] is designed to measure the fracture toughness of the rock and replicates situations encountered downhole with a borehole and bi-wing fracture.

A phase fields model for hydraulic fracturing

A variational model of fracture in a poroelastic medium

Consider a reservoir consisting of a perfectly brittle isotropic homogeneous linear poroelastic material with A the Hooke's law tensor and G c the critical energy release rate Chapter 3. A phase-field model for hydraulic fracturing in low permeability reservoirs: propagation of stable fractures occupying a domain Ω ⊂ R n , n = 2 or 3 in its reference configuration. The domain is partially cut by a sufficiently regular crack set Γ ⊂ Ω with Γ ∩ ∂Ω = ∅. A uniform pressure denoted by p applies on both faces of the fracture lips i.e. Γ = Γ + ∪ Γ -and pore pressure denoted by p p applies in the porous material which follows the Biot poroelastic coefficient λ. The sound region Ω \ Γ is subject to a time independent boundary displacement ū(t) = 0 on the Dirichlet part of its boundary ∂ D Ω and time stress dependent g(t) = σ • ν on the remainder ∂ N Ω = ∂Ω \ ∂ D Ω, where ν denotes the appropriate normal vector. For the sake of simplicity body forces are neglected such that at the equilibrium, the stress satisfies, div σ = 0

where the Cauchy stress tensor follows Biot's theory [START_REF] Biot | General theory of three-dimensional consolidation[END_REF], i.e. σ = σλp p I, σ being the effective stress tensor. The infinitesimal total deformation e(u) is the symmetrical part of the spatial gradient of the displacement field u,

e(u) = ∇u + ∇ T u 2 .
The stress-strain relation is σ = Ae(u), so that,

σ = A e(u) - λ 3κ p p I ,
where 3κ is the material's bulk modulus. Those equations can be rewritten in a variational form, by multiplying the equilibrium by the virtual displacement v ∈ H 1 0 (Ω \ Γ; R n ) and using Green's formula over Ω \ Γ. After calculation, we get that,

Ω\Γ σ : e(v) dx - ∂ N Ω g(t) • v dH n-1 - Γ p v • ν dH n-1 = 0 (3.1)
where H n-1 denotes the n -1-dimensional Hausdorff measure, i.e. its aggregate length in 2 dimensions and surface area in 3 dimensions. Finally, we remark that the above equation (3.1) can be seen as the Euler-Lagrange equation for the minimization of the elastic energy,

E(u, Γ) = Ω\Γ 1 2 A e(u) - λ 3κ p p I : e(u) - λ 3κ p p I dx - ∂ N Ω g(t) • u dH n-1 - Γ p u • ν dH n-1 (3.2) amongst all displacement fields u ∈ H 1 (Ω \ Γ; R n ) such that u = 0 on ∂ D Ω.

A phase fields model for hydraulic fracturing

Remark 2 Of course, fluid equilibrium mandates continuity of pressure so that p p = p along Γ. Our choice to introduce two pressure fields is motivated by our focus on lowpermeability reservoirs. In this situation, assuming very small leak-off, it is reasonable to assume that for short injection time, the pore pressure is "almost" constant away from the crack, hence that p = p p .

We follow the formalism of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF] and propose a time-discrete variational model of crack propagation. To any crack set Γ ⊂ Ω and any kinematically admissible displacement field u, we associate the fracture energy,

E(u, Γ) = Ω\Γ 1 2 A e(u) - λ 3κ p p I : e(u) - λ 3κ p p I dx - ∂ N Ω g(t) • u dH n-1 - Γ p u • ν dH n-1 + G c H n-1 (Γ) (3.3)
Considering then a time interval [0, T ] and a discrete set of time steps 0 = t 0 < t 1 < • • • < t N = T , and denoting by p i , p p i and g i , the crack pressure, pore pressure and external stress at time t i (i > 0), we postulate that the displacement and crack set (u i , Γ i ) are minimizers of E amongst all kinematically admissible displacement fields u and all crack sets Γ satisfying a growth condition Γ j ⊂ Γ for all j < i, with Γ 0 possibly representing pre-existing cracks.

It is worth emphasizing that in this model, no assumptions are made on the crack geometry Γ i . As in Francfort and Marigo's pioneering work [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], minimization of the total fracture energy is all that is needed to fully identify the crack geometry (path) and topology (nucleation, merging, branching).

Variational phase-field approximation

Several techniques have been proposed for the numerical implementation of the fracture energy E, the main difficulty being to handle discontinuous displacements along unknown surfaces. In recent years, variational phase-field models, originally devised in [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF][START_REF] Ambrosio | Approximation of functionals depending on jumps by elliptic functionals via Γ-convergence[END_REF], and extended to brittle fracture [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] have become very popular.

We follow this approach by introducing a regularization length , an auxiliary field α with values in [0, 1] representing the unknown crack surface, and the regularized energy.

E (u, α) = Ω 1 2 A (1 -α)e(u) - λ 3κ p p I : (1 -α)e(u) - λ 3κ p p I dx - ∂ N Ω g(t) • u dH n-1 + Ω pu • ∇α dx + 3G c 8 Ω α + |∇α| 2 dx (3.4)
where α = 0 is the undamaged state material and α = 1 refers to the broken part. One can recognize the AT 1 model introduced in the Chapter 1 which differs from one used in [START_REF] Chukwudozie | Application of the Variational Fracture Model to Hydraulic Fracturing in Poroelastic Media[END_REF].

Chapter 3. A phase-field model for hydraulic fracturing in low permeability reservoirs: propagation of stable fractures

At each time step, the constrained minimization of the fracture energy E is then replaced with that of E , with respect to all (u i , α i ) such that u i is kinematically admissible and 0 ≤ α i-1 ≤ α i ≤ 1.

The Γ-convergence of (3.4) to (3.3), which constitutes the main justification of variational phase-field models is a straightforward extension of [START_REF] Chambolle | An approximation result for special functions with bounded variations[END_REF][START_REF] Chambolle | Addendum to "An Approximation Result for Special Functions with Bounded Deformation[END_REF], or [START_REF] Iurlano | A density result for gsbd and its application to the approximation of brittle fracture energies[END_REF]. It is quite technical and not quoted here. The form of the regularization of the surface energy in (3.4) is slightly different from the one originally proposed in [START_REF] Bourdin | A variational approach to the modeling and numerical simulation of hydraulic fracturing under in-situ stresses[END_REF][START_REF] Bourdin | A variational approach to the numerical simulation of hydraulic fracturing[END_REF] but this choice is motivated by the work of [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF][START_REF] Bourdin | Morphogenesis and propagation of complex cracks induced by thermal shocks[END_REF].

In the context of poro-elasticity, the regularization of the elastic energy of the form of,

Ω 1 2 A (1 -α)e(u) - λ 3κ p p I : (1 -α)e(u) - λ 3κ p p I dx
is different from that of [START_REF] Mikelic | A quasistatic phase field approach to fluid filled fractures[END_REF] and follow-up work, or [START_REF] Miehe | Minimization principles for the coupled problem of Darcy-Biot-type fluid transport in porous media linked to phase field modeling of fracture[END_REF][START_REF] Wilson | A phase-field model for fracture in piezoelectric ceramics[END_REF] which use a regularization of the form

Ω 1 2 (1 -α) 2 A e(u) - λ 3κ p p I : e(u) - λ 3κ p p I dx.
This choice is consistent with the point of view that damage takes place at the sub-pore scale, so that the damage variable α should impact the Cauchy stress and not the effective poro-elastic stress. Note that as → 0, both expressions will satisfy Γ-convergence to E.

A fundamental requirement of hydraulic fracturing modeling is volume conservation, that is the sum of the fracture volume and fluid lost to the surrounding reservoir must equal the amount of fluid injected denoted V . In the K-regime, the injected fluid is inviscid and no mass is transported since the reservoir is impermeable. Of course, reservoir impermeability means no fluid loss from fracture to reservoir and this lack of hydraulic communication means that the reservoir pressure p p and fracture fluid pressure p are two distinct and discontinuous quantities. Furthermore, the zero viscosity of the injected fluid is incompatible with any fluid flow model, leaving global volume balance as the requirement for computing the unknown fracturing fluid pressure p. In the sequel we set aside the reservoir pressure p p and consider this as a hydrostatic stress offset in the domain, which can be recast by applying a constant pressure on the entire boundary of the domain.

Numerical implementation

The numerical implementation of the variational phase-field model is well established. In the numerical simulations presented below, we discretized the regularized fracture energy using linear or bilinear finite elements. We follow the classical alternate minimizations approach of [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and adapt to volume-driven fractures where main steps are: i. For a given (α, p) the minimization of E with respect to u is an elastic problem with the prescribed boundary condition. To solve this, we employed preconditioned conjugate gradient methods solvers.

3.2. Numerical verification case of a pressurized single fracture in a two and three dimensions

ii. The minimization of E with respect to α for fixed (u, p) and subject to irreversibility (α ≥ α i-1 ) is solved using variational inequality solvers provided by PETCs [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | PETSc Web page[END_REF].

iii. For a fixed (u, α), the total volume of fluid can be computed, such that,

V = - Ω u • ∇α dx.
The idea is to rescale the fluid pressure using the secant method (a root-finding algorithm) based on a recurrence relation.

A possible algorithm to solve volume-driven hydraulic fracturing is to use nested loops. The inner loop solves the elastic problem i. and rescale the pressure iii. until the error between the target and the computed volume is below a fixed tolerance. The outer loop is composed of ii. and the previous procedure and the exit is triggered once the damage has converged. This leads to the following Algorithm 2 where δ V and δ α are fixed tolerances. Remark that the inner loop solves a linear problem, hence, finding the pressure p associated to the target volume V should converge in strictly less than four iterations. All computations were performed using the open source mef90 1 .

In-situ stresses play a huge role in hydraulic fracture propagation and the ability to incorporate them in a numerical model is an important requirement for robust hydraulic fracturing modeling. Our numerical model easily accounts for these compressive stresses on boundaries of the reservoir. However in-situ stresses simulated cannot exceeded the maximum admissible stress of the material given by σ c = 3EG c /8 . We run a series of two-and three-dimensions computations to verify our numerical model and investigate stability of fractures.

Numerical verification case of a pressurized single fracture in a two and three dimensions

Using the Algorithm 2 a pressurized line and penny shape fractures have been respectively simulated in two-and three-dimensions, and their results compared with the closed form solutions. Both problems have a symmetric axis, i.e. its aggregate a reflexion axis in 2d and a rotation in 3d, leading to a invariant geometry drawn on Figure 3.1. Also, all geometric and material parameters are identically set up for both problems and summarized in the Table 3.1. The closed form solutions provided by Sneddon in [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF][START_REF] Sneddon | The opening of a griffith crack under internal pressure[END_REF] are recalled in the Appendix C and assume an infinite domain with vanishing stress and displacement at the boundary. To satisfy those boundary conditions we performed simulations on a huge domain clamped at the boundary, where the reservoir size is 100 times larger than the pre-fracture length as reported in the Table 3.1. To moderate the number of elements in the domain, a casing (W, H) with a constant refined mesh size of resolution h is encapsulated around the fracture. Outside the casing a coarsen mesh is spread out see Figure 3.1.

Chapter 3. A phase-field model for hydraulic fracturing in low permeability reservoirs: propagation of stable fractures Algorithm 2 Volume driven hydraulic fracturing algorithm at the step i 1: Let j = 0 and α 0 := α i-1 2: repeat 3:

Set, p k-1 i = p k i and V k-1 i = V k i 4:
p k+1 i := p k i -V k i (p k i -p k-1 i )/(V k i -V k-1 i ) 6:
Compute the equilibrium,

u k+1 := argmin u∈C i E (u, α j ) 7:
Compute volume of fractures,

V k+1 i := - Ω u k+1 • ∇α j dx 8: k := k + 1 9: until V k i -V i L ∞ ≤ δ V 10:
Compute the damage,

α j+1 := argmin α∈D i α≥α i-1 E (u j+1 , α)
11:

j := j + 1 12: until α j -α j-1
L ∞ ≤ δ α 13: Set, u i := u j and α i := α j refined mesh, size h coarsen mesh symmetry axis 3.1: Parameters used for the simulation of a single fracture in two and three dimensions.

A loading cycle is preformed by pressurizing the fracture until propagation, then, pumping all the fluid out of the crack. The pre-fracture of length l 0 is measured by a isovalues contour plot for α = .8 before refilling the fracture of fluid again. The reason of this is we do not have an optimal damage profile at the fracture tips, leading to underestimate the critical pressure p c . Similar issues have been observed during the nucleation process in [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] where G c is overshoot due to the solution stability. Snap-shots of the damage before and after the loading cycle in the Figure 3.3 illustrate differences between damage profiles at the crack tips. Since the critical crack pressure is a decreasing function with respect to the crack length, the maximum value is obtained at the loading point when the crack initiates (for the pre-fracture). One can see on the Figure 3.4 that the penny shape fracture growth is not necessary symmetrical with respect to the geometry but remains a disk shape which is consistent with the invariant closed form solution.

We know from prior work see [START_REF] Bourdin | The variational approach to fracture[END_REF] that the "effective" numerical toughness is quantified by (G c ) eff = G c (1 + 3h/(8 ) ) in two dimensions. However, for the penny shape crack (G c ) eff = G c (1 + 3h/(8 ) + 2h/l ), where 2h is the thickness of the crack and l the radius. The additional term of 2h/l comes from the lateral surface contribution which becomes negligible for thin fractures.

The fluid pressure p and the fracture length l closed form solution with respect to the total injected volume of fluid V is provided by [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF] and is recalled in the Appendix C. Figure 3.2 shows a perfect match between the numerical results and the closed solution for the line fracture and penny shape crack. In both cases as long as the V ≤ V c the crack does not grow, and since V > V c the pressure drop as p ∼ V -1/3 (line fracture) and p ∼ V -1/5 (penny shape crack). Notice that the pressure decreases when the crack grows, therein a pressure driven crack is necessary unstable, indeed there is no admissible pressure over the maximum value p c .

Remark 3 The Griffith regime requires σ

c = 3E G c /(8 ) ≥ πE G c /(4l) = p c
in two dimensions, leading to l ≥ 2π /3. Therefore, the pre-fracture must be longer than twice the material internal length to avoid any size effects phenomena as reported in Chapter 2.

Those simulations show that the variational phase field model to hydraulic fracturing recovers Griffith's initiation and propagation for a single pressurized crack. Even if this can be seen as a toy example because the fracture propagation is rectilinear, without any changes on the implementation multi-fracking can be simulated as illustrated in the Figure 3.5. Fracture paths are obtained by total energy minimization and satisfies Griffith's propagation criterion.

Multi fractures in two dimensions

Multi fractures in two dimensions

One of the most important features of our phase field hydraulic fracturing model is its ability to handle multiple fractures without additional computational or modeling effort than is required for simulating single fracture. This capability is highlighted in the following study of the stimulation of a network of parallel fractures. All cracks are subject to the same pressure and we control the total amount of fluid injected into cracks, i.e. fluid can migrate from a crack to another via a wellbore.

The case where all fractures of a parallel network propagate (multi fracking scenario) is often postulated. However, considering the variational structure of Griffith leads to a different conclusion. For the sake of simplicity consider only two parallel fractures. A virtual extension of one of the cracks (variational argument) induces a drop of pressure in both fractures. Consequently the shorter fracture is sub-critical and remains unchanged since the pressure p < p c . Moreover the longer fracture requires less pressure to propagate than the shorter because the critical pressure decreases with the crack length. Finally the longer crack continues to propagate. This non restrictive situation can be extended to multiple fractures (parallel and the same size). In the sequel, we propose to revisit the hypothesis of multi-fracking by performing numerical simulations using the Algorithm 2. Consider a network of infinite parallel cracks with the same pressure p where their individual length is l and the spacing between cracks is δ drawn in the Figure 3.6 (left). At the initial state all pre-cracks have the same length denoted l 0 and no in-situ stresses is applied on the reservoir domain.

Multi-fracking closed form solution

This network of parallel cracks is a duplication of an invariant geometry, precisely a strip domain Ω = (-∞, +∞) × [-δ, δ] cut in the middle by a fracture Γ = [-l, l] × {0}. An asymptotic solution of this cell domain problem is provided by Sneddon in [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF]

V (ρ) = 8pδ 2 E π ρ 2 f (ρ), (3.5)
where the density of fractures ρ = lπ/(2δ) and

f (ρ) = 1 -ρ 2 /2 + ρ 4 /3 + o(ρ 6 ).
The Taylor series of f (ρ) in 0 provided by Sneddon differs from one given in the reference [START_REF] Murakami | Handbook of stress intensity factors[END_REF] where f (ρ) = 1ρ 2 /2 + 3ρ 4 /8 + o(ρ 6 ). The latter is exactly the first three terms of the expansion of

f (ρ) = 1 1 + ρ 2 . (3.6)
The critical pressure satisfying Griffith propagation for this network of fractures problem is

p(ρ) = E G c δ(ρ 2 f (ρ)) (3.7)
Of course the closed form expression consider that all cracks grow by symmetry. It is convenient for numerical reason to consider an half domain and an half fracture (a crack lip) of the reference geometry such that we have (Ω 1 , Γ 1 ) and by symmetry expansion (Ω 2 , Γ 2 ), (Ω 4 , Γ 4 ),..,(Ω 2n , Γ 2n ) illustrated in the Figure 3.6 (right).

Numerical simulation of multi-fracking by computation of unit cells construction

The idea is to reproduce numerically multi-fracking scenario, thus simulation is performed on stripes of length 2L with pre-fractures of length 2l 0 such that, geometries considered are:

Ω 2n = [-L, L] × [0, (2n -2)δ] Γ 0, 2n = [-l 0 , l 0 ] × n k=1 {2(k -1)δ} (3.8)
for n ≥ 1, n being the number of crack lips. Naturally a crack is composed of two lips. The prescribed boundary displacement on the top-bottom extremities is u y (0) = u y (2(n -1)δ) = 0, and on the left-right is u(±L) = 0. All numerical parameters used are set up in the Table 3.2.

h

L δ l 0 E ν G c 0.005 10 1 0.115 1 0 1 3h Table 3.2: Parameters used in the numerical simulation for infinite cracks Using the same technique of loading cycle as in section 3.2 and after pumping enough fluid into the system of cracks we observed in all simulations performed that only one fracture grows, precisely the one at the boundary as illustrated in the Figure 3.7. By

Multi fractures in two dimensions

using reflexion symmetry we have a periodicity of large fractures of 1/n. We notice that simulations performed never stimulate middle fracture. Indeed, by doing so after reflexions this will lead to an higher periodicity cases. The pseudo-color blue is for undamaged material and turns white when (α ≤ .01) (visibility reason). The full colors correspond to numerical simulations cells domain see table 3.2, and opacity color refers to the rebuild solution using symmetries. In all simulations only one crack propagates in the domain. Using the multiplicity pictures from left to right we obtain a fracture propagation periodicity denoted period. of 6/6, 3/6, 1.5/6 and 1/6.

To compare the total injected fluid V between simulations, we introduce the fluid volume density i.e the fluid volume for a unit geometry cell given by 2V /n. The evolution of normalized pressure, volume of fluid per cell and length are plotted in Figure 3.8 and show that the multi-fracking situation (one periodic) match perfectly with the close form solution provided by the equations (3.7),(3.5) and (3.6). Also, one can see that Sneddon approximation is not accurate for dense fractures. We can observe from simulations in Figure 3.8 that a lower periodicity (1/n) of growing cracks implies a reduction in pressure evolution. Also notice that the rate of pressure drop increases when the number of long cracks decrease, so that rapid pressure drop may indicate a poor stimulation. Also this loss of multi fracking stimulation decreases the fracture surface are for resource recovery. All cracks propagating simultaneously case is not stable in the sense that there exits a lower energy state with fewer growing crack. However as we will be discussed in the section 3.3.3 multi fracking may work for low fracture density since their interactions are negligible.

Chapter 3. A phase-field model for hydraulic fracturing in low permeability reservoirs: propagation of stable fractures

Multi-fracking for dense fractures

In the following we investigate critical pressure with respect to the density of fracture for different periodicity. Colored plot are numerical results for different domain sizes Ω 1 , Ω 2 , Ω 4 and Ω 6 . The solid black line is the closed form solution and gray the approximated solution given by Sneddon [START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF].

Let us focus on fractures propagation when their interactions become stronger i.e. higher fracture density ρ = lπ/(2δ). We start by normalizing the pressure relation for multi-fracking equation (3.7) with p = E G c /(lπ) which is a single fracture problem studied in section 3.2.

r p (ρ) = 2ρ (ρ 2 f (ρ)) = 2(ρ 2 +1) 3/2 ρ 2 +2
.

(3.9)

Remark that r p (0) = 1 means that critical pressure for largely spaced fractures are identical to a line fracture in a infinite domain problem, thus cracks behave without interacting each other. We run a set of numerical simulations using the same set of parameters than previously recalled in the Table 3.2 except that δ varies. For a high fractures density we 3.4. Fracture stability in the burst experiment with a confining pressure discovered another loss of symmetry shown on Figure 3.9 such that the fracture grows only in one direction. Figure 3.9: Domains in the deformed configuration for respectively Ω 2 and Ω 4 with 2δ = .5. The pseudo-color blue is for undamaged material and turns white when (α ≤ .01) (visibility reason). The full colors correspond to numerical simulation domain and opacity color refers to the rebuild solution using symmetries. In all simulations only one crack tip propagates in one direction in the simulated domain.

We report pressures obtained numerically depending on the fracture density in the Figure 3.10 and comparison with (3.9). One can see that the closed form solution is in good agreement with numerical simulation for the periodicity one and also lower periodicity obtained by doing ρ ← ρ/n in the equation (3.9). We see that for low fractures density the critical pressure is equal to a line fracture in a infinite domain. For higher fractures density, interactions become stronger and propagating all fractures require a high pressure compare to grow only one of them. As an example, a network of pre-fractures of length l = 6.36m and spaced δ = 10m thus ρ = 1, in this situation the required pressure is equal to r(1)K Ic / √ lπ with r(1) = 1.4 to propagate all cracks together compare to r(1) = 1 for only one single fracture. Naturally the system bifurcate to less fractures propagation leading to a drop of the fluid pressure.

Fracture stability in the burst experiment with a confining pressure

This section focuses on the stability of fractures propagation in the burst experiment. This laboratory experiment was conceived to measure the resistance to fracturing K Ic (also called the fracture toughness) of rock under confining pressure which is a critical parameter to match the breakdown pressure in mini-frac simulation. The idea is to provide a value of K Ic for hydraulic fracturing simulations in the K-regime [START_REF] Detournay | Propagation regimes of fluid-driven fractures in impermeable rocks[END_REF]. However, past experimental studies suggest that the fracture toughness of rock is dependent on the confining pressure under which the rock is imposed. Various methodologies exist for the Black dash line is r p (ρ/n) with 1/n the periodicity. Colored line is numerical results for respectively a periodicity 6/6, 3/6 and 1.5/6. measurement of K Ic under confining pressure and results differ in each study. The most accepted methodology in petroleum industry is the so called burst experiment, which was proposed by Abou-Sayed [START_REF] Abou-Sayed | An Experimental Technique for Measuring the Fracture Toughness of Rocks under Downhole Stress Conditions[END_REF], as the experimental geometry replicates a situation encountered downhole with a borehole and bi-wing fracture. Under linear elastic fracture mechanics, stable and unstable crack growth regime have been calculated depending on the confining pressure and geometry. During unstable crack propagation the phase-field models for hydraulic fracturing do not bring information. Instead we perform a Stress Intensity Factor (SIF) analysis along the fracture path to determine propagation stability regimes, herein this section is different from the phase-field sprite of the dissertation. However at the end we will verify the ability of the phase-field model to capture fracture stability transition from stable to unstable.

The burst experiment

The effect of confining pressure on the fracture toughness was first studied by Schmidt and Huddle [START_REF] Schmidt | Effect of Confining Pressure on Fracture Toughness of Indiana Limestone[END_REF] on Indiana limestone using single-edge-notch samples in a pressure vessel. In their experiments, increase in the fracture toughness up to four fold have been reported. Other investigations to quantify the confining pressure dependency were performed on the three point bending [START_REF] Müller | Brittle crack growth in rocks[END_REF][START_REF] Vásárhelyi | Influence of pressure on the crack propagation under mode i loading in anisotropic gneiss[END_REF], modified ring test [START_REF] Thiercelin | Fracture Toughness and Hydraulic Fracturing[END_REF], chevron notched Brazillian disk [START_REF] Roegiers | Rock fracture tests in simulated downhole conditions[END_REF], cylinder with a partially penetrating borehole [START_REF] Holder | Measurements of effective fracture toughness values for hydraulic fracture: Dependence on pressure and fluid rheology[END_REF][START_REF] Sitharam Thallak | The pressure dependence of apparent hydrofracture toughness[END_REF], and thick wall cylinder 3.4. Fracture stability in the burst experiment with a confining pressure with notches [START_REF] Abou-Sayed | An Experimental Technique for Measuring the Fracture Toughness of Rocks under Downhole Stress Conditions[END_REF][START_REF] Chen | Laboratory measurement and interpretation of the fracture toughness of formation rocks at great depth[END_REF] and without notches [START_REF] Stoeckhert | Mode i fracture toughness of rock under confining pressure[END_REF]. Published results on Indiana limestone are shown in Figure 3.11 and the data suggest the fracture toughness dependency on the confining pressure with a linear relationship. Provided increasing reports on confining pressure dependent fracture toughness, theoretical works to describe the mechanisms focus mainly on process zones ahead of the fracture as a culprit of the "apparent" fracture toughness including Dugdale type process zone [START_REF] Zhao | Determination of in situ fracture toughness[END_REF][START_REF] Sato | Cohesive crack analysis of toughness increase due to confining pressure[END_REF], Barenblatt cohesive zone model [START_REF] Allan M Rubin | Tensile fracture of rock at high confining pressure: implications for dike propagation[END_REF], and Dugdale-Barenblatt tension softening model [START_REF] Hashida | Numerical simulation with experimental verification of the fracture behavior in granite under confining pressures based on the tension-softening model[END_REF][START_REF] Fialko | Numerical simulation of high-pressure rock tensile fracture experiments: Evidence of an increase in fracture energy with pressure?[END_REF]. The burst experiment developed by Abou-Sayed [START_REF] Abou-Sayed | An Experimental Technique for Measuring the Fracture Toughness of Rocks under Downhole Stress Conditions[END_REF] is one of the most important methods to determine the critical stress intensity factor of rocks subject to confining pressure in the petroleum industry as the geometry closely represents actual downhole conditions of hydraulic fracturing stimulation (Figure 3.12). A hydraulic internal pressure is applied on a jacketed borehole of the thick-walled cylinder with pre-cut notches. Also, a confining pressure is applied on the outer cylinder. The inner and the outer pressures increase keeping a constant ratio of the outer to the inner pressure until the complete failure of the sample occurs and the inner and outer pressures will equilibrate to the ambient pressure abruptly. This test has great advantages in sample preparation, no fluid leak off to the rock, and easeness of measurement with various confining pressures. In this section, we firstly revisit the derivation of the stress intensity factor and analyze stabilities of fracture growth from actual burst experiment results. Subsequent analytical results indicate that fracture growth is not necessarily unstable and can have a stable phase in our experiments. In fact, stable fracture propagation has been observed also in past studies with PMMA samples [START_REF] Clifton | Determination of the critical-stress-intensity factor KIc from internally pressurized thick-walled vessels -The critical-stress-intensity factor of thick-walled rock specimens is determined from the pressure at failure[END_REF] and sandstone and shale rocks without confining pressure [START_REF] Zhixi | Determination of rock fracture toughness and its relationship with acoustic velocity[END_REF].

Evaluation and computation of the stress intensity factor for the burst experiment

Under Griffith's theory and for a given geometry (a, b, L) see Figure 3.12, the fracture stability is governed by,

K I (P i , L, b, a, P o ) ≤ K Ic
where K Ic is a material property named the critical fracture toughness. The stress intensity factor (SIF) denoted K I is such that, K I < 0 when crack lips interpenetrate and K I ≥ 0 otherwise. Let us define dimensionless parameters as,

w = b a , l = L b -a , r = P o P i (3.10)
Hence, the dimensionless crack stability becomes

K * I (1, l, w, r) ≤ K Ic (P i √ aπ) (3.11)
where K * I (1, l, w, r) = K I (1, l, w, r)/ √ aπ. Necessarily, the inner pressure must be positive

P i > 0 to propagate the crack.
For a given thick wall ratio w and pressure confinement r, we are able to evaluate the fracture toughness of the material by computing K * I if the experiment provides a value of the inner pressure P i and the crack length L at the time when the fracture propagates. The difficulty is to measure the fracture length in-situ during the experiment whose technique is yet to be established. However the burst experiment should be designed for unstable crack propagation. The idea is to maintain the crack opening by keeping the tensile load at the crack tips all along the path, so that the sample bursts (unstable crack propagation) after initiation. Therefore the fracture toughness is computed for the pre-notch length and the critical pressure measured.

3.4. Fracture stability in the burst experiment with a confining pressure Let us study the evolution of K * I (1, l, w, r) with the crack length l for the parameter analysis (w, r) to capture stability crack propagation regimes.

Using Linear Elastic Fracture Mechanics (LEFM) the burst problem denoted (B) is decomposed into the following elementary problems: a situation where pressure is applied only on the inner cylinder called the jacketed problem (J) and a problem with only a confining pressure applied on the outer cylinder problem named (C). This decomposition is illustrated in Figure 3.13. Therefore, the SIF for (B) can then be superposed as

K B * I (1, l, w, r) = K J * I (1, l, w) -rK C * I (1, l, w) (3.12)
where K C * I (1, l, w) is positive for negative applied external pressure P o . In Abou-Sayed [START_REF] Abou-Sayed | An Experimental Technique for Measuring the Fracture Toughness of Rocks under Downhole Stress Conditions[END_REF] the burst problem is decomposed following the Figure 3.14 such that, the decomposition is approximated by the jacketed problem (J) and the unjacketed problem (U) in which the fluid pressurized all internal sides. We get the following SIF,

K B * I (1, l, w, r) ≈ K J * I (1, l, w) -rK U * I (1, l, w) (3.13)
where K U * I (1, l, w) ≥ 0 for a positive P o applied in the interior of the geometry. Note that in our decomposition, no pore pressure (P p ) is considered in the sample, i.e. a drain evacuates the embedded pressure in the rock.

L 2a 2b L L 2a 2b L L 2a 2b L L 2a 2b L = + L 2a 2b L + = = L 2a 2b L L 2a 2b L L 2a 2b L - (B) (J) (C) (U) (B) (J) (B) (M)
+ = = L 2a 2b L L 2a 2b L L 2a 2b L - (B) (J) (C) (U) (B) (J) (B) (M)
Normalized stress intensity factor for the jacketed and unjacketed problems have been derived in Clifton [START_REF] Clifton | Determination of the critical-stress-intensity factor KIc from internally pressurized thick-walled vessels -The critical-stress-intensity factor of thick-walled rock specimens is determined from the pressure at failure[END_REF]. The Figure 3.15 shows a good agreement between our results (computational SIF based on the G θ methods) and one provided by Clifton [START_REF] Clifton | Determination of the critical-stress-intensity factor KIc from internally pressurized thick-walled vessels -The critical-stress-intensity factor of thick-walled rock specimens is determined from the pressure at failure[END_REF]. The G θ technique [START_REF] Destuynder | Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile[END_REF][START_REF] Suo | On the application of g (θ) method and its comparison with de lorenzi's approach[END_REF] is an estimation of the second derivatives of the potential energy with respect to the crack length, i.e. to make a virtual perturbation of the domain (vector θ) in the crack propagation direction. Then, the SIF is calculated using Irwin formula K I = EG/(1ν 2 ) based on the computed G.

Influence of the confinement and wall thickness ratio on stability of the initial crack

Based on the above result we compare K C * I with K U * I (Figures 3.13 and 3.14), and we found out their relative error is less than 15% for l ∈ [.2, .8] and w ∈ {3, 7, 10}. So, in a first approximation both problems are similar.

For the burst experiment, the fracture propagation occurs when (3.11) becomes an equality, thus we have P i = K Ic /(K B * I √ aπ). A decreasing K B * I induces a growing P i , a contrario a growing K B * I implies to decrease the inner pressure which contradicts the burst experiment set up (monotonic increasing pressure). Consequently the fracture growth is unstable (brutal) for a growing K B * I , and vice versa. In the Figure 3.16 we show different evolutions of the stress intensity factor with the crack length for various wall thickness ratio and confinement. We observe that when the confining pressure r increases fractures propagation are contained and the same effect is noticed for larger thick wall ratio w.

Depending where the pre-fracture tip is located we can draw different fracture regime summarized in three possible evolutions (see Figure 3.17 (a) For this evolution K B * I is strictly increasing thus for any pre-fracture length l 0 the sample will burst. The idea is the fracture initiates once the pressure is critical, then propagates along the sample until the failure. A sudden drop of the pressure is measured signature of the initiation pressure. By recording this pressure P i the fracture toughness K Ic is calculated using equation (3.11).

(b) By making a pre-fracture l 0 ≥ l SU , this leads to the same conclusion than (a).

However for l U S ≤ l 0 ≤ l SU the fracture propagation is stable. To get an estimation of the fracture toughness, we need to track the fracture and to measure its length otherwise is vain. A risky calculation is to assume the fracture initiation length be at the inflection point l SU before the burst. Reasons are the critical point can be a plateau shape leading to imprecise measure of l SU , secondly, since the rock is not a perfect brittle materials the l SU can be slightly different.

(c) For Griffith and any cohesive models which assume compressive forces in front of the notch tips, the fracture propagation is not possible. Of course others initiation criterion are possible as critical stress as an example.

Application to sandstone experiments

A commercial rock mechanics laboratory provided fracture toughness results for different pressure ratios on sandstones and the geometries summarized in the Table 3.3. As their end-caps and hardware are built for 0.25' center hole diameter with 2.5" diameter sample, w values are restricted to 9. Considering no pore pressure and applying stricto sensu the following equation . l U S is a critical point from unstable to stable crack propagation, vice versa for l SU . The fracture does not propagates at stop point denoted l ST by taking l equals to the dimensionless pre-notch length l 0 and the critical pressure recorded P i = P ic , we obtain that the fracture toughness K Ic is influenced by the confining pressure r as reported in the last column of the Table 3.3. However, the evolutions of K B * I with respect to l in the Figure 3.18 (right) shows that all confining experiments (Id 1-5) have a compressive area in front of the fracture tips. Moreover pre-fractures are located in the stable propagation regime, in fine the sample cannot break according to Griffith's theory.

P i √ aπK B * I (1, l, w, r) = K Ic , Chapter 3.

Sample ID 2a

[in] w The wall thickness cylinder w and the confining pressure ratio r play a fundamental role in the crack stability regime, to obtain a brutal fracture propagation after initiation smaller (w, r) is required. A possible choice is to take w = 3 for r = {1/8, 1/6} as shown in Figure 3. [START_REF] Bažant | Scaling of Structural Strength[END_REF].

P ic [Psi] r l 0 K Ic [Psi √ in] Id 0 0.
A stable-unstable regime is observed for (r = 1/6, w = 5). We performed a numerical simulation with the phase-field model to hydraulic fracturing to verify the ability of the simulation to capture the bifurcation point. For that we fix K Ic = 1, the geometric parameters (a = 1, b = 5, l 0 = .15, r = 1/5) and the internal length = 0.01. Then, by pressuring the sample (driven-pressure) damage grows until the critical point. After this 3.4. Fracture stability in the burst experiment with a confining pressure loading, the damage jumps to the external boundary and break the sample. The normalized SIF is computed using the K Ic /(P i √ aπ) for different fracture length and reported in the Figure 3.19

Remark 4 Stability analysis can be also done by volume-driven injection into the inner cylinder using phase-field models. This provides stable fracture propagation, and normalized stress intensity factor can be rebuild using simulations outputs.

Conclusion

Through this chapter we have shown that the phase-field models for hydraulic fracturing is a good candidate to simulate fractures propagation in the toughness dominated regime. The verification is done for a single-fracture and multi-fracking propagation scenario.

Simulations show that the multi-fractures propagation is the worst case energetically speaking contrary to the growth of a single fracture in the network which is the best total energy minimizer. Moreover the bifurcation to a loss of symmetries (e.g. single fracture tip propagation) is intensified by the density of fractures in the network. The pressure-driven burst experiment focuses on fracture stability. The confining pressure and the thickness of the sample might contain fractures growth. By carefully selecting those two parameters (confinement pressure ratio and the geometry) the experiment can be designed to calculate the fracture toughness for rocks.

In short those examples illustrate the potential of the variational phase-field models for hydraulic fracturing associated with the minimization principle to account for stable volume-driven fractures. The loss of symmetry in the multi-fracking scenario is a relevant example to illustrate the concept of variational argument. Same results is confirmed by coupling this model with fluid flow as detailed in Chukwudozie [START_REF] Chukwudozie | Application of the Variational Fracture Model to Hydraulic Fracturing in Poroelastic Media[END_REF].

Substituting (3.19) into (3.14), the fluid pressure is obtained.

p = 3 2 G 2 c E π V (3.20)
Similarly, the fracture length during propagation is obtained by substituting (3.16) into (3.14).

l = 3 E V 2 4π G c (3.21)
Penny-Shaped Fracture (3d domain):

For a penny-shaped fracture in a 3d domain, the fracture volume is

V = 16pl 3 3E (3.22)
where l denotes the radius, while the critical fluid pressure [START_REF] Sneddon | The opening of a griffith crack under internal pressure[END_REF] is

p c = πG c E 4l 0 (3.23)
For an initial fracture radius l 0 , the critical volume is,

V c = 64πl 5 0 G c 9E (3.24)
If one follows a procedure similar to that for the line fracture, we will obtain the following relationships for the evolution of the fluid pressure and fracture radius

p c = 5 π 3 G 3 c E 2 12V l = 5 9E V 2 64πG c (3.25)
Chapter 4. Variational models of perfect plasticity Definition 7 (Generalized standard plasticity models)

i. A choice of independent states variables which includes one or multiple internal variables.

ii. Define a convex set where thermodynamical forces lie in.

Concerning i. we choose the plastic strain tensor (symmetric) p and the infinitesimal total deformation denoted e(u). The total strain is the symmetrical part of the spatial gradient of the displacement u, i.e.

e(u) = ∇u + ∇ T u 2 .
The kinematic admissibility is the sum of the plastic and elastic strains denoted ε, given by, e(u) = ε + p.

For ii. consider a free energy density ψ a differentiable convex state function which depends on internal variables. Naturally, thermodynamical forces are defined from the free energy by

σ = ∂ψ ∂e (e, p), τ = - ∂ψ ∂p (e, p). (4.1)
Commonly the free energy takes the form of ψ(e, p) = 1 2 A(e(u)p) : (e(u)p), where A is the Hooke's law tensor. It follows that, σ = τ = A(e(u)p). However for clarity we continue to use τ instead. Internal variables (e, p) and their duals (σ, τ ) are second order symmetric tensors and become n × n symmetric matrices denoted M n s after a choice of an orthonormal basis and space dimension of the domain (Ω ⊂ R n ). To complete the second statement ii., let K be a non empty closed convex subset of M n s where τ lies in. This subset is called elastic domain for τ . Assume that K is fixed and time independent, such that its boundary is the convex yield surface f Y : M n s → R, defined by,

τ ∈ K = {τ * ∈ M n s : f Y (τ * ) ≤ 0} (4.2)
Precisely, for any τ that lies in the interior of K denoted int(K) the yield surface is strictly negative. Otherwise, τ belongs to the boundary noted ∂K and the yield function vanishes:

f Y (τ ) < 0, τ ∈ int(K) f Y (τ ) = 0, τ ∈ ∂K . (4.3)
Let us apply the normality rule on it to get the plastic evolution law. In the case where ∂K is differentiable the plastic flow rule is defined as, 4.1. Ingredients for generalized standard plasticity models

ṗ = η ∂f Y ∂τ (τ ), with η = 0 if f Y (τ ) < 0 ≥ 0 if f Y (τ ) = 0 (4.4)
where η is the Lagrange multiplier. Sometimes the convex K has corners and the outer normal cannot be defined (f Y is not differentiable), thus, the normality rule is written using Hill's principle, also known as maximum dissipation power principle, i.e.,

τ ∈ K, (τ -τ * ) : ṗ ≥ 0, ∀τ * ∈ K. (4.5)
This is equivalent to say that ṗ lies in the outer normal cone of K in τ ,

ṗ ∈ N K (τ ) := { ṗ : (τ * -τ ) ≤ 0 ∀τ * ∈ K}. (4.6)
However we prefer to introduce the indicator function of τ ∈ K, and write equivalently the normality rule as, ṗ lies in the subdifferential set of the indicator function. For that, the indicator function is,

I K (τ ) = 0 if τ ∈ K +∞ if τ / ∈ K (4.7)
and is convex by construction. The normality rule is recovered by applying the definition of subgradient, such that, ṗ is a subgradient of I K at a point τ ∈ K for any τ * ∈ K, given by,

τ ∈ K, I K (τ * ) ≥ I K (τ ) + ṗ : (τ * -τ ), ∀τ * ∈ K ⇔ ṗ ∈ ∂I K (τ ), τ ∈ K (4.8)
where the set of all sub-gradients at τ is the sub-differential of I K at τ and is denoted by ∂I K (τ ). At this stage of the analysis, Hill's principle is equivalent to convex properties of the elastic domain K and the normality plastic strain flow rule.

For τ ∈ K, Hill ⇔ ṗ ∈ N K (τ ) ⇔ ṗ ∈ ∂I K (τ ) (4.9)

Dissipation of energy during plastic deformations

All ingredients are settled, such as, we have the variable set (u, p) and their duals (σ, τ ) which lie in the convex set K. Also, the plastic evolution law is given by ṗ ∈ ∂I K (τ ).

It is convenient to compute the plastic dissipated energy during a plastic deformation process. For that, the dissipated plastic power density can be constructed from the Clausius-Duhem inequality. To construct such dissipation energy let us define first the support function

H(q), q ∈ M 3 s → H(q) := sup τ ∈K {τ • q} ∈ (-∞, +∞] (4.10)
The support function is convex, 1-homogeneous, Chapter 4. Variational models of perfect plasticity

H(λq) = λH(q), ∀λ > 0, ∀q ∈ M n s (4.11)
and it follows the triangle inequality, i.e., H(q 1 + q 2 ) ≤ H(q 1 ) + H(q 2 ), for every q 1 , q 2 ∈ M n s .

(4.12)

The support function of the plastic strain rate H( ṗ) is null if the plastic flow is zero, non negative when 0 ∈ K, and takes the value +∞ when K is not bounded. Using Clausius-Duhem inequality for an isotherm transformation, the dissipation power is defined by

D = σ : ė -ψ, (4.13)
and the second law of thermodynamics enforce the dissipation to be positive or null,

D = τ : ṗ ≥ 0. (4.14)
Using Hill's principle, the definition of the support function and some convex analysis, one can show that the plastic dissipation is equal to the support function of the plastic flow.

D = H( ṗ) (4.15)
The starting point to prove (4.15) is the Hill's principle or equivalently the plastic strain flow rule.

For τ ∈ K, τ : ṗ ≥ τ * : ṗ, ∀τ * ∈ K. By passing the right term to the left and taking the supremum over all ṗ ∈ M n s , we get,

sup ṗ∈M n s {τ : ṗ -H( ṗ)} ≥ 0. (4.18)
Since K is a non empty close convex set, H( ṗ) is convex and lower semi continuous, we have built the convex conjugate function of H( q) in the sense of Legendre-Fenchel. Moreover, one observes that the conjugate of the support function is the indicator function, given by,

I K (τ ) := sup ṗ∈M n s {τ : ṗ -H( ṗ)} = 0 if τ ∈ K + ∞ if τ / ∈ K (4.19)
Hence, the following equality holds for τ ∈ K,

D = τ : ṗ = H( ṗ). (4.20)

Variational formulation of perfect plasticity models

Remark 5 The conjugate subgradient theorem says that, for τ ∈ K a non empty closed convex set,

ṗ ∈ ∂I K (τ ) ⇔ D = τ : ṗ = H( ṗ) + I K (τ ) ⇔ τ ∈ ∂H( ṗ)
Finally, once the plastic dissipation power defined, by integrating over time [t a , t b ] for smooth evolution of p, the plastic dissipated energy is,

D(p; [t a , t b ]) = t b ta H( ṗ(s)) ds (4.21)
This problem is rate independent because the dissipation does not depend on the strain rate , i.e. D( ė, ṗ) = D( ṗ) and is 1-homogeneous.

Variational formulation of perfect plasticity models

Consider a perfect elasto-plastic material with a free energy ψ(e, p) occupying a smooth region Ω ⊂ R n , subject to time dependent boundary displacement ū(t) on a Dirichlet part ∂ D Ω of its boundary. For the sake of simplicity the domain is free of stress and no body force applies on it, such that, σ • ν = 0 on the complementary portion ∂ N Ω = ∂Ω \ ∂ D Ω, where ν denotes the appropriate normal vector. Assume the initial state of the material being (e 0 , p 0 ) = (0, 0) at t = 0. Internal variables e(u) and p are supposed to be continuous-time solution of the quasi-static evolution problem. At each time the body is in elastic equilibrium with the prescribed loads at that time, such as it satisfies the following equations,

                 σ = ∂ψ ∂e (e, p) in Ω τ = - ∂ψ ∂p (e, p) ∈ ∂H( ṗ) in Ω div(σ) = 0 in Ω u = ū(t) on ∂ D Ω σ • ν = 0 on ∂ N Ω
We set aside problems where plasticity strain may develop at the interface ∂ D Ω. The problem can be equivalently written in a variational formulation, which is based on two principles, i. Energy balance ii. Stability condition Let the total energy density be defined as the sum of the elastic energy and the dissipated plastic energy, E t (e(u), p) = ψ(e(u), p)ψ(e 0 , p 0 ) + D(p; [0, t])

Energy balance

The concept of energy balance is related to the evolution of state variables in a material point, and enforce the total energy rate be equal to the mechanical power energy at each time, i.e. Ėt = σ t : ėt .

(

The total energy rate is,

Ėt = ∂ψ ∂e (e t , p t ) : ėt + ∂ψ ∂p (e t , p t ) : ṗt + H( ṗt ), (4.23)
and using the definition of τ = -∂ψ/∂e and σ = ∂ψ/∂e, we obtain,

τ t • ṗt = sup τ ∈K {τ : ṗt } (4.24)

Stability condition for the plastic strain

The stability condition for p is finding stable p t ∈ M n s for a given loading deformation e t . We propose to approximate the continuous time evolution by a time discretization, such that,

0 = t 0 < • • • < t i < • • • < t N =
t b and at the limit max i |t it i-1 | → 0. At the current time t i = t, let the material be at the state e t i = e and p t i = p and the previous state (e t i-1 , p t i-1 ). The discretized plastic strain rate is ṗt (p-p t i-1 )/(t-t i-1 ). During the laps time from t i-1 to t the increment of plastic energy dissipated is t t i-1 H( ṗt )ds H(pp t i-1 ). Hence taking into account all small previous plastic dissipated energy events, the total dissipation is approximated by,

D(p) := H(p -p t i-1 ) + D(p t i-1 ) (4.25)
At the current time, a plastic strain perturbation is performed for a fixed total strain changing the system from (e, p) to (e, q). The definition of the stability condition adopted here is written as a variation of the total energy between this two states, p stable, e given ⇔ ψ(e, q) + H(qp t i-1 ) ≥ H(pp t i-1 ) + ψ(e, p), ∀q ∈ M 3 s (4.26)

We wish to highlight the stability definition adopted, which is for infinitesimal transformations the flow rule.

H(qp t i-1 ) ≥ H(pp t i-1 ) -ψ(e, q)ψ(e, p) qp : (qp), ∀q ∈ M n s , q = p (4.27)

Consider small variations of the plastic strain p in the direction p for a growing total energy, such that for some h > 0 small enough and p + hp ∈ M n s we have, Using the Legendre transform, we get,

τ ∈ ∂H(p -p t i-1 ) ⇔ (p -p t i-1 ) ∈ ∂I K (τ ). (4.29)
To recover the continuous-time evolution stability for p, divide by δt = tt i-1 and pass δt to the limit. We recover the flow rule ṗ ∈ ∂I K (τ ), or equivalently in the conjugate space τ ∈ ∂H( ṗ).

Let us justify the definition adopted of the stability by showing that there is no lowest energy that can be found for a given e t . Without loss of any generality assume a continuous straight smooth path p(t) starting at p(0) = p and finishing at p(1) = q, such as, (4.31)

t ∈ [0, 1] → p(t) = (1 -t)p + tq, ∀q ∈ M n s (4.
The right hand side is path independent, by taking the infimum over all plastic strain paths, we get,

inf t →p(t) p(0)=p, p(1)=q 1 0 H( ṗ(s))ds ≥ τ * : (q -p) (4.32)
The left hand side does not depends on τ * , taking the supremum for all τ * ∈ K, and applying the triangle inequality for any p t i-1 , one obtains,

inf t →p(t) p(0)=p, p(1)=q 1 0 H( ṗ(s))ds ≥ H(q -p) ≥ H(q -p t i-1 ) -H(p -p t i-1 ). (4.33)
which justifies the a posteriori adopted definition of the stability.

The stability condition for the displacement is performed on the first chapter and we simply recover the equilibrium constitutive equations for the elastic problem with the prescribed boundary conditions.

Numerical implementation and verification of perfect

elasto-plasticity models

- ∂ N Ω g(t) • u dH n-1 ,
where H n-1 denotes the Hausdorff n -1-dimensional measure of the boundary. Typical plastic yields criterion used for metal are Von Mises or Tresca, which are well known to have only a bounded deviatoric part of the stress, thus they are insensitive to any stress hydrostatic contributions. Consequently, the plastic strain rate is also deviatoric ṗ ∈ dev(M n s ) and it is not restrictive to assume that p ∈ dev(M n s ). For being more precise but without going into details, existence and uniqueness is given for solving the problem in the stress field, σ ∈ L 2 (Ω; M n s ) (or e(u) ∈ L 2 (Ω; M n s ) ) with a yield surface constraint σ ∈ L ∞ (Ω; dev(M n s )). Experimentally it is observed that plastic strain deformations concentrate into shear bands, as a macroscopic point of view this localization creates sharp surface discontinuities of the displacement field. In general the displacement field cannot be solved in the Sobolev space, but find a natural representation in a bounded deformation space u ∈ BD(Ω) when the plastic strain becomes a Radon measure p ∈ M(Ω ∪ ∂ D Ω; dev(M 3 s )). The problem of finding (u, p) minimizing the total energy and satisfying the boundary conditions is solved by finding stable states variables trajectory i.e. stationary points. This quasi-static evolution problem, is numerically approximated by solving the incremental time problem, i.e. for a given time interval [0, T ] subdivided into (N + 1) steps we have,

0 = t 0 < t 1 < • • • < t i-1 < t i < • • • < t N = T .
The discrete problem converges to the continuous time evolution provided max i (t it i-1 ) → 0, and the total energy at the time t i in the discrete setting is,

E t i (u i , p i ) = Ω 1 2 A(e(u i ) -p i ) : (e(u i ) -p i ) + D i (p i ) dx - ∂ N Ω g(t i ) • u i dH n-1
where, 4.3. Numerical implementation and verification of perfect elasto-plasticity models

D i (p i ) = H(p i -p i-1 ) + D i-1 (4.34)
for a prescribed u i = ūi on ∂ D Ω. Let i be the the current time step, the problem is finding (u i , p i ) that minimizes the discrete total energy, i.e (u i , p i ) := argmin

u∈C i p∈M(Ω∪∂ D Ω;dev(M 3 s )) E t i (u, p) (4.35)
where p = (ū iu) • ν on ∂ D Ω and C i is the set of admissible displacement,

C i = {u ∈ H 1 (Ω) : u = ūi on ∂ D Ω}.
The total energy E(u, p) is quadratic and strictly convex in u and p separately. For a fixed u or p, the minimizer of E(•, p) or E(u, •) exists, is unique and can easily be computed. Thus, a natural algorithm technique employed is the alternate minimization detailed in Algorithm 3, where δ p is a fixed tolerance.

More precisely, at the loading time t i , for a given p j i , let find u j i that minimizes E(u, p j i ), notice that the plastic dissipation energy does not depend on the strain e(u), thus,

u j i := argmin u∈C i Ω 1 2
A(e(u)p j i ) : (e(u)p j i ) dx -

∂ N Ω g(t) • u dH n-1 (4.36)
This is a linear elastic problem. Then, for a given u j i let find p on each element cell, such as it minimizes E(u j i , p). This problem is not easy to solve in the primal formulation,

p j i := argmin p∈M(Ω∪∂ D Ω;dev(M n s )) 1 2
A(e(u j i )p) : (e(u j i )p) + H(pp i-1 )

but from the previous analysis, the stability condition of this problem is

A(e(u j i )p) ∂H(pp i-1 ). Using the Legendre-transform, the stability of the conjugate problem is given by

(p -p i-1 ) ∈ ∂I K (A(e(u j i ) -p)).
One can recognize the flow rule in the discretized time. This is the stability condition of the problem,

p j i := argmin p∈M(Ω∪∂ D Ω;dev(M n s )) A(e(u j i )-p)∈K 1 2 A(p -p i-1 ) : (p -p i-1
).

The minimization with respect to u is a simple linear problem solved using preconditioned conjugated gradient while minimization with respect to p can be reformulated Solve the equilibrium,

u j+1 := argmin u∈C i E i (u, p j ) 4:
Solve the plastic strain projection on each cell, p j+1 := argmin j := j + 1 6: until p jp j-1 L ∞ ≤ δ p 7: Set, u i := u j and p i := p j

Numerical verifications

A way to do a numerical verification is to recover the closed form solution of a bi-axial test in 3D provided in [START_REF] Gilles A Francfort | A case study for uniqueness of elasto-plastic evolutions: The bi-axial test[END_REF].

In the fixed orthonormal basis (e 1 , e 2 , e 3 ), consider a domain Ω = (-d/2, d/2) × (-l/2, l/2) × (0, l), (d < l), with the boundary conditions:

   σ 11 = 0 on x 1 = ±d/2 σ 22 = g 2
on x 1 = ±l/2 σ 13 = σ 23 = 0 on x 3 = 0, l and, add

u 3 = 0 on x 3 = 0 u 3 = tl on x 3 = l.
Considering the classical problem to solve,

   div(σ) = 0 in Ω σ = Ae(u)
in Ω e(u) = (∇u + ∇ T u)/2 in Ω constrained by a Von Mises plasticity yield criterion,

3 2 dev(σ) : dev(σ) ≤ σ p
It is shown in [START_REF] Gilles A Francfort | A case study for uniqueness of elasto-plastic evolutions: The bi-axial test[END_REF] that the domain remains elastic until the plasticity is triggered at a critical loading time t c as long as

0 ≤ g 2 ≤ σ p / √ 1 -ν + ν 2 , t c = 1 2E (1 -2ν)g 2 + 4σ 2 p -3g 2 2
where (E, ν) denote respectively the Young's modulus and the Poisson ratio. For 0 ≤ t ≤ t c the elastic solution stands for

       σ(t) = g 2 e 2 ⊗
                             σ(t) =g 2 e 2 ⊗ e 2 + σ3 e 3 ⊗ e 3 , σ3 = 1 2 g 2 + 4σ 2 p -3g 2 2 e(t) = -ν(1 + ν) g 2 E e 1 ⊗ e 1 + (1 -ν 2 ) g 2 E e 2 ⊗ e 2 + t(-νe 1 ⊗ e 1 -νe 2 ⊗ e 2 + e 3 ⊗ e 3 ) p(t) =(t -t c ) - g 2 + σ3 2 σ3 -g 2 e 1 ⊗ e 1 + 2g 2 -σ3 2 σ3 -g 2 e 2 ⊗ e 2 + e 3 ⊗ e 3 u(t) = -ν(1 + ν) g 2 E -νt c - g 2 + σ3 2 σ3 -g 2 (t -t c ) x 1 e 1 + (1 -ν 2 ) g 2 E -νt c + 2g 2 -σ3 2 σ3 -g 2 (t -t c )
x 2 e 2 + tx 3 e 3 (4.38) A numerical simulation has been performed on a domain parametrized by l = .5 and d = .2, pre-stressed on opposite faces by g 2 = .5 with the material parameters E = 1,σ p = 1 and a Poisson ratio set to ν = .3. For those parameters, numerical results and exact solution have been plotted see Figure 4.1, and matches perfectly.

One difficulty is to get closed form for different geometry and plasticity criterion. Alternate minimization technique converge to the exact solution on this example for Von Mises in 3D

Conclusion

The adopted strategy to model a perfect elasto-plastic material is to prescribe the elastic stress domain set (closed convex) with plastic yields functions without dealing with corners and approximate the continuous evolution problem by discretized time steps. The 99 Chapter 4. Variational models of perfect plasticity implemented algorithm solves alternately the elastic problem and the plastic projection onto the yield surface. Hence, there is no difficulty to implement other perfect plastic yield criteria. A verification is performed on the biaxial test for Von Mises plastic yield criteria.

Chapter 5

Variational phase-field models of ductile fracture by coupling plasticity with damage Phase-field models referred to as gradient damage models of brittle fracture are very efficient to predict cracks initiation and propagation in brittle and quasi-brittle materials [START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | Numerical implementation of a variational formulation of quasi-static brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]. They were originally conceived as an approximation of Francfort Marigo's variational formulation [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] which is based on Griffith's idea of competition between elastic and fracture energy. Their model inherits a fundamental limitation of Griffith's theory which is a discontinuity of the displacement belongs to the damage localization strip, and this is not observed during fractures nucleation in ductile materials. Moreover, they cannot be used to predict cohesive-ductile fractures since no permanent deformations are accounted for. Plasticity models [START_REF] Suquet | Sur les équations de la plasticité: existence et régularité des solutions[END_REF][START_REF] Salençon | Elasto-plasticité[END_REF][START_REF] Halphen | Sur les matériaux standard généralisés[END_REF][START_REF] Maso | Quasistatic crack growth in elasto-plastic materials: The two-dimensional case[END_REF][START_REF] Babadjian | Quasi-static evolution in nonassociative plasticity: the cap model[END_REF] are widely used to handle with the aforementioned effects by the introduction of the plastic strain variable. To capture ductile fracture patterns the idea is to couple the plastic strain coming from plasticity models with the damage in the phase-field approaches to fracture.

The goal of this chapter is to extend the Alessi-Marigo-Vidoli work [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF] by considering any associated perfect plasticity and to provide a general algorithm to solve the problem for any dimensions. We provide a qualitative comparison of crack nucleation in various specimen with published experimental results on metals material. We show capabilities of the model to recover cracks patterns characteristics of brittle and ductile fractures. After the set of parameters being adjusted to recover ductile fracture we focus solely on such regime to study cracks nucleation and propagation phenomenology in mild notched specimens.

The chapter is organized as follow: Section 5.1.1 starts by aggregating some experiments illustrating mechanisms of ductile fracture which will constitute basis of numerical comparisons provided in the last part of this chapter. Section 5.1.2 is devoted to the introduction of variational phase-field models coupled with perfect plasticity and to recall some of their properties. Section 5.1.3 focuses on one dimension bar in traction to provide the cohesive response of the material and draw some fundamental properties similarly 5.1. Phase-field models to fractures from brittle to ductile to [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF]. A numerical implementation technique to solve such coupled models is provided in section 5.2. For the remainder we investigate ductile fracture phenomenology by performing simulations on various geometries such as, rectangular specimen, a mild notch 2d plane strain and 3d round bar respectively exposed in sections 5. Numerous experimental evidences show a common phenomenology of fracture nucleation in a ductile materials. To illustrate this, we have selected relevant experiments showing fracture nucleation and propagation in a plate and in a round bar.

For instance in [START_REF] Spencer | The influence of iron content on the plane strain fracture behaviour of aa 5754 al-mg sheet alloys[END_REF] the role of ductility with the influence of the iron content in the formation of shear band have been investigated. Experiments on Aluminum alloy AA 5754 Al -Mg show fractures nucleation and evolution in the thickness direction of the plate specimen illustrated in Figure 5.2.

The tensile round bar is another widely used test to investigate ductile fractures. However, tracking fractures nucleation inside the material is a challenging task and requires special equipment like tomography imaging to probe. Nevertheless Benzerga [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF][START_REF] Amine Benzerga | Synergistic effects of plastic anisotropy and void coalescence on fracture mode in plane strain[END_REF] and Luu [START_REF] Luu | Déchirure ductile des aciers à haute résistance pour gazoducs (X100)[END_REF] results show pictures of cracks nucleation and propagation inside those types of samples see Figure 5.19. A simpler method is the fractography which consists in studding fracture surfaces of materials after failure of the samples. Typical ductile fractures 5.1. Phase-field models to fractures from brittle to ductile and powerful approach to study theoretically and solve numerically those problems. The coupling between both models is done at the proposed total energy level. We start by recalling some important properties of variational phase-field models interpreted as gradient-damage models and variational perfect plasticity.

Consider an elasto-plastic-damageable material with A the Hooke's law tensor occupying a region Ω ⊂ R n in the reference configuration. The region Ω is subject to a time dependent boundary displacement ū(t) on a Dirichlet part of its boundary ∂ D Ω and time stress dependent g(t) = σ • ν on the remainder ∂ N Ω = ∂Ω \ ∂ D Ω, where ν denotes the appropriate normal vector. A safe load condition is required for g(t) to set aside issues in plasticity theory. For the sake of simplicity body forces are neglected such that at the equilibrium, the stress satisfies, div(σ) = 0 in Ω

The infinitesimal total deformation e(u) is the symmetrical part of the spatial gradient of the displacement field u, i.e. e(u) = ∇u + ∇ T u 2

Since the material has permanent deformations, it is usual in small deformations plasticity to consider the plastic strain tensor p (symmetric) such that the kinematic admissibility is an additive decomposition, e(u) = ε + p where ε is the elastic strain tensor.

The material depends on the damage variable denoted α which is bounded between two extreme states, α = 0 is the undamaged state material and α = 1 refers to the broken part. Let the damage deteriorate the material properties by making an isotropic modulation of the Hooke's law tensor a(α)A, where the stiffness function a(α) is continuous and decreasing such that a(0) = 1, a(1) = 0. In linearized elasticity the recoverable energy density of the material stands for, ψ(e(u), α, p) := 1 2 a(α)A(e(u)p) : (e(u)p)

Consequently the relation which relates the stress tensor σ to the strain is, One can recognize the Hill's principle by applying the definition of subdifferential and the indicator function. Since b(α)K is a none empty closed convex set, using Legendre-Fenchel, the conjugate of the plastic flow is σ ∈ b(α)∂H( ṗ), where the plastic dissipation potential H(q) = sup τ ∈K {τ : q} is convex, subadditive, positively 1-homogeneous for all q ∈ M n×n s . The dissipated plastic energy is obtained by integrating the plastic dissipation power over time, such that,

φ p := t 0 b(α)H( ṗ(s)) ds (5.1)
This dissipation is not unique and we have to take into account the surface energy produced by the fracture. Inspired by the phase-field models to brittle fracture [START_REF] Ambrosio | On the approximation of free discontinuity problems[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Pham | Stability of homogeneous states with gradient damage models: Size effects and shape effects in the three-dimensional setting[END_REF] we define the surface dissipation term as,

φ d := t 0 σ 2 c 2Ek w (α) α + 2 ∇α • ∇ α + b (α) α t 0 H( ṗ(s)) ds dt (5.2)
where the first term is the classical approximated surface energy in brittle fracture and the last term is artificially introduced to be combined with φ p . Precisely, after summation of the free energy ψ(e(u), α, p), the work force, the dissipated plastic energy φ p and the dissipated damage energy φ d , the total energy has the following form,

E t (u, α, p, p) = Ω 1 2 a(α)A(e(u) -p) : (e(u) -p) dx - ∂ N Ω g(t) • u dH n-1 + Ω b(α) t 0 H( ṗ(s)) ds dx + σ 2 c 2Ek Ω w(α) + 2 |∇α| 2 dx (5.3)
where p = t 0 ṗ(s) ds is the cumulated plastic strain which is embedded in the cumulated plastic dissipation energy t 0 H( ṗ(s)) ds. The surface dissipation potential w(α) is a continuous increasing function such that w(0) = 0 and up to a rescaling, w(1) = 1. Since the damage is a dimensionless variable, the introduction of ∇α enforce to have > 0 a regularized parameter which has a dimension of the length. Note that the total energy (5.3) is composed of two dissipations potentials ϕ p and ϕ d coupled where,

ϕ p = Ω b(α) t 0 H( ṗ(s)) ds dx, ϕ d = σ 2 c 2Ek Ω w(α) + 2 |∇α| 2 dx.
(5.4)

5.1. Phase-field models to fractures from brittle to ductile Taking p = 0 in (5.3), the admissible stress space is bounded by,

A -1 σ : σ ≤ σ 2 c Ek max α w (α) c (α)
where E is the Young's modulus, the compliance function is c(α) = 1/a(α) and let k = max α w (α) c (α) . Therefore, without plasticity in one dimensional setting an upper bound of the stress is σ c .

A first conclusion is the total energy (5.3) is composed of two coupled dissipation potentials associated with two yields surfaces and their evolutions will be discussed later.

In the context of smooth triplet state variable ζ = (u, α, p) and since the above total energy (5.3) must be finite, we have α ∈ H 1 (Ω) and e(u), p belong to L 2 (Ω). However, experimentally it is observed that plastic strain concentrates into shear bands. In our model since ṗ ∈ b(α)K, the plastic strain concentration is driven by the damage localization and both variables intensifies on the same confined region denoted J(ζ), where J is a set of "singular part" which a priori depends on all internal variables. Also, the damage is continuous across the normal surfaces of J(ζ) but not the gradient damage term which may jump. Accordingly, the displacement field cannot be solved in the Sobolev space, but find a natural representation in special bounded deformation space SBD if the Cantor part of e(u) vanishes, so that the strain measure can be written as,

e(u) = e(u) + u ν H n-1 on J(ζ(x))
where e(u) is the Lebesgue continuous part and denotes the symmetrized tensor product. For the sake of simplicity, consider the jumps set of the displacement being a smooth enough surface, i.e the normal ν is well defined, and there is no intersections with boundaries such that J(ζ)∩∂Ω = ∅. The plastic strain turns into a Dirac measure on the surface J(ζ). Without going into details, the plastic strain lies in a non-conventional topological space for measures called Radon space denoted M.

Until now, the damage evolution have not been set up and the plastic flow rule is hidden in the total energy adopted. Let us highlight this by considering the total energy be governed by three principles; damage irreversibility, the stability of E t (u, α, p) with respect to all admissible variables (u, α, p) and the energy balance.

We focus on the time-discrete evolution, by considering a time interval [0, T ] subdivided into (N + 1) steps such that, 0 = t

0 < t 1 < • • • < t i-1 < t i < • • • < t N = T .
The following discrete problem converges to the continuous time evolution provided max(t it i-1 ) → 0. At any time t i , the sets of admissible displacement, damage and plastic strain fields respectively denoted C i , D i and Q i are: 107 Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage

C i = u ∈ SBD(Ω) : u = ū(t i ) on ∂ D Ω , D i = α ∈ H 1 (Ω) : α i-1 ≤ α < 1 in Ω , Q i = p ∈ M( Ω; M n×n s ) such that, p = u ν on J ζ(x) (5.5)
and because plastic strains may develop at the boundary, we know from prior works on plasticity [START_REF] Dal Maso | Quasistatic evolution problems for linearly elastic-perfectly plastic materials[END_REF] that we cannot expect the boundary condition to be satisfied, thus we will have to set up p = (ū(t i )u) ν on ∂ D Ω.

It is convenient to introduce Ω ⊃ Ω a larger computational domain which includes the jump set and ∂ D Ω, this will become clearer. Note that the damage irreversibility is in the damage set D i . The total energy of the time-discrete problem is composed of (5.3) on the regular part and b(α)D ( u ν, [0, t i ]) on the singular part, such that,

E t i (u, α, p) = Ω\J(ζ) 1 2 a(α)A(e(u) -p) : (e(u) -p) dx - ∂ N Ω g(t i ) • u dH n-1 + Ω b(α)D i (p) dx + σ 2 c 2Ek Ω\J(ζ) w(α) + 2 |∇α| 2 dx (5.6)
where

D i (p) = H(p -p i-1 ) + D i-1 (5.7)
the total energy is defined over the regular and singular part of the domain, and the evolution is governed by, Definition 8 (Time discrete coupled plasticity-damage evolution by local minimization)

At every time t i find stable variables trajectory (u i , α i , p i ) ∈ C i × D i × Q i that satisfies the variational evolution: i. Initial conditions: u 0 = 0, α 0 = 0 and p 0 = 0 ii. Find the triplet ζ i = (u i , α i , p i ) which minimizes the total energy,

E t i (u, α, p)
iii. Energy balance,

E t i (u i , α i , p i ) =E t 0 (u 0 , α 0 , p 0 ) + i k=1 ∂ D Ω (σ k ν) • (ū k -ūk-1 ) dH n-1 - ∂ N Ω (g(t k ) -g(t k-1 )) • u k dH n-1
(5.8)

The damage and plasticity criterion are obtained by writing the necessary first order optimality condition of the minimizing problem E t i (u, α, p). Explicitly, there exists h > 0 small enough, such that for

(u i + hv, α i + hβ, p i + hq) ∈ C i × D i × Q i , E t i (u i + hv, α i + hβ, p i + hq) ≥ E t i (u i , α i , p i ) (5.9)
Consider that the displacement at u i in the direction v might extend the jump set of J(v). The variation of the total energy

E t i (u i + hv, α i + hβ, p i + hq) is equal to, Ω\(J(ζ i )∪J(v)) 1 2 a(α i + hβ)A e(u i + hv) -(p i + hq) : e(u i + hv) -(p i + hq) dx - ∂ N Ω g(t i ) • (u i + hv) dH n-1 + Ω\(J(ζ i )∪J(v)) b(α i + hβ)D i (p i + hq) dx + J(ζ i )∪J(v) b(α i + hβ)D i (( u i + hv) ν) dH n-1 + σ 2 c 2Ek Ω\(J(ζ i )∪J(v)) w(α i + hβ) + 2 |∇(α i + hβ)| 2 dx
(5.10) Note that the plastic dissipation term is split over the regular part and the singular part and for simplicity we set aside the plastic strain localization on the Dirichlet boundary.

Equilibrium and kinematic admissibility:

Take β = 0 and q = 0 in (5.9) and (5.10) such that E t i (u i +hv, α i , p i ) ≥ E t i (u i , α i , p i ). Using (5.7) we just have to deal with the current plastic potential H which is subadditive and 1-homogeneous. Hence, the fourth term in (5.10) becomes, 5.1. Phase-field models to fractures from brittle to ductile 2. The damage yield criteria in the bulk,

f D (σ t , α t (x), pt (x)) := - 1 2 c (α t (x)) E σ 2 t + σ 2 c 2kE w (α t (x)) -2 2 α t (x) + b (α t (x))σ p pt (x) ≥ 0 (5.31) 3. The damage yield criteria on x 0 , b (α t (x 0 )) u(x 0 ) σ p - 2 σ 2 c kE α t (x 0 ) ≥ 0 (5.32)
4. The damage yield criteria on ±L,

α t (-L) ≥ 0, α t (L) ≤ 0 (5.33)
5. The plastic yield criteria in the bulk and on the jump, We restrict our study to r Y = σ c /σ p > 1, meaning that the plastic yield surface is below the damage one. Consequently after the elastic stage, the bar will behave plastically. During the plastic stage, the cumulation of plastic strain decreases f D until the damage yield criteria is reached. On the third stage both damage and plasticity evolves simultaneously such that f D = 0 and f Y = 0 on the jumps x 0 . Of course there is no displacement jump on the bar before the third stage. Let expose the solution (u, α, p) Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage for the elastic, plastic and plastic damage stages.

f Y (σ t , α t (x)) := |σ t | -b(α t (x))σ p ≤ 0 ( 5
The elastic response of the bar ends once the tension reached u t = σ p /E. During this regime the damage and plastic strain remain equal to zero. After this loading point, the plasticity stage begins and we have a uniform p = p = u tσ p /E and α = 0 in Ω. Since b (α) < 0 and p increases during the plastic stage, the damage yield criteria f D decreases until the inequality (5.31) becomes an equality. At this loading time both criterion are satisfied, such that, f Y = 0 and f D = 0. Hence, plugging the equation (5.34) into (5.31), we get,

-b (α t (x))p t (x) = σ p E - 1 2 c (α t (x))b 2 (α t (x)) + r 2 Y 2k w (α t (x)) -2 2 α t (x) (5.39)
By taking α t (x) = 0 in the above equation, we get the condition when the plastic stage ends, for a uniform plastic strain,

p = u t - σ p E = σ p (-b (0))E r 2 Y 2k w (0) - 1 2 c (0)b 2 (0) (5.40)
The last stage is characterized by the evolution of the damage. For a given x 0 take L long enough to avoid any damage perturbation at the boundary such that, the damage remains equal to zero at the extremities of the bar α(±L) = 0 and assume being maximum at x 0 , α(x 0 ) = β. Let α ≥ 0 over [-L, x 0 ) with α (-L) = 0, multiplying the equation (5.31) by 2α and integrate over [-L, x 0 ), we get,

- 2E σ p x 0 -L b (α t (x))α t (x)p t (x) dx = c(β) -c(0) σ 2 t σ 2 p + r 2 Y k w(β) -2 β 2 (5.41)
A priori, the cumulated plastic strain evolves along the part of the bar [-L, x 0 ), but since the maximum damage value β is reached on x 0 and the stress is uniform in the bar we have σ t (x) ≤ b(β)σ p . In other words the plasticity does not evolve anymore in the bar except on x 0 , and p is equal to (5.40). We obtain a first integral of the form of,

2 β 2 = k r 2 Y c(β) -c(0) b 2 (β) + w(β) + 2 b(β) -b(0) p Ek σ p r 2 Y (5.42)
We know that on the jump set, we have,

b (β) u(x 0 ) σ p - 2 σ 2 c kE β = 0 (5.43)
Since β is known, the stress on the bar and the displacement jump on x 0 can be computed. We define the energy release rate as the energy dissipated by the damage process, 5.2. Numerical implementation of the gradient damage models coupled with perfect plasticity

G t := Ω\{x 0 } σ c 2kE w(α t (x)) + 2 α 2 t (x) + b(α t (x))σ p p dx + b(α t (x 0 ))σ p u(x 0 )
and the critical value is given for complete damage localization once σ = 0.

Let us recall some fundamental properties for a(α), b(α) and w(α) to satisfy. Naturally the stiffness function must satisfy a (α) < 0, a(0) = 1 and a(1) = 0, and the damage potential function w (α) > 0, w(0) = 0 and up to a rescaling w(1) = 1. The required elastic phase is obtained for α → -a 2 (α)w (α)/a (α) is strictly increasing. The coupling function b (α) < 0 ensure that the damage yield surface decreases with the cumulated plastic strain and b(0) = 1. For numerical reason (a, b, w) must be convex with respect to α which is not the case for the provided closed from solution in [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF] for AT k see Table 5.1. Consequently, we prefer the model named AT 1 where a 1d computed solution example (dark lines) is compared with the numerical simulation (colored lines) see Figure 5.3. The numerical implementation is detailed in the following section 5.2. For this 1d example, we see the three phases described below in the stress-displacement plot, precisely the stress softening leads to a localization of the damage in which a cohesive response is obtained at the center.

Name

a(α) w(α) b(α) AT 1 (1 -α) 2 α a(α) + η b AT k 1 -w(α) 1 + (c 1 -1)w(α) 1 -(1 -α) 2 (1 -w(α)) c 2

Numerical implementation of the gradient damage models coupled with perfect plasticity

In the view to numerically implement the gradient damage model coupled with perfect plasticity it is common to discretized in time and space. For the time discretization evolution we refer to the Definition 8. However in the numerical implantation we do not enforce energy balance condition justified by following the spirit of [START_REF] Bourdin | The variational formulation of brittle fracture: numerical implementation and extensions[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Functions space are discretized through standards finite elements methods over the domain. Both damage and displacement fields are projected over linear Lagrange elements. Whereas the plastic strain tensor is approximated by piecewise constant element. By doing so we probably use the simplest finite element to approximate the evolution problem. Conversely, the chosen finite element space cannot describe the jumps set of u and the localization of p, however it might be possible to account such effects by using instead discontinuous Galerkin methods. Nevertheless, as you will see on numerical simulations performed, the plasticity concentrates in a strip of few elements once the damage localizes. Numerically 5.2. Numerical implementation of the gradient damage models coupled with perfect plasticity we are not restricted to Von Mises plastic criterion only but any associated plasticity. Since a(α), b(α) and w(α) are convex the total energy is separately convex with respect to all variables (u, α, p) but that is not convex. A proposed algorithm to solve the evolution is alternate minimization which guarantees a path decreasing of the energy, but the solution might not be unique. At each time step t i , the minimization for each variables are performed as follows:

i. For a given (α, p) the minimization of E with respect to u is an elastic problem with the prescribed boundary condition. To solve this we employed preconditioned conjugate gradient methods solvers.

ii. The minimization of E with respect to α for fixed (u, p) and subject to irreversibility (α ≤ α i-1 ) is solved using variational inequality solvers provided by PETCs [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF][START_REF] Balay | PETSc users manual[END_REF][START_REF] Balay | PETSc Web page[END_REF].

iii. For a fixed (u, α) the minimization of E with respect to p is not straight forward the raw formulation however reformulated as a constraint optimization problem turns being a plastic strain projection onto a convex set which is solved using SNLP solvers provided by the open source snlp 1 . Boundaries of the stress elastic domain is constrained by a series of yields functions describing the convex set without dealing with none differentiability issues typically corners.

The retained strategy to solve the evolution problem is to use nested loops. The inner loop solves the elasto-plastic problem by alternate i. and iii. until convergence. Then, the outer loop is composed of the previous procedure and ii., the exit is triggered once the damage has converged. This leads to the following Algorithm 4, where δ α and δ p are fixed tolerances. Argument in favor of this strategy is the elasto-plastic is a fast minimization problem, whereas compute ii. is slow, but changing loops orders haven't be tested. All computations were performed using the open source mef90 2 . Verifications of the numerical implementation have been performed on the elastodamage problem and elasto-plasticity problem separately considering three and two dimensions cases. The plasticity is verified with the existence and uniqueness of the bi axial test for elasto-plasticity in [START_REF] Gilles A Francfort | A case study for uniqueness of elasto-plastic evolutions: The bi-axial test[END_REF]. The implementation of the damage have been checked with propagation of fracture in Griffith regime, the optimal damage profile in 2d and many years of development by Bourdin. The verification of the coupling is done by comparison with the one dimensional setting solution in section 5.1.3. Solve the equilibrium,

u k+1 := argmin u∈C i E t i (u, α j , p k ) 6:
Solve the plastic strain projection on each cells,

p k+1 := argmin p∈M n s a(α j )A(e(u k+1 )-p)∈b(α j )K 1 2 A(p -p i-1 ) : (p -p i-1 ) 7: k := k + 1 8: until p k -p k-1 L ∞ ≤ δ p 9:
Set, u j+1 := u k and p j+1 := p k 10:

Compute the damage,

α j+1 := argmin α∈D i α≥α i-1 E t i (u j+1 , α, p j+1 ) 11:
j := j + 1 12: until α jα j-1 L ∞ ≤ δ α 13: Set, u i := u j , α i := α j and p i := p j 5.3 Numerical simulations of ductile fractures

Plane-strain ductility effects on fracture path in rectangular specimens

The model offer a large variety of possible behaviors depending on the choice of functions a(α), b(α), w(α) and the plastic yield function f Y (τ ) considered. From now, the presentation is limited to AT 1 in Table 5.1 and Von Mises plasticity such that,

f Y (σ) = ||σ|| eq -σ p
where ||σ|| eq = n n-1 dev(σ) : dev(σ) and dev(σ) denotes the deviatoric stresses.

Considering an isotropic material, the set of parameters to calibrate is (E, ν, σ p , σ c , ) where the Young's modulus E, the Poisson ratio ν and the plastic yield stress σ p can be easily characterized by experiments. However, σ c and are still not clear but in brittle fracture nucleation they are estimated by performing experiments on notched specimen 5.3. Numerical simulations of ductile fractures see [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. Hence, a parameter analysis for our model is to study influences of the ratio r Y = σ c /σ p and independently.

Consider a rectangular specimen of length (L = 2) and width (H = 1) in plane strain setting, made of a sound material with the set up E = 1, ν = . Let first performed numerical simulations by varying the stress ratio of initial yields surfaces r Y ∈ [. [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF] with an internal length equal to = .02 smaller than the geometric parameters (L, H) and let others parameter unchanged.

The damage fields obtained after failure of samples are summarized on the Figure 5.5. A transition from a straight to a slant fracture for an increasing r Y is observed similarly to the Ti glass alloy in the Figure 5.1. A higher initial yields stress ratio induces a larger plastic strain accumulation leading to a thicker damage localization strip. The measure of the fracture angle reported in Figure 5.5 does not take into account the turning crack path profile around free surfaces caused by the damage condition ∇α • ν = 0. Clearly, for the case σ c < σ p the fracture is straight and there is mostly no accumulation of plastic strain. However due to plasticity, damage is triggered along one of shears bands, resulting of a slant fracture observation in both directions but never two at the same time. Now, let us pick up one of this stress ratio r Y = 5 for instance and vary the internal length ∈ [0.02, 0.2]. The stress vs. displacement is plotted in Figure 5.6 and shows various stress jumps amplitude during the damage localization due to the snap-back intensity. This effect is well known in phase-field models to brittle fracture and pointed out by [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF]. A consequence of this brutal damage localization is a sudden drop of the stress, when this happens the energy balance is not satisfied. Continuous and discontinuous energies evolution is observed for respectively = 0.2 and = 0.02 plotted on Figure 5.7.

The attentive reader may notice that the plastic energy decreases during the damage localization which contradicts the irreversibility hypothesis of the accumulation of dissipated plastic energy. Actually the plotted curve is not accurately representative of the dissipated plasticity energy but it is a combination of damage and plasticity such that a part of this energy is transformed into a surface energy contribution. Hence, those dis- Snap-shots of damage, accumulated plastic strain and damage in a deformed configuration fields are illustrated in Figure 5.9 for different loading time (a, b, c, d) shown in the Figure 5.6. The cumulated plastic strain is concentrated in few mesh elements across the surface of discontinuity (fracture center). Because damage and plasticity evolve together along this strip it is not possible to dissociate mechanism coming from pure plasticity or damage independently. It can be interpreted as a mixture of permanent deformation and voids growing with mutual cause and effects relationship.

Plane-strain simulations on two-dimensional mild notched specimens

In the sequel we restrict our scope to study fractures nucleation and propagation in ductile regime (r Y = σ c /σ p large enough) for a mild notched specimen. Experimentally this design shape samples favor fractures around the smallest cross section size. Necking is a well known instability phenomena during large deformations of a ductile material.

A consequence of the necking on a specimen is a cross sectional reductions which implies a curved profile to the deformed sample. Since we are in small deformations setting, necking cannot be recovered, thus we artificially pre-notch the geometry (sketched in Figure 5.10 with the associated Table 5.2) to recover a plastic strain concentrations. For more realistic numerical simulations and comparisons with pictures of the experiments on Aluminum alloy AA 5754 Al -Mg in Figure 5.2, we set material properties (see Table 5.3) such that the internal length is in the range of grain size, σ c is chosen to recover 7% elongation and (E, ν, σ p ) are given. We assume that the material follows Von Mises perfect plasticity criteria and the elastic stress domain shrinks from σ p to the lower limit of 15% of σ p . The experiments are built such that displacements are controlled at the extremities of the plate and observations are in the sheet thickness direction. Hence, the 2d plane strain theory is adopted for numerical simulations. Also we have studied two types of boundaries conditions, clamped and rollers boundary condition respectively named set-up A and set-up B. .10: Specimen geometry with nominal dimensions, typical mesh are ten times smaller that the one illustrated above. Note that meshes in the center area of the geometry are refined with a constant size h. Also a linear growing characteristic mesh size is employed from the refined area to a coarsen mesh at the boundary. considered, such that, the set-up B provides a slant fracture shear dominating with nucleation at the center and propagation along one of the shear band, and for the set-up A, the fracture nucleates at the center, propagates along the specimen section and bifurcate following shear bands. Final crack patterns are pure shear configuration and a slant-flatslant path. Again some snap shots of damage, cumulated plastic strain and damage in deformed configuration are presented in Figure 5.12 and Figure 5.13 for respectively the set-up A and B. Time loadings highlighted by letter are reported in the stress vs. strain plot in Figure 5.11. Main phenomenon are: (a) during the pure plastic phase there is no damage and the cumulated plastic strain is the sum of two large shear bands where the maximum value is located at the center, (b) the damage is triggered on the middle and develops following shear bands as a "X" shape, (c) a macro fracture nucleates at the center but stiffness remained and the material is not broken, (d) failure of the specimen with the final crack pattern. Close similarities between pictures of ductile fracture nucleations from simulations and experimental observations can be drawn. However, we were not able to capture cup-cones fractures. To recover the desired effect we introduced a perturbation in the geometry such that the parabola shape notch is no more symmetric along the shortest cross section axis, i.e. an eccentricity is introduced by taking ρ < 1 see the Figure 5.10. In a sense there is no reason that necking induces a perfectly symmetric mild notch specimen. Leaving all parameters unchanged and taking ρ = .9 we observed two cracks patterns: a shear dominating and cup-cones for respectively set-up B and set-up A illustrated in Figure 5.14. This type of non-symmetric profile with respect to the shortest cross section axis implies a different stress concentration between the right and the left side of the sample which consequently leads to unbalance the plastic strain concentrations intensity on both parts. Since damage is guided by the dissipated plastic energy we have recovered this cup cones fracture with again a macro fracture has nucleated at the center. Also the set-up B with ρ = .9 is not significantly perturbed to get a new crack path but still in the shear dominating mode.

Ductile fracture in a round notched bar

A strength of the variational approach is that it will require no modification to perform numerical simulations in three dimensions. Also this part is devoted to recover common observations made on ductile fracture in a round notched bar such as cup-cones and shear dominating fractures shapes. The ductile fracture phenomenology for low triaxility (defined as the ratio of the hydrostatic over deviatoric stresses) have been investigated by Benzerga [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF], relevant pictures of cracks nucleation and propagation into a round bar Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage with none destructive techniques is summarized in the Figure 5.19 . Since we focus on the fracture phenomenology we do not attribute physical values to material parameter but give attentions to the yield stress ratio r Y and the internal length . The internal length governs the thickness of the localization which has to be small enough compared to the specimen radius to observe a distinct fracture. In the other sides, drives the characteristics mesh size, typically ∼ 3h which constraint the numerical cost. For clarity the cumulated plastic strain will not be shown anymore since it does not provide further information on the fracture path than the damage. Based on the above results, boundary conditions play a fundamental role in our simulations so we will consider two cases: an eccentric mild notched shape (ρ = .7) specimens in the set-up A and B respectively associated to clamped and rollers boundary conditions. Both geometries are solids of revolution (tensile axis revolution) based on the sketch Figure 5.10 and Table 5 Those simulations were performed with 48 cpus during 48 hours on a 370 000 mesh nodes for 100 time steps with the provided resources of high performance computing of Louisiana State University3 . Results of numerical simulations are shown on the Figures 5.17 The ductile fracture phenomenology is presented by Benzerga-Leblond [START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF], and shows the voids growing and coalescence during the early stage of stress softening, then a macro fracture nucleates at the center end propagates following shear lips formations. Numerical simulations at the loading time (a) for the set-up A and B show a diffuse damage in the middle of the specimen which is exactly a loss of stiffness in the material. This can be interpreted as an homogenization of voids density. A sudden macro crack appears around the loading time (b) which corresponds to the observation made. From (b) to (c) the crack follows shear lips formation in a shear dominating or cup-cones crack patterns depending on the prescribed boundary conditions clamped (set-up A) or rollered (set-up B). These numerical examples suggest that variational phase-field models of ductile fracture are capable of predicting crack nucleation and propagation in low triaxiality specimen for the 2d plane strain specimen and round bar for a simple model considered.

Conclusion

In contrast with most of literature on ductile fracture we proposed a variational model by coupling gradient damage models and perfect plasticity following seminal papers of [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF]. In this chapter, we have investigated crack nucleation and propagation in multiple geometries in simple case of functions under Von Mises perfect plasticity. We confirmed observations reported elsewhere in the literature that fracture nucleates at the center of the specimen and propagates following shear bands before reaching free surfaces for low triaxiality configuration in ductile materials. Our numerical simulations also highlight that crack patterns observed is strongly dependent of the prescribed boundary conditions and geometry which leads to a plastic dissipated energy concentrations path. The strength of the proposed phase-field model is the ability to handle with both ductile and brittle fractures which mostly have been separated like oil and water. The key parameter to capture this transition is the ratio of initial yields surfaces of damage over plastic one. We show that variational phase-field models are capable of qualitative predictions of crack nucleation and propagation in a mild notch range of geometries including two and three dimensions, hence, this model is a good candidate to address the aforementioned issues. Also, the energy balance is preserved since the fracture evolution is smooth driven by and internal length. Of course, there are still many investigations to performed before claiming the superiority of the model such that, fracture nucleation at a notch of a specimen (high triaxiality) which due to the unbounded hydrostatics pressure for the plasticity criteria (Von Mises for instance) leads to hit the damage yield surface first, consequently a brittle response is attended. To get a cohesive response a possible choice of plastic yield surface is to consider a cap model closing the hydrostatic pressure in the stress space domain.

Chapter 6

Concluding, remarks and recommended future work

In this dissertation, we studied the phenomena of fracture in various structures using phase-field models.

The phase-field models have been derived from Francfort Marigo's variational models to fracture which have been conceived as an approximation of Griffith's theory. In Chapter 1 we exposed a complete overview and main properties of the model. In Chapter 2, we applied the phase-field models to study fracture nucleation in a V-and U-notches geometries. Supported by numerous validation we have demonstrated the ability of the model to make quantitative prediction of crack nucleation in mode I. The model is based on general energy minimization principle and does not require any ad-hoc criteria, just to adjust the internal length. Moreover the model properly accounts for size effects that cannot be recovered from Griffith-based theory. In Chapter 3 we have shown that the extended model to hydraulic fracturing satisfies Griffith's propagation criterion and there is no issues to handle with multi-fracking scenario. The fracture path is dictated by the minimization principle of the total energy. A loss of crack symmetry is observed in the case of a pressurized network of parallel fractures. In Chapter 4, we solely focused on the perfect elasto-plasticity models and we started by the classical approach to its variational formulation. A verification of the alternated algorithm technique is exposed. The last chapter was devoted to combine models exposed in the first and the fourth chapter to perform cohesive and ductile fractures. Our numerical simulations have shown the capability of the model to retrieve main features of ductile fractures in a mild notch specimen, precisely nucleation and propagation phenomenon. Also, we have observed that crack paths are sensitive to the geometry and boundary conditions applied on it.

In short, we have demonstrated that variational phase-field models address some of vexing issues associated with brittle fractures: scale effects, nucleation, existence of a critical stress and path prediction. By a simple coupling with the well known perfect plasticity theory, we recovered phenomenology of ductile fractures patterns.

Of course, there are still remaining issues that need to be addressed. Our numerical simulations do not enforce energy balance as indicated by a drop of the total energy upon crack nucleation without string singularities illustrated in Chapter 2. Perhaps extensions into phase field models of dynamic fracture will address this issue.

Also fracture in compression remains an issue in variational phase-field models. It is not clear of either of this models is capable of simultaneously accounting for nucleation under compression and self-contact.

A recommended future work is to study ductile fractures following the spirit of Chapter 2. The idea is by varying the yields stress ratio recover first the brittle initiation criterion and then study the ductile fracture for different notch angles.

  Primal feasibility l ≥ l(t) 3. Dual feasibility λ ≥ 0 4. Complementary slackness λ(ll(t)) = 0

10 )

 10 Splitting the second integral over ∂Ω = (∂ N Ω ∪ ∂ D Ω) \ Γ(l) into the Dirichlet and the remainder Neumann boundary part and using the condition v = 0 on ∂ D Ω, the Gateaux derivative of E becomes,

Figure 1 . 1 :

 11 Figure 1.1: Sketch on the left shows the evolution of the crack (red curve) for a strict decreasing function G(1, l) subject to the irreversibility (l ≥ l 0 ) and G(1, l) ≤ G c /t 2 . Picture on the right shows the crack evolution for a local minimality principle (red curve) and for a global minimality (blue curve) without taking into account the energy balance.

w 1 . 3 .Remark 1

 131 Limit of the damage energy The above expression (1.68) is invariant by a change of variable x = x, thus β(x) = β( x)

Chapter 1 .

 1 Combining this two statements, we deduce that there exists a , b , c in I such that, a ≤ b ≤ c , and lim →0 α (a ) = lim →0 α (c ) = 0 and lim →0 α (b ) = 1 thus, I w(α ) + (α ) 2 dx = b a w(α ) + (α ) 2 dx + c b w(α ) + (α ) 2 dx (1.78) Again using the identity, a 2 + b 2 ≥ 2|ab|, we have that, Variational phase-field models of brittle fracture where (x) := x 0 w(s)ds. Using the substitution rule we then get, b a w(α ) + (α ) 2 dx ≥ 2 | (b ) -(a )| , (1.80) and since (0) = 0 and (1) = c w , we obtain, ) + (α ) 2 dx ≥ 2c w ,

  82) and α = 1 for |x| ≤ b , we get that, b -b

0 H n- 1 2 x∈ΩChapter 1 .

 0121 ({x; d(x) = y}) dy = 2H n-1 (J(u)) (1.100) and for the second term, b ≤d(x)≤δ w(α (d(x))) + |∇α (d(x))| 2 dx = w(α (d(x))) + α (d(x))∇d(x) 2 dH n-1 (x) dy = δ b x∈Ω w(α (y)) + α (y)∇d(x) 2 dH n-1 ({x; d(x) = y}) dy = δ b w(α (y)) + α (y) dH n-1 ({x; d(x) = y}) dy (1.101) Making the change of variable y = x , Variational phase-field models of brittle fracture b ≤d(x)≤δ w(α (d(x))) + |∇α (d(x))| 2 dx = δ b w(α (y)) + α (y) 2 H n-1 ({x; d(x) = y}) dy = δ b w(α (y)) + α (y) 2 s (y) dy

) 1 . 4 .Definition 5 (

 145 Numerical implementationand the discrete time evolution problem is given by, Damage discrete evolution by local minimizers)

  Figure 2.1(left) shows the outcome of a surfing experiment on a rectangular domain Ω = [0, 5] × [-1 2 , 1 2

Figure 2 . 1 :

 21 Figure 2.1: Mode-I "surfing" experiment along straight (left) and circular (right) paths. Dependence of the crack length and elastic energy release rate on the loading parameter for multiple values of .

Figure 2 . 2 :

 22 Figure 2.2: Pac-man geometry for the study of the crack nucleation at a notch. Left: sketch of the domain and notation. Right: relation between the exponent of the singularity λ and the notch opening angle ω determined by the solution of equation (2.10). For any opening angle ω we apply on ∂ D Ω the displacement boundary condition obtained by evaluating on ∂ D Ω the asymptotic displacement (2.12) with λ = λ(ω).

  The mode-I Pac-Man test Consider a Pac-Man-shaped 3 domain with radius L and notch angle ω as in Figure 2.2(left). In linear elasticity, a displacement field associated with the stress field (2.7) is

Figure 2 . 3 :σ

 23 Figure 2.3: Pac-Man test with the AT 1 model, L = 1, = 0.015, ω = 0.7π, and ν = 0.3.From left to right: typical mesh (with element size ten times larger than that in typical simulation for illustration purpose), damage field immediately before and after the nucleation of a crack, and plot of the energies versus the loading parameter t. Note the small damaged zone ahead of the notch tip before crack nucleation, and the energetic signature of a nucleation event.

Figure 2 . 4 :

 24 Figure 2.4: Identification of the generalized stress intensity factor: σ θθ (r,0) (2π r) λ-1 along the domain symmetry axis for the AT 1 (left) and AT 2 (right) models with undamaged notch conditions, and sub-critical loadings. The notch aperture is ω = π/10

Figure 2 . 5 :

 25 Figure 2.5: Critical generalized critical stress intensity factor at crack nucleation as a function of the internal length for ω 0 (left) and ω π/2 (right). AT 1 -U, AT 1 -D, AT 2 -U, and AT 2 -D refer respectively to computations using the AT 1 model with damaged notch and undamaged notch boundary conditions, and the AT 2 model with damaged notch and undamaged notch boundary conditions. (K Ic ) eff := G eff E 1-ν 2 denotes the critical mode-I stress intensity factor modified to account for the effective toughness G eff .

Figure 2 . 6 :

 26 Figure 2.6: Critical generalized stress intensity factor k for crack nucleation at a notch as a function of the notch opening angle ω. Results for the AT 1 and AT 2 models with damaged -D and undamaged -U notch lips conditions. The results are obtained with numerical simulations on the Pac-Man geometry with (K Ic ) eff = 1 and = 0.01 so that σ c = 10 under plane-strain conditions with a unit Young's modulus and a Poisson ratio ν = 0.3.

  .m 1-λ ]

Figure 2 . 7 :

 27 Figure 2.7: Critical generalized stress intensity factor k c vs notch angle.Comparison between numerical simulations with the AT 1 and AT 2 models and damaged and undamaged boundary conditions on the notch edges with experiments in steel from[START_REF] Strandberg | Fracture at V-notches with contained plasticity[END_REF] (top-left), and Duraluminium (top-right) and PMMA (bottom) from[START_REF] Seweryn | Brittle fracture criterion for structures with sharp notches[END_REF].

Figure 2 . 8 :

 28 Figure 2.8: Critical generalized stress intensity factor k c vs notch angle and depth in PVC foam samples from [94]. Numerical simulations with the AT 1 model with damaged and undamaged notch conditions (left), and AT 2 model with damaged and undamaged notch conditions (right).

Figure 2 . 10 :

 210 Figure 2.10: Critical generalized stress intensity factor k c vs notch angle for Al 2 O 3 -7%ZrO 2 (left) and PMMA (right). The black markers represents all experimental results. The numerical results are obtained through the Pac-Man test using the AT 1 model. See Tables 2.8-2.9 in the Appendix B for the raw data.

Figure 2 .

 2 Figure 2.12: DENT geometry

2. 3 .AT 2 -

 32 Size effects in variational phase-field models 1 -U ρ = 0.5AT 1 -U ρ = 1.25 AT 1 -U ρ = 2.5 U ρ = 0.5 AT 2 -U ρ = 1.25 AT 2 -U ρ = 2.5

Figure 2 . 13 :h = ρ 2 ah h = R 100 Figure 2 . 14 :

 213100214 Figure 2.13: Crack nucleation at U-notches. Comparison between experimental data of [92] and numerical simulations using the AT 1 (top) and AT 2 (bottom) models.

Figure 2 . 15 :

 215 Figure 2.15: Damage field at the boundary of the hole in the elastic phase 0 < t < t e (left), the phase with partial damage t e < t < t c (center), and after the nucleation of a crack t > t c (right). Blue: α = 0, red: α = 1. The simulation is for ρ = 1.0 and a/ = 5.

Figure 2 . 16 :

 216 Figure 2.16: Normalized applied macroscopic stress t e /σ c at damage initiation as a function of the aspect ratio ρ for a/ = 1 (left) and of the relative defect sizes a/ for ρ = 1 and ρ = 0.1 (right).

52 2. 3 .Figure 2 . 17 :

 523217 Figure2.17: Normalized applied macroscopic stress t c /σ e at crack nucleation for an elliptic cavity in an infinite plate. Left: shape effect for cavities of size much larger than the internal length (a/ = 48); the solid line is the macroscopic stress at the damage initiation t e (see also Figure2.16) and dots are the numerical results for the AT 1 model. Right: size effect for circular (ρ = 1.0) and highly elongated (ρ = 0.1) cavities.

Figure 2 . 53 Chapter 2 .Figure 2 . 18 :

 2532218 Figure 2.18: Initiation of a crack of length 2a in a plate of finite width 2W . The numerical results (dots) are obtained with the AT 1 model for = W/25. The strength criterion and the Griffith's criterion (2.18).

Figure 3 . 1 :

 31 Figure 3.1: Sketch of the geometry (invariant). The symmetry axis being a reflection for 2d and a revolution axis in 3d.

3. 2 .

 2 Numerical verification case of a pressurized single fracture in a two and three dimensions

Chapter 3 .Figure 3 . 2 :

 332 Figure 3.2: Evolutions of normalized p, V and l for the line fracture (left column figures) and penny shape crack (right column figures). Colored dots refer to numerical results and solid black lines to the closed form solution given in Appendix C. For the line fracture, V c = 4πl 3 0 (G c ) eff /E and p c = E (G c ) eff /(πl 0 ), where E = E/(1ν 2 ) in plane strain theory and E = E in plane stress. For the penny shape crack, V c = 8/3 πl 5 0 (G c ) eff /E and p c = πE(G c ) eff /(4l 0 ).

3. 2 .Figure 3 . 3 :

 233 Figure 3.3: Snap-shots of damage for the line fracture example at different loadings, such that, before the loading cycle (top), before refilling the fracture (middle) and during the propagation (bottom). The red color is fully damage material and blue undamaged. We see the casing mesh which encapsulates the fracture.

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Snap shots (view from above) of fracture damage (α ≥ .99) for the penny shape crack example at different loadings, that is before refilling the fracture (left) and during the propagation (right). The solid black lines are the limit of the casing.

Figure 3 . 6 :

 36 Figure 3.6: Infinite network of parallel cracks domain (left). Domain duplications form the smallest invariant domain (right).

Figure 3 . 7 :

 37 Figure 3.7: Domains in the deformed configuration for respectively Ω 1 , Ω 2 , Ω 4 and Ω 6 .The pseudo-color blue is for undamaged material and turns white when (α ≤ .01) (visibility reason). The full colors correspond to numerical simulations cells domain see table3.2, and opacity color refers to the rebuild solution using symmetries. In all simulations only one crack propagates in the domain. Using the multiplicity pictures from left to right we obtain a fracture propagation periodicity denoted period. of 6/6, 3/6, 1.5/6 and 1/6.

6 Figure 3 . 8 :

 638 Figure 3.8: Plots of normalized variables such that crack pressure, average fracture length and energy density (per Ω) vs. fluid volume density (per Ω) respectively on the (top-left) and (top-right) and (bottom-right). The aperture of the longest crack for 2V /(nV c ) = 13.Colored plot are numerical results for different domain sizes Ω 1 , Ω 2 , Ω 4 and Ω 6 . The solid black line is the closed form solution and gray the approximated solution given by Sneddon[START_REF] Sneddon | Crack problems in the classical theory of elasticity[END_REF].

Chapter 3 . 6 Figure 3 . 10 :

 36310 Figure 3.10: Ratio of critical pressures (multi-fracking over single fracture) vs. the inverse of the fracture density (hight density on the left x-axis and low density on the right side).Black dash line is r p (ρ/n) with 1/n the periodicity. Colored line is numerical results for respectively a periodicity 6/6, 3/6 and 1.5/6.

Figure 3 . 11 :

 311 Figure 3.11: Fracture toughness vs. confining pressure for the Indiana limestone

Chapter 3 .LFigure 3 . 12 :

 3312 Figure 3.12: Schematic of burst experiment for jacketed bore on the (left). Pre-(middle) and Post-(right) burst experiment photos.

Figure 3 . 13 :

 313 Figure 3.13: Rigorous superposition of the burst problem.L

Figure 3 . 14 :

 314 Figure 3.14: Superposition of the burst problem applied in Abou-Sayed (1978).

Figure 3 . 15 :

 315 Figure 3.15: Comparison of the normalized stress intensity factor for the jacketed and unjacketed problems receptively denoted K J * I and K U * I vs. the normalized crack length l. Numerical computational SIF based on G θ method (colored lines) overlay plots provided by Clifton in [53].

10 Figure 3 . 16 :

 10316 Figure 3.16: Computed normalized SIF vs. normalized crack length l for two confining pressure ratios r = 1/8 (dash lines) and r = 1/6 (solid lines) and various w = {3, 4, 7, 10} (colored lines).

Figure 3 . 17 :

 317 Figure 3.17: Three possible regime for K B * I denoted (a), (b) and (c). l U S is a critical point from unstable to stable crack propagation, vice versa for l SU . The fracture does not propagates at stop point denoted l ST

Figure 3 . 18 :

 318 Figure 3.18: Computed normalized SIF vs. normalized crack length for the unconfined (left) and confined (right) burst experiments according to the Table 3.3.

5 r = 1 /6 w = 5 Figure 3 . 19 :

 515319 Figure 3.19: Colored lines are computed normalized SIF vs. normalized crack length for unstable propagation (l 0 ≥ .5). Red markers are time step results obtained using the phase-field model.

(4. 16 )

 16 By applying the supremum function for all τ * ∈ K, it comes that, for τ ∈ K, τ : ṗ ≥ sup τ * ∈K {τ * : ṗ}.(4.17)

4. 2 .

 2 Variational formulation of perfect plasticity models q = p + hp, ∀p ∈ M n s Plug this into (4.27) and send h → 0, then using the definition of Gateaux derivative and the subgradient, the stability definition leads to τ = -∂ ψ ∂p (e, p) = lim h→0 -ψ(e, p + hp) -ψ(e, p) hp ∈ ∂H(pp t i-1 ). (4.28)

30 )τ

 30 For any given τ * a fixed element of K, * : p(s) ds = τ * : (qp).

Chapter 4 .Algorithm 3 1 :

 431 Variational models of perfect plasticity as a constraint optimization problem implemented using SNLP solvers provided by the open source snlp 1 . All computations were performed using the open source mef90 2 . Elasto-plasticity alternate minimization algorithm for the step i Let j = 0 and p 0 := p i-1 2: repeat 3:

  p i-1 ) : (pp i-1 ) 5:

  line segment [( d/2, l/2, 0), (d/2, l/2, l)]

Figure 4 . 1 :

 41 Figure 4.1: The closed-form solution equations (4.37),(4.38) are denoted in solid-lines and dots referred to numerical results is in dots. (Top-left) and (-right) figures show respectively the hydrostatics evolution of stresses and plastic strains with the time loading. The figure in the bottom shows displacements for t = 2.857 along the lineout axis [(-d/2, -l/2, 0) × (d/2, l/2, l)]

  3.1, 5.3.2 and 5.3.3. 5.1 Phase-field models to fractures from brittle to ductile 5.1.1 Experimental observations of ductile fractures It is common to separate fractures into two categories; brittle and ductile fractures with different mechanisms. However relevant experiments [110] on Titanium alloys glass show a transition from brittle to ductile fractures response (see Figure 5.1) by varying only one parameter: the concentration of Vanadium. Depending on the Vanadium quantity, they observed a brutal formation of a straight crack, signature of brittle material response for low concentrations. Conversely a smooth stress softening plateau is measured before failure for higher concentrations. The post mortem samples show a shear dominating fracture characteristic of ductile behaviors. / Plot of uniaxial tension test data with optical images of dogbone specimens post-failure. (top eering stress as a function of engineering strain is plotted from tensile tests on dogbone samples o series composites. Samples were loaded until failure at a constant strain rate of 0.2 mm/min. Curv n the x-axis to highlight differences in plastic deformation behavior between alloys. (top right) h of complete V0 dogbone sample after failure in tension. (bottom) Optical microscope images at ilure in deformed alloys V2-V10 and DV1.

Figure 5 . 1 :

 51 Figure 5.1: Pictures produced by [110] show post failure stretched specimens of Ti-based alloys V x --Ti 53 -x/2 Zr 27 -x/2 Cu 5 Be 15 V x . From left to right: transition from brittle to ductile with a concentration of Vanadium respectively equal to 2%, 6% and 12%.

  α)A(e(u)p) Plasticity occurs in the material once the stress reaches a critical value defined by the plastic yield function f Y : M n×n s → R convex such that f Y (0) < 0. We proposed to couple the damage with the admissible stress set through the coupling function b(α) such that, the stress is constrained by σ ∈ b(α)K, where K := {τ ∈ M n×n s s.t. f Y (τ ) ≤ 0} is a non empty close convex set. The elastic stress domain is subject to isotropic transformations by b(α) a state function of the damage. Naturally to recover a stress-softening response, Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage the coupling function b(α) is continuous decreasing such that b(0) = 1 and b(1) = η b , where η b is a residual. By considering associated plasticity the plastic potential is equal to the yield function and the plastic flow occurs once the stress hits the yield surface, i.e. σ ∈ b(α)∂K. At this moment, the plastic evolution is driven by the normality rule such that the plastic flow lies in the subdifferential of the indicator function denoted I at σ, written as, ṗ ∈ ∂I b(α)K (σ)

.34) 6 . 7 . 2 c( 5 . 36 ) 8 .

 6725368 Plastic flow rule in the bulk, b(α t (x))σ p | ṗt (x)|σ t ṗt (x) = 0(5.35) The damage consistency in the bulk, jump and boundary,f D (α t (x), p t (x), pt (x)) αt (x) = 0 b (α t (x 0 )) u(x 0 ) σ p -2 σ kE α t (x 0 ) αt (x) = 0 α t (±L) αt (±L) = 0The energy balance at the boundary, α t (±L) αt (±L) = 0 (5.37)9. The irreversibility which applies everywhere in Ω, 0 ≤ α t (x) ≤ 1, αt (x) ≥ 0 (5.38)

Chapter 5 .Figure 5 . 3 :

 553 Figure 5.3: Comparisons of the computed solution (dark lines) for AT 1 see Table5.1 with the numerical simulation (colored lines) for parameters E = 1, σ p = 1, = 0.15, σ c = 1.58, L = .5 and η b = 0. The (top-left) picture shows the stress-displacement evolution, the (top-right) plot is the displacement jump vs. the stress during the softening behavior. The (bottom-left) figure shows the damage profile during the localization for three different loadings. The (bottom-right) is the evolution of the energy release vs. the displacement jump also known as the cohesive law(Barenblatt).

Chapter 5 .Algorithm 4 1 :: repeat 3 :

 5413 Variational phase-field models of ductile fracture by coupling plasticity with damage Alternate minimization algorithm at the step i Let, j = 0 and α 0 := α i-1 , p 0 := p i-1 2Let, k = 0 and p 0 := p j

3 Figure 5 . 4 :

 354 Figure 5.4: Rectangular specimen in tensile with rollers boundary condition on the leftright extremities and stress free on the remainder. The characteristic mesh size is h = /5.

Figure 5 . 5 :

 55 Figure 5.5: Shows fracture path angle vs. the initial yields stress ratio r Y . Transition form straight to slant crack characteristic of a brittle -ductile fracture transition.

Figure 5 . 6 : 9 . 5 . 3 .Figure 5 . 7 :Figure 5 . 8 :

 569535758 Figure 5.6: Stress vs. displacement plot for σ c /σ p = 5, shows the influence of the internal length on the stress jump amplitude signature of the snap back intensity. Letters on the curve = .1 referees to loading times when snap-shots of α and p are illustrated in Figure 5.9.

123 Chapter 5 .Figure 5 . 9 :

 123559 Figure 5.9: Rectangular stretched specimen with rollers boundary displacement for parameters σ c /σ p = 5 and = .1, showing snap-shots of damage, cumulated plastic strain and damage in deformed configuration (the displacement magnitude is 1%) at different loading time refereed to the plot 5.6 for (a, b, c, d). The cumulated plastic strain defined as p = t 0 || ṗ(s)||ds has a piecewise linear color table with two pieces, [0, 14] for the homogeneous state and [14, 600] for visibility during the localization process. Moreover the maximum value is saturated.

Figure 5

 5 Figure 5.10: Specimen geometry with nominal dimensions, typical mesh are ten times smaller that the one illustrated above. Note that meshes in the center area of the geometry are refined with a constant size h. Also a linear growing characteristic mesh size is employed from the refined area to a coarsen mesh at the boundary.

Figure 5 . 11 :

 511 Figure 5.11: Plot of the stress vs. strain (tensile axis component) for the mild notch specimen with clamped and rollers interfaces conditions respectively set-up A and set-up B

126 5. 3 .Figure 5 . 12 :Figure 5 . 13 :

 1263512513 Figure 5.12: Zoom in the center of mild notched stretched specimen with clamped boundary displacement (set-up A) showing snap-shots of damage, cumulated plastic strain and damage in deformed configuration (the displacement magnitude is 1) at different loading time refers to Figure 5.11 for (a, b, c, d). The cumulated plastic strain color table is piecewise linear with two pieces, [0, .35] for the homogeneous state and [.35, 2.5] for visibility during the localization process. Moreover the maximum value is saturated. The pseudo color turns white when (α ≥ 0.995) for the damage on the deformed configuration figure.

5. 3 .Figure 5 . 14 :

 3514 Figure 5.14: Zoom in the center of eccentric mild notched stretched specimen (ρ = .9) showing snap-shots of damage, cumulated plastic strain and damage in deformed configuration (the displacement magnitude is 1) at the failure loading time, for the set-up A and B. The cumulated plastic strain color table is piecewise linear with two pieces, [0, .35] for the homogeneous state and [.35, 2.5] for visibility during the localization process. Moreover the maximum value is saturated. The pseudo color turns white when (α ≥ 0.995) for the damage on the deformed configuration figure.

  and 5.18 were fractures patterns are similar to one observed in the literature see pictures 5.15 and 5.16. An overview of the fracture evolution in round bar are exposed in the Figures 5.19

  .

Figure 5 . 15 :

 515 Figure 5.15: Photo produced by [107] showing cup cones fracture in a post mortem rounded bar.

Figure 5 . 16 :

 516 Figure 5.16: Photo produced by [107] showing shear dominating fracture in a post mortem rounded bar.

Figure 5 .

 5 Figure 5.17: Snap-shot of the damage in deformed configuration for the set-up A after failure, two pieces next to each other.

Figure 5 .

 5 Figure 5.18: Snap-shot of the damage in deformed configuration for the set-up B after failure, two pieces next to each other.

Chapter 5 .Figure 5 . 19 :

 5519 Figure 5.19: Picture in Benzerga-Leblond[START_REF] Benzerga | Ductile fracture by void growth to coalescence[END_REF] shows the phenomenology of ductile fracture in round notched bars of high strength steel: damage accumulation, initiation of macroscopic crack, crack growth and shear lip formation. Numerical simulations shows the overlapped stress vs. displacement blue and orange curves for respectively set-up A and setup B, and snap shots of damage slices in the deformed round bar. The hot color table illustrates the damage, the red color turns white for α ≥ 0.95 which correspond to less than 0.25% of stiffness.

  2.11: Critical load in the three-and four-point bending experiments of a Al 2 O 3 -7%ZrO 2 sample (left) and four-point bending of a PMMA sample (right) from[START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] compared with numerical simulations using the AT 1 model and undamaged notch boundary conditions. Due to significant variations in measurements in the first set of experiments, each data point reported in[START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] is plotted. For the PMMA experiments, average values are plotted. See Table2.10 and 2.11 in the Appendix B for raw data.

Table 2 .

 2 3: Critical generalized stress intensity factor k for crack nucleation at a notch as a function of the notch opening angle ω from Figure2.5. Results for the AT 1 and AT 2 models with damaged -D and undamaged -U notch lips conditions. The results are obtained with numerical simulations on the Pac-Man geometry with (K Ic ) eff = 1 and = 0.01 so that σ c = 10 under plane-strain conditions with a unit Young's modulus and a Poisson ratio ν = 0.3.

	500 1.292	1.084	1.349	1.284
	10 0°0.500 1.308	1.091	1.328	1.273
	20 0°0.503 1.281	1.121	1.376	1.275
	30 0°0.512 1.359	1.186	1.397	1.284
	40 0°0.530 1.432	1.306	1.506	1.402
	50 0°0.563 1.636	1.540	1.720	1.635
	60 0°0.616 2.088	1.956	2.177	2.123
	70 0°0.697 2.955	2.704	3.287	3.194
	80 0°0.819 4.878	4.391	5.629	5.531
	85 0°0.900 6.789	5.890	7.643	7.761
	89 9°0.998 9.853	8.501	9.936	9.934

Table 2 .

 2 4: Generalized critical stress intensity factors as a function of the notch aperture in soft annealed tool steel, (AISI O1 at -50 • C). Experimental measurements from[START_REF] Strandberg | Fracture at V-notches with contained plasticity[END_REF] using SENT and TPB compared with Pac-Man simulations with the AT 1 model.

			Experiments Undamaged notch	Damaged notch
	2ω	Mat	k	(exp) c	stdev k c (num)	rel. error k	(num) c	rel. error
	0°H80	0.14	0.01	0.18	22.91 %	0.15	5.81 %
		H100 0.26	0.02	0.34	24.62 %	0.28	7.61 %
		H130 0.34	0.01	0.44	29.34 %	0.36	5.09 %
		H200 0.57	0.02	0.74	47.60 %	0.61	6.53 %
	90°H80	0.20	0.02	0.22	12.65 %	0.21	4.73 %
		H100 0.36	0.02	0.41	12.29 %	0.38	4.10 %
		H130 0.49	0.05	0.54	11.33 %	0.50	0.50 %
		H200 0.81	0.08	0.91	20.54 %	0.83	2.21 %
	140°H80	0.53	0.06	0.53	0.37 %	0.48	9.26 %
		H100 0.89	0.04	0.92	3.43 %	0.84	5.91 %
		H130 1.22	0.10	1.25	2.95 %	1.13	7.48 %
		H200 2.02	0.14	2.07	4.92 %	1.89	6.80 %
	155°H80	0.86	0.07	0.83	3.63 %	0.75	14.36 %
		H100 1.42	0.08	1.42	0.14 %	1.29	10.63 %
		H130 1.90	0.10	1.95	2.82 %	1.76	8.06 %
		H200 3.24	0.15	3.23	0.89 %	2.92	11.02 %
	Table 2.5: Generalized critical stress intensity factors as a function of the notch aper-
	ture in Divinycell® PVC foam. Experimental measurements from [94] using four point
	bending compared with Pac-Man simulations with the AT 1 model.

Table 2 .

 2 6: Generalized critical stress intensity factors as a function of the notch aperture in Duraluminium. Experimental measurements from[START_REF] Seweryn | Brittle fracture criterion for structures with sharp notches[END_REF] using single edge notch tension compared with Pac-Man simulations with the AT 1 model.

			Experiments Undamaged notch	Damaged notch
	ω	Type	k c (exp)	stdev k c (num)	rel. error k	(num) c	rel. error
	10°DENT 1.87	0.03	2.50	25.29 %	2.07	10.03 %
	20°DENT 1.85	0.03	2.53	26.89 %	2.13	12.97 %
	30°DENT 2.17	0.03	2.65	18.17 %	2.33	6.92 %
	40°DENT 2.44	0.02	3.07	20.65 %	2.73	10.70 %
	50°DENT 3.06	0.05	3.94	22.31 %	3.54	13.63 %
	60°DENT 4.35	0.18	5.95	26.97 %	5.41	19.69 %
	70°DENT 8.86	0.18	11.18	20.74 %	10.10	12.26 %
	80°DENT 28.62	0.68	27.73	3.20 %	24.55	16.56 %
	90°DENT 104.85 10.82 96.99	8.11 %	85.37	22.82 %
	Table 2.7: Generalized critical stress intensity factors as a function of the notch aper-
	ture in PMMA. Experimental measurements from [165] using single edge notch tension
	compared with Pac-Man simulations with the AT 1 model.

Table 2 .

 2 8: Generalized critical stress intensity factors as a function of the notch aperture in Aluminium oxide ceramics. Experimental measurements from[START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] using three and four point bending compared with Pac-Man simulations.

		Experiments Undamaged notch	Damaged notch
	2ω	a/h k c (exp)	stdev k c (num)	rel. error k c (num)	rel. error
	60°0.1 1.41	0.02	1.47	4.5%	1.29	9.3%
		0.2 1.47	0.04	1.47	0.4%	1.29	14.0%
		0.3 1.28	0.03	1.47	13.0%	1.29	0.4%
		0.4 1.39	0.04	1.47	5.8%	1.29	7.8%
	90°0.1 2.04	0.02	1.98	3.0%	1.81	12.9%
		0.2 1.98	0.01	1.98	0.0%	1.81	9.6%
		0.3 2.08	0.03	1.98	5.1%	1.81	15.2%
		0.4 2.10	0.03	1.98	5.9%	1.81	16.1%
	120°0.1 4.15	0.02	3.87	7.3%	3.63	14.3%
		0.2 4.03	0.06	3.87	4.2%	3.63	11.0%
		0.3 3.92	0.18	3.87	1.4%	3.63	8.0%
		0.4 3.36	0.09	3.87	13.0%	3.63	7.4%
	Table 2.9: Generalized critical stress intensity factors as a function of the notch aperture
	in PMMA. Experimental measurements from [71] using three and four point bending
	compared with Pac-Man simulations. The value a/h refers to the ratio depth of the
	notch over sample thickness. See Figure 2.9 for geometry and loading.

Table 2 .

 2 10: Critical load reported in[START_REF] Yosibash | Failure criteria for brittle elastic materials[END_REF] using three-and four-point bending experiments of an Al 2 O 3 -7%ZrO 2 sample compared with numerical simulations using the AT 1 model and undamaged notch boundary conditions. TPB and FPB refer respectively to three point bending and four point bending. See Figure2.9 for geometry and loading.

	2ω	a/h P c (exp)	[N] stdev P	(num) c	[N] rel. error
	60°0.1 608.50	6.69	630.81	3.5%
		0.2 455.75	12.48 451.51	0.9%
		0.3 309.00	8.19	347.98	11.2%
		0.4 258.75	6.61	268.69	3.7%
	90°0.1 687.33	5.19	668.69	2.8%
		0.2 491.00	2.94	491.41	0.1%
		0.3 404.33	5.44	383.33	5.5%
		0.4 316.00	4.24	297.48	6.2%
	120°0.1 881.75	4.60	822.22	7.2%
		0.2 657.25	9.36	632.32	3.9%
		0.3 499.60	25.41 499.50	0.0%
		0.4 336.25	9.09	386.87	13.1%

Table 2 .

 2

11: Load at failure reported in

Table 3 . 3

 33

		288 10.07 1310	0	0.218	365
	Id 1	0.279 8.93	9775 1/8 0.2025	1462
	Id 2	0.258 9.65 14907 1/8 0.2060	1954
	Id 3	0.273 9.12 11282 1/6 0.2128	1023
	Id 4	0.283 8.82 17357 1/6 0.1102	2550
	Id 5	0.257 9.70 18258 1/6 0.2022	1508

: Rock specimen dimensions provided by the commercial laboratory and calculated fracture toughness.

  where A is the Hooke's law tensor. The domain is subject to time dependent stress boundary condition σ • ν = g(t) on ∂ N Ω. A safe load condition g(t) is prescribed to prevent issues in plasticity theory. The total energy is formulated for every x ∈ Ω and every t by E t (u, p) =

	4.3.1 Numerical implementation of perfect plasticity models
	Consider the same problem with stress conditions at the boundary and a free energy of
	the form of,	ψ(e(u), p) =	1 2	A(e(u) -p) : (e(u) -p),
	Ω	1 2	A(e(u) -p) : (e(u) -p) +	0	t	sup τ ∈K {τ : ṗ(s)}ds dx

  e 2 + νg 2 e 3 ⊗ e 3 + tEe 3 ⊗ e 3 ⊗ e 2 + t(-νe 1 ⊗ e 1νe 2 ⊗ e 2 + e 3 ⊗ e 3 ) x 2 e 2 + t(-νx 1 e 1νx 2 e 2 + x 3 e 3 )

	e(t) = -ν(1 + ν) e 2 u(t) = -g 2 E e 1 ⊗ e 1 + (1 -ν 2 ) g 2 E ν(1 + ν) E g 2 x 1 e 1 + 1 -ν 2 E g 2 (4.37)
	After the critical loading, permanent deformation takes place in the structure and the
	solution is

Table 5 .

 5

1: Variety of possible models, where c 1 , c 2 are constants.

Table 5 .

 5 2: Specimen dimensions. All measures are in [mm]. The internal length is specified in Table 5.3. We observed two patterns of ductile fractures depending on the boundary condition 125 Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage

	E	ν	σ p	σ c
	[GPa]		[MPa] [GPa] [µm]
	70	.33	100	2	400

Table 5 .

 5 3: Material parameters used for AA 5754 Al -Mg.

Table 5 .

 5 .[START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF]. The smallest damageable plastic yield surface is given for 5% of σ p . 4: Specimen dimensions. For the internal length refer to the Table5.5.

	L	H	W	r	D	d	l	h
	4.5 2.2 1.05 .5 1.09 0.98 0.82	/2.5
			E ν σ p r Y
			1 .3 1	12 .03

Table 5 .

 5 5: Parameters used for 3d simulations.

Karush-Kuhn-Tucker

available at https://www.bitbucket.org/bourdin/mef90-sieve

available at https://www.bitbucket.org/bourdin/mef90-sieve

available at https://bitbucket.org/cmaurini/gradient-damage

https://en.wikipedia.org/wiki/Pac-Man

available at https://www.bitbucket.org/bourdin/mef90-sieve

available at http://abs-5.me.washington.edu/snlp/ and at https://bitbucket.org/bourdin/ snlp

available at https://www.bitbucket.org/bourdin/mef90-sieve

http://www.hpc.lsu.edu

Remerciements

problem. Perhaps extensions into phase field models of dynamic fracture will address this issue.

Fracture in compression remains an issue in variational phase-field models. Although several approaches have been proposed that typically consist in splitting the strain energy into a damage inducing and non damage inducing terms, neither of the proposed splits are fully satisfying (see [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Lancioni | The variational approach to fracture: A practical application to the french Panthéon[END_REF][START_REF] Li | Gradient Damage Modeling of Dynamic Brittle Fracture[END_REF] for instance). In particular, it is not clear if either of this models is capable of simultaneously accounting for nucleation under compression and self-contact.

Finally, even though a significant amount of work has already been invested in extending the scope of phase-field models of fracture beyond perfectly brittle materials, to our knowledge, none of the proposed extensions has demonstrated its predictive power yet.

Appendix C

Single fracture in a infinite domain Line Fracture (2d domain):

The volume of a line fracture in a 2d domain is

where E = E/(1ν 2 ) in plane strain and E = E in plane stress theory. Before the start of propagation, l = l 0 and the fluid pressure in this regime is

If we consider an existing line fracture with an initial length of l 0 . Prior to fracture propagation, the fracture length does not change so that l = l 0 . Since fracture length at the onset of propagation is l 0 , the critical fluid pressure [START_REF] Sneddon | The opening of a griffith crack under internal pressure[END_REF] is

The critical fracture volume at the critical fluid pressure is obtained by substituting (3.16) into (3.14)

During quasi-static propagation of the fracture, l ≥ l 0 and the fracture is always in a critical state so that (3.16) applies. Therefore, the fluid pressure and fracture length in this regime are

Chapter 4

Variational models of perfect plasticity

Elasto-plasticity is a branch of solid mechanics which deals with permanent deformation in a structure once the stress reached a critical value at a macroscopic level. This topic is a vast research area and it is impossible to cover all contributions. We will focus on recalling basic mathematical and numerical aspects of perfect elasto-plasticity in small strain theory under quasi-static evolution problems. The perfect elasto-plastic materials fall into the theory of generalized standard materials developed by [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF][START_REF] Suquet | Sur les équations de la plasticité: existence et régularité des solutions[END_REF][START_REF] Salençon | Elasto-plasticité[END_REF][START_REF] Marigo | From clausius-duhem and drucker-ilyushin inequalities to standard materials[END_REF][START_REF] Mielke | A mathematical framework for generalized standard materials in the rate-independent case[END_REF].

Recently, a modern formalism of perfect plasticity arose [START_REF] Dal Maso | Quasistatic evolution problems for linearly elastic-perfectly plastic materials[END_REF][START_REF] Solombrino | Quasistatic evolution problems for nonhomogeneous elastic plastic materials[END_REF][START_REF] Babadjian | Quasi-static evolution in nonassociative plasticity: the cap model[END_REF][START_REF] Francfort | Small-strain heterogeneous elastoplasticity revisited[END_REF][START_REF] Gilles A Francfort | A case study for uniqueness of elasto-plastic evolutions: The bi-axial test[END_REF], the idea is to discretize in time and find local minimizers of the total energy. In this chapter we focus only on perfect elasto-plasticity materials and set aside the damage. We start with concepts of generalized standard materials in the section 4.1. Then using some convex analysis [START_REF] Ekeland | Convex analysis and variational problems[END_REF][START_REF] Temam | Mathematical problems in plasticity[END_REF] we show the equivalence with the variational formulation presented in the section 4.2. The last part 4.3 presents an algorithm to solve perfect elasto-plasticity materials evolution problems. A numerical verification example is detailed at the end of the chapter.

Ingredients for generalized standard plasticity models

For the moment we set aside the evolution problem and we focus on main ingredients to construct standard elasto-plasticity models [START_REF] Germain | Continuum thermodynamics[END_REF][START_REF] Quoc | Stability and nonlinear solid mechanics[END_REF][START_REF] Suquet | Sur les équations de la plasticité: existence et régularité des solutions[END_REF]. This theory requires a choice of internal variables, a recoverable and a dissipation potentials energies where both functionals are convex. The driving forces (conjugate variables) usually the stress and the thermodynamical force lie respectively in the elastic and dissipation potential energies. For smooth evolutions of the internal variables, the material response is dictated by the normality rule of the dissipation potential convex set (flow law rule). By doing so, it is equivalent to find global minimizers of the total energy sum of the elastic and dissipation potential energies.

Consider that our material has a perfect elasto-plastic response and can be modeled by the generalized standard materials theory, which is based on two statements. In all of these experiments main observations of fractures nucleation reported are: (i) formations of shear bands in "X" shape intensified by necking effects, (ii) growing voids and coalescence, (iii) macro-crack nucleation at the center of the specimen, (iv ) propagation of the macro crack, straightly along the cross section or following shear bands depending on the experiment and (v ) failure of the sample when the fracture reaches external free surfaces stepping behind shear bands path. Observed fracture shapes are mostly cup-cones or shear dominating.

The aforementioned ductile features examples will be investigated through this chapter by considering similar geometries such as, rectangular samples, round notched specimens in plane strain condition and round bars.

Pioneers to model ductile fractures are Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [START_REF] Barenblatt | The mathematical theory of equilibrium of cracks in brittle fracture[END_REF] with their contributions on cohesive fractures following Griffith's idea. Later on, a modern branch focused on micro voids nucleations and convalescence as the driven mechanism of ductile fracture. Introduced by Gurson [START_REF] A L Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I -Yield Criteria and Flow Rules for Porous Ductile Media[END_REF] a yield surface criterion evolves with the micro-void porosity density. Then, came different improved and modified versions of this criterion, Gurson-Tvergaard-Needleman (GTN) [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF][START_REF] Needleman | An analysis of ductile rupture in notched bars[END_REF], Rousselier [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF], Leblond [START_REF] Leblond | An improved gurson-type model for hardenable ductile metals[END_REF] to be none exhaustive. The idea to couple phase-field models to brittle fracture with plasticity to recover cohesive fractures is not new and have been developed theoretically and numerically in [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Conti | Phase field approximation of cohesive fracture models[END_REF][START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Crismale | Viscous approximation of quasistatic evolutions for a coupled elastoplastic-damage model[END_REF][START_REF] Miehe | Phase field modeling of ductile fracture at finite strains. a variational gradient-extended plasticity-damage theory[END_REF][START_REF] Wadier | Mécanique de la rupture fragile en présence de plasticité : modélisation de la fissure par une entaille[END_REF][START_REF] Miehe | Phase field modeling of ductile fracture at finite strains. a variational gradient-extended plasticity-damage theory[END_REF].

Gradient damage models coupled with perfect plasticity

Our model is settled on the basis of perfect plasticity and gradient damage models which has proved to be efficient to predict cracks initiation and propagation in brittle materials. Both mature models have been developed separately and are expressed in the variational formulation in the spirit of [START_REF] Mielke | Evolution of rate-independent systems[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Piero | Variational Analysis and Aerospace Engineering, volume 33 of Optimization and Its Applications[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Dal Maso | Quasistatic evolution problems for linearly elastic-perfectly plastic materials[END_REF] which provides a fundamental Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage

(5.11)

Passing E t i (u i , α i , p i ) to left and dividing by h and letting h → 0 at the limit, we obtain,

By integrating by part the integral term in e(v) over Ω \ (J(ζ i ) ∪ J(v)), we get,

(5.13)

where σ i = a(α i )A e(u i )p i . Without plasticity there is no cohesive effect, hence, σ i ν = 0 and the non-interpenetration condition leads to u i • ν ≥ 0 on J(ζ i ), however for a general cohesive model we do not have information for σ i ν on J(ζ i ). So, to overcome this issue we restrict our study to material with tr (p i ) = 0, consequently on the jump set J(ζ i ) we have tr (

The material can only shear along J(ζ i ) which is commonly accepted for Von Mises and Tresca plasticity criterion. Thus, we have v • ν = 0 on J(ζ i ) and naturally σ i ν = 0 on J(v). The last term of (5.13) stands for,

Combining the above equation, (5.12) and (5.13), considering J(v) = ∅ and by a standard localization argument i.e. taking v concentrated around H n-1 and zero 5.1. Phase-field models to fractures from brittle to ductile almost everywhere, we obtain that all the following integrals must vanish,

which leads to the equilibrium and the prescribed boundary conditions,

Note that the normal stress σ i ν is continuous across J(ζ i ) but the tangential component might be discontinuous.

2. Plastic yield criteria on the jump set: Since the above equation (5.15) holds, for J(v) = ∅ in (5.12) we have,

Thus, on each point of the jump set

The right hand side of the above inequality,

Considering Von Mises criterion we get on the left hand side,

Taking the maximum for all ν = 1, and letting σ i = a(α i )ς i we obtain that (5.17) becomes,

This condition is automatically satisfied for Von Mises since a(α i )/b(α i ) ≤ 1. We refer the reader to [START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Francfort | The elastoplastic exquisite corpse: A suquet legacy[END_REF][START_REF] Gilles A Francfort | A case study for uniqueness of elasto-plastic evolutions: The bi-axial test[END_REF] for more details.

Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage 3. Damage yield criteria in the bulk: Taking v = 0 and q = 0 thus J(v) = ∅ in the optimality condition (5.9), such that E t i (u i , α i + hβ, p i ) ≥ E t i (u i , α i , p i ), then dividing by h and passing to the limit, we get, after integrating by parts the ∇α•∇β term over Ω \ J(ζ i ),

The above equation holds for any β ≥ 0, hence, all contributions must be positive, such that in Ω \ J(ζ i ), we have,

The damage yield criterion is composed of the classical part from gradient damage models and a coupling part in b (α). When the material remains undamaged and plasticity occurs, the cumulation of dissipated plastic energy combined with the property that b (α) < 0 leads to decrease the left hand side which becomes an equality up to a critical plastic dissipation. At this moment the damage is triggered.

4. Damage yield criteria in the jump set: From (5.18) we have,

The gradient damage is discontinuous across the jump set J(ζ i ) due to plastic strain concentration and vice versa.

5. Damage boundary condition: From (5.18) we have,

6. Plastic yield criteria in the bulk: Take v = 0 and β = 0 thus J(v) = ∅ in the optimality condition (5.9) such that

where

Since ψ is differentiable by letting h → 0 and applying the subgradient definition to (5.22), we get -∂ψ/∂p i ∈ b(α i )∂H(p ip i-1 ).

We recover the stress admissible constraint provided by the plastic yield surface.

The damage state decreases the plastic yield surface leading to a stress softening property.

7. Flow rule in the bulk: Applying the convex conjugate (Legendre-Fenchel) to the above equation we get,

which is the flow rule in a discrete settings, by letting max(t it i-1 ) → 0 we get the time continuous one.

Damage consistency:

The damage consistency is recovered using the energy balance condition which is not fully exposed here. However the conditions obtained are:

Damage irreversibility in the domain:

The damage irreversibility constraint is,

All of this conditions are governing laws of the problem. The evolution of the yields surfaces are given by the equations (5.19) and (5.23).

Chapter 5. Variational phase-field models of ductile fracture by coupling plasticity with damage

Application to a 1d setting

The goal of this section is to apply the gradient damage model coupled with perfect plasticity in 1d setting by considering a bar in traction. Relevant results are obtained through this example such as, the evolutions of the two yields functions, the damage localization process and the role of the gradient damage jump term which governs the displacement jump set. We refer the reader to Alessi-Marigo [START_REF] Alessi | Variational approach to fracture mechanics with plasticity[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity: variational formulation and main properties[END_REF][START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF][START_REF] Alessi | Coupling damage and plasticity for a phase-field regularisation of brittle, cohesive and ductile fracture: One-dimensional examples[END_REF] for a complete exposition of this 1d application.

In the sequel, we consider a one-dimensional evolution problem of an homogeneous elasto-plastic-damageable bar Ω = [-L, L] stretched by a time controlled displacements at boundaries where damage remains equal to zero. Assume that a unique displacement jump may occur on the bar located at the coordinate x 0 , thus the admissible displacement, damage and plastic strain sets are respectively,

The state variables of the sound material is at the initial condition (u 0 , α 0 , p 0 ) = (0, 0, 0). In one dimensional setting the plastic yield criteria is |τ | ≤ σ p , thus the plastic potential power is given by,

By integrating over the process, the dissipated plastic energy density is σ p p where the cumulated plastic strain is p = t 0 | ṗs |ds. Since no external force is applied, the total energy of the bar is given by,

where E is the Young's modulus and (•) = ∂(•)/∂x. The quadruple state variables (u t , α t , p t , pt ) ∈ C t ×D ×Q×M(Ω, R) is solution of the evolution problem, if the following conditions holds:

1. The equilibrium,

The stress is constant along the bar hence it is only function of time.

Titre : Modèles variationnels à champ de phase pour la rupture de type fragile et ductile: nucléation et propagation

Mots clefs : Modèles à champ de phase pour la rupture, nucléation de fissure, effet d'échelle dans les matériaux fragiles, modèles d'endommagement à gradient, fracturation hydraulique, stabilité des fissures, modèles de plasticités, approche variationnelle, rupture ductile.

Résumé :

Les simulations numériques des fissures de type fragile par les modèles d'endommagement à gradient deviennent maintenant très répandues. Les résultats théoriques et numériques montrent que dans le cadre de l'existence d'une pré-fissure la propagation suit le critère de Griffith. Alors que pour le problème à une dimension la nucléation de la fissure se fait à la contrainte critique, cette dernière propriété dimensionne le paramètre de longueur interne.

Dans ce travail, on s'attarde sur le phénomène de nucléation de fissures pour les géométries communément rencontrées et qui ne présentent pas de solutions analytiques. On montre que pour une entaille en U-et Vl'initiation de la fissure varie continument entre la solution prédite par la contrainte critique et celle par la ténacité du matériau. Une série de vérifications et de validations sur différents matériaux est réalisée pour les deux géométries considérées. On s'intéresse ensuite à un défaut elliptique dans un domaine infini ou très élancé pour illustrer la capacité du modèle à prendre en compte les effets d'échelles des matériaux et des structures.

Dans un deuxième temps, ce modèle est étendu à la fracturation hydraulique. Une première phase de vérification du modèle est effectuée en stimulant une pré-fissure seule par l'injection d'une quantité donnée de fluide. Ensuite on étudie la simulation d'un réseau parallèle de fissures. Les résultats obtenus montrent qu'une seule fissure est activée dans ce réseau et que ce type de configuration vérifie le principe de moindre énergie. Le dernier exemple se concentre sur la stabilité des fissures dans le cadre d'une expérience d'éclatement à pression imposée pour l'industrie pétrolière. Cette expérience d'éclatement de la roche est réalisée en laboratoire afin de simuler les conditions de confinement retrouvées lors des forages.

La dernière partie de ce travail se concentre sur la rupture ductile en couplant le modèle à champ de phase avec les modèles de plasticité parfaite. Grâce à la structure variationnelle du problème on décrit l'implémentation numérique retenue pour le calcul parallèle. Les simulations réalisées montrent que pour une géométrie légèrement entaillée la phénoménologie des fissures ductiles comme par exemple la nucléation et la propagation sont en concordances avec ceux reportées dans la littérature.

Title : Variational phase-field models from brittle to ductile fracture: nucleation and propagation Keywords: Phase-field models of fracture, crack nucleation, size effects in brittle materials, validation & verification, gradient damage models, hydraulic fracturing, crack stability, plasticity model, variational approach, ductile fracture Abstract : Phase-field models, sometimes referred to as gradient damage, are widely used methods for the numerical simulation of crack propagation in brittle materials. Theoretical results and numerical evidences show that they can predict the propagation of a pre-existing crack according to Griffith's criterion. For a one-dimensional problem, it has been shown that they can predict nucleation upon a critical stress, provided that the regularization parameter is identified with the material's internal characteristic length.

In this work, we draw on numerical simulations to study crack nucleation in commonly encountered geometries for which closed-form solutions are not available. We use U-and V-notches to show that the nucleation load varies smoothly from the one predicted by a strength criterion to the one of a toughness criterion when the strength of the stress concentration or singularity varies. We present validation and verification of numerical simulations for both types of geometries. We consider the problem of an elliptic cavity in an infinite or elongated domain to show that variational phase field models properly account for structural and material size effects.

In a second movement, this model is extended to hydraulic fracturing. We present a validation of the model by simulating a single fracture in a large domain subject to a control amount of fluid. Then we study an infinite network of pressurized parallel cracks. Results show that the stimulation of a single fracture is the best energy minimizer compared to multi-fracking case. The last example focuses on fracturing stability regimes using linear elastic fracture mechanics for pressure driven fractures in an experimental geometry used in petroleum industry which replicates a situation encountered downhole with a borehole called burst experiment.

The last part of this work focuses on ductile fracture by coupling phase-field models with perfect plasticity. Based on the variational structure of the problem we give a numerical implementation of the coupled model for parallel computing. Simulation results of a mild notch specimens are in agreement with the phenomenology of ductile fracture such that nucleation and propagation commonly reported in the literature.