Tout d'abord, je souhaite remercier mon directeur de thèse, le professeur Julien Yvonnet pour son encadrement et son investissement durant ma thèse, j'ai apprécié sa pédagogie et son humilité.

Il a pu me montrer, au-delà de mon travail de thèse, l'envers du décor de la vie académique au sein du laboratoire Modélisation et Simulation Multi Echelle (MSME) de l'université Gustave Eiel.

De la même manière, il a pu me faire découvrir le processus d'organisation et de préparation de conférences scientiques nationales et internationales. Ce fut très captivant.

Je remercie mon second directeur, le professeur Johann Guilleminot de

Background

This thesis is part of a research eort developed within, and funded by, the LabEx MMCD (which is a French acronym for Multi-scale Modeling and Experimentation of Materials for Sustainable Construction) at Université Paris-Est, France. The objective of this LabEx is to advance and disseminate fundamental knowledge for the understanding of the physical phenomena that govern the properties of construction materials and geomaterials. In their natural environment, these materials are in permanent contact with humidity causing swelling and hydric shrinkage, inducing many problems in civil engineering, for example during severe drought. In this context, the study of crack propagation in clay materials with humidity is of particular interest. With the goal of better understanding the degradation mecanisms on civil engineering structures, and specically here the failure and leakage in underground radioactive storage facilitiessuch as the structures located in Bure, France (see Fig. 1.1). The design and construction of such structures involve various types of safety specications, including certied structural integrity with respect to cracking mechanisms driven by aging and humidity. In order to assess the quality of containment and its evolution over time, extensive experimental studies were therefore carried out to characterize the nature and behavior of the geological layers (with Callovo-Oxforian argillaceous rocks being collected on site) [Wang, 2012].

Figure 1.1 Nuclear waste storage structure in Bure, France (in French) [Vanlerberghe, 2017].

From a simulation standpoint, a deep understanding of the underlying physical phenomena is required to enhance the accuracy and the predictiveness of the computational fracture models.

Local mechanisms, in particular, depend on the nature of the constitutive materials (which can be, e.g., cementitious materials and rocks in civil engineering applications), and investigations at ne scales are often crucial to develop predictive material failure models. In this case, the ne scale randomness that is typically exhibited by geological materials can strongly aect simulation outcomes and must thus be integrated in the numerical setting.

Following the above discussion, the overarching goal of this thesis is to develop a modeling framework for micro-cracking in clay materials that combines computational modeling, enhanced by multiscale and stochastic aspects, with experimental characterization.

In the rest of this chapter, we rst review the state of the art regarding numerical approaches for modeling crack propagation.

Research Objectives and Thesis Outline

The thesis is organized as follows. After an extensive review of the cracking models commonly used in the literature, Chapter 2 focuses on the construction of a probabilistic model in the context of random heterogeneous material cracking problems. Then chapters 2, 3 and 4 focus on the aspects of hydric shrinkage of clays with respectively a study on modeling and numerical implementation, followed by an experimental study and a simulation section. The last part Chapter 2 presents the construction of a probabilistic model in the case of a heterogeneous environment. The formulation used here aims at identifying the material parameters of a multi-scale cracking problem using the formulation of the phase eld method. The identication of the problem is based on Monte-Carlo simulations. The mesocopic material properties are computed on the microstructure through the sliding window approach, and the cracking properties are obtained by an inverse problem based on the mean maximum response between the micro and meso scales.

Chapter 3 raises the issue of hydric shrinkage. This part tackles the modelling aspects with the extension in the general framework to large deformations as well as the coupling with a hydric shrinkage model. The problem formulation, the numerical implementation and the benchmarks are described.

Chapter 4 exposes the experimental part of the thesis, where wet clay and rigid inclusion samples are followed during the drying process. The setup as well as the experimental conditions are presented. The objective here is to identify and quantify the mechanisms during the drying process. For this purpose, recent image analysis techniques have been used in order to access to the deformation map, providing valuable and precise information on the complex physics of this phenomenon.

Chapter 5 aims to reproduce, through numerical simulation, the dynamics observed experimentally. The observation is focused on the crack propagation pattern, as well as on the local deformations of the samples. The comparison of these simulations with the experiments will allow us to analyse the relevance of the model and to identify what has not been taken into account. On the basis of the available information, new strategies will be proposed for future improvements. Molecular Dynamics is a well-known simulation tool (see e.g. [Rapaport, 2004]) to simulate complex phenomena at the atomistic scale. This technique is restricted for scales where all atoms are explicitly described. Non-reversible phenomena, dislocations and other microscopic mechanisms leading to crack activation are easily described, as well as post-fracture frictions, by simply considering the interactions between atoms as nonlinear functions in a Newtonian Dynamics approach (discrete mechanics). The technique consists in solving the dynamic equilibrium (second Newton's law) of atom positions based on a Hamiltonian formulation:

Computational Modeling of Crack Propagation

M d 2 q(t) dt 2 = -∇V (q(t)) + f ext , (1.1)
where M is the mass of each particle and V is a potential (see Fig. 1.2) whose general model is given in the form:

V (q) = v (2) ij (q ij ) + v (3) ijk (q ij , q ik ) + v (p)
ijk(...) (q ij , q ik , q i (...)) .

(1.2)

Figure 1.2 MD particles scheme -MD/FEM coupled.

The empirical potentials are dened as functions of the number of interacting particles p and the nature of the atoms involved. For example, p = 2 typically corresponds to axial forces of repulsion or attraction (for, e.g., steric repulsion, Coulomb or dipolar charge attraction/repulsion), while p = 3 and p ≥ 3 can be associated with, e.g., torsional forces and other types of contributions (e.g., metallic bounds), respectively. The MD does not require complex models and can be easily implemented. For the specic case of cracking, the crack initiation and propagation is naturally taken into account without specic treatment: when atoms are far from each other's, using appropriate potentials the interaction forces decrease and atoms do not have any more interactions (see Fig. 1.3). [Souguir, 2018] works (with LabEX MMCD) on adapted potentials so-called reactive potentials that are relevant over large ranges of interatomic distances. However, the major drawback lies in the drastic restriction of space and time domains for which the simulations are conducted, inducing intractable computational times and memory requirements for engineering scale problems. For example, a volume V = 1 µm 3 contains roughly 10 11 atoms for a silica crystal (10 11 roughly corresponds to 6.5 Tera-bytes of memory). Another strong limitation lies in the stable time steps related to such small scales, which are of the order of femtoseconds (∆t ∝ 10 -15 s). To decrease these limitations, bridging techniques between atomistic and continuum approaches have been developed, where the crack phenomenon is localized in a small discrete (atomistic) area, while the rest of the domain (uncracked) is modeled with continuum approaches. In this context, we can mention the Bridging Domain approach [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF] or the Quasi-continuum approach [START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF] (see Fig. 1.2).

In the literature several works concerning cracking at the atomistic scale focus on the phenomena of pre-crack instability [START_REF] Abraham | Instability dynamics of fracture: a computer simulation investigation[END_REF], Souguir, 2018] and the propagation velocity [START_REF] Zhou | Dynamic crack processes via molecular dynamics[END_REF] stress.

Figure 1.3 MD snapshots and zoom of brittle fracture for one critical tensile strain step load ε [START_REF] Wang | High damage tolerance of electrochemically lithiated silicon[END_REF]. 1.3.1.2 Discrete Element Method (DEM) The Discrete Element method was introduced by [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] in rocks mechanics and granular media. Such approach is similar to molecular dynamics, but particles are here considered at a larger scale: they can constitute physical discrete particles such as sand grains, or represent a portion of matter. These particles are modeled by rigid elements interacting through contact and adhesion forces F (i) c , where i is the index of the particle, varying from 1 to N, with N is the number of interacting neighbors (friction, plasticity, adhesion, cohesive contact) having all their own model. The solution to this problem consists in adding to each particle F c all these contributions and in using an integration method to calculate the next position δu of each particle of the system following the MD scheme (see Fig. 1.4). The advantage of this method is to be able to faithfully represent the crack phenomena at the particle scale (see Fig. 1.5). However, when considering continuous media, the denition of particles and interaction models becomes delicate. Indeed, the same scaling rules as in MD prevents to apply the technique to macroscale problems, as a realistic structure might encompass an unaordable number of particles. As in MD, the method can be coupled with continuum approaches, as e.g. in the "Combined Finite-Discrete Element" (CFDE) method widely used to deal with multi scale problems more eciently. Such technique has been applied to ceramics, rock, powder, impact and cracking [START_REF] Munjiza | Combined single and smeared crack model in combined nite-discrete element analysis[END_REF].

Figure 1.5 Crack path obtained by DEM in vitreous biopolymer material [START_REF] Hedjazi | Application of the discrete element method to crack propagation and crack branching in a vitreous dense biopolymer material[END_REF]. 1.3.1.3 Lattice Spring Model (LSM) The Lattice Spring model (LSM) is another discrete approach, aiming at modeling continuum matter with discrete approaches, to facilitate the description of crack propagation. This method is inspired by condensed matter physics and is frequently used to simulate deformation and fracture. It has been shown that in the assumption of linear elasticity, LSM can be equivalent to the Finite Element Method (see [START_REF] Ostoja-Starzewski | Spring network models in elasticity and fracture of composites and polycrystals[END_REF]). In this approach the nodes are associated with non-volumetric particles connected by springs. The deformation of the spring is calculated between two particles position (x i , x j ) as ε (ij) = E(x i , x j ) , (1.3) where E is a deformation function. Under mechanical stress, each node moves and change of position induces a nodal force such as:

F i = v=1,...,vmax K(x i , x v ) , (1.4)
where K is a function associating the spring stiness between the i-th node of interest, the v-th neighbor node and theirs positions. In this method, the static equilibrium is invoked

F i + F ext = 0 (1.5)
with F ext is the imposed external force. In the application, when a spring exceeds a critical threshold ε > ε c , the crack starts to initiate. The broken springs no longer transmit forces to their neighbors, the system enters an unsteady state, the new rearrangement is done by solving, in each node linked to the broken springs (see Fig. 1.7), the following equation:

F i = v=1,...,v max-k K(x i , x j ) = 0 .
(1.6) Figure 1.6 Lattice crack process.

Above, k is the number of broken springs for a node. This can cause a variation of deformation which leads to other ruptures. This chain process, which is solved by the equilibrium of local forces at each node, ends when all induced deformations return below the critical threshold.

Numerically this method imposes a structured mesh, which can drastically limit the approach when considering structures with complex geometries.

In the literature [START_REF] Vogel | Studies of crack dynamics in clay soil: Ii. a physically based model for crack formation[END_REF], Malthe-Sørenssen et al., 1998] have shown the eciency of these models by comparing its numerical modeling with experiments on homogeneous clay materials under shrinkage conditions. The model shows good crack pattern reproduction, bifurcation angles and model parameters having a physical sense (randomness of heterogeneous friction) on the nal state of the sample. However, such approach induces ad-hoc models of interactions and cannot be easily extended to more complex phenomena (plasticity, heterogeneous materials, complex geometries). 1.3.1.4 Peridynamics This technique originally proposed by [Silling, 2000] is closely related to MD or DEM. It assumes that the material is described by discrete particles which represent the continuous matter. The dynamic equilibrium is solved according to: ρ(x)ü(x, t) = R f (u(x , t) -u(x, t), x -x, x)dV x + b(x, t) ,

(1.7) Figure 1.7 Desiccation crack modelled by LSM [START_REF] Vogel | Studies of crack dynamics in clay soil: Ii. a physically based model for crack formation[END_REF].

where ρ is the density, u the displacement, f a density force vector function that integrates the displacements and position of its neighbors x . The right-hand term of this equation is a convolution term, similar to nonlocal elasticity models. Once discretized, it matches the Discrete Element Method above, except that the number of interaction is larger, as one point is related to a large number of surrounding points. Again, in this method, the rupture criterion is based on a threshold elongation (u(x, t) -u(x , t)) between two particles. When this link breaks up there is no more interaction and the forces are redistributed on the neighborhood links, as in LSM. It then inherits the benets of DEM and MD to easily model complex cracks patterns (see Fig. 1.8) but has the same drawbacks (ambiguous denition and calibration of interaction models to reproduce general mechanical behavior, spatial convergence issues, etc.) Figure 1.8 Peridynamic simulation of impact [START_REF] Littlewood | Strong local-nonlocal coupling for integrated fracture modeling[END_REF].

Cohesize zone models 1.3.2.1 Cohesive elements

Cohesive models have been introduced in the work of [Dugdale, 1960] and [Barenblatt, 1961] for ductile and brittle materials. In these models, cracks are described by cohesive laws on the lips of the crack (called cohesive zones see Fig. 1.9). The energy of the system with cohesive element can be written as follows :

E = Ω Ψ(ε)dΩ + Γ Ψ I ([[u]])dΓ , (1.8)
where [[δu]] = u + -u -denotes the displacement jump across Γ. Variation of the energy with respect to the displacement eld gives the weak form:

Ω σ(ε) : ε(δu)dΩ + Γ f ([[u]]) : [[δu]]dS = ∂Ω F F ext • δudS .
(1.9)

In this formulation, f ([[u]]) = ∂Ψ I /∂ [[u]] is interpreted as the cohesive force (depending on a chosen cohesive law) between the lips of Γ depending of the value of the displacement jump [[u]],

Chapter 1. Introduction dening the local eect of the interface on the global response (see Fig. 1.9). Numerically, the cohesive model is applied between the boundary of elements whose nodes are doubled. It has been shown that this technique has shortcomings, including (i) non-convergence of the energy with mesh renement; (b) strong mesh-dependency (the cracks are constrained to follow the boundary of the elements) and (iii) additional nonphysical compliance in the material due to the surface spring layer model added into the energy. The reader is referred to [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF] crack zone scheme.

0 δ lin δ CZM ||σ|| F σ c
Linear law exponential law basic material traction force response under dierent cohesive laws.

Figure 1.9 Cohesive zone model approach.

for an extensive review on commonly used cohesive zone models. We note that [START_REF] Vo | Numerical modelling of desiccation cracking of clayey soil using a cohesive fracture method[END_REF] used this approach to investigate the desiccation cracking process of clay soils with experimental results.

Augmented Finite Element Method (A-FEM)

This method, introduced by [START_REF] Hansbo | A nite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF], aims at introducing discontinuity in the model at the level of the numerical discretization. In other words, it does not constitute a crack propagation model but a numerical technique to introduce discontinuities in a convenient manner in regular meshes, even though it can be combined with, e.g., a cohesive model or a crack propagation criterion. The main idea is summarized below. In this technique, a discontinuous enrichment is induced in the FEM displacement approximation by separating an element in several sub-entities called mathematical element (ME) dened on Ω α as illustrated on Fig. 1.10 α = {1, 2}. On each ME, the displacement u(x) eld is interpolated in the FEM sense, using classical nodes u i and "ghost" nodes u i that correspond to split nodes:

∀x ∈ Ω i , u α (x) = N α (x) • {u i,α } , (1.10) with u i = (u i , u i ) T .
The forms of elementary matrices are provided as

K α • {u i,α } = f Γα + f α , (1.11) where K α = Ωα B T α CB α dΩ; f Γα = Γα N T α • t α dS; f α = ∂Ωα N T α • F α dS .
(1.12)

Figure 1.10 AFEM inter-element superposition.

Above, C denotes the elasticity matrix representation, F α the surfaces forces, and t α is the internal force related to the discontinuity where the cohesive law can be used (t α = t α (δ)) [START_REF] Mergheim | A nite element method for the computational modelling of cohesive cracks[END_REF]. All the interpolation and its derivatives matrix are modied from FEM and explicited in [START_REF] Ling | An augmented nite element method for modeling arbitrary discontinuities in composite materials[END_REF]. Numerically this method does introduce neither additional degrees of freedom nor remeshing. However, severe issues come with mesh dependency, illconditioned matrix systems, and non-convergence with respect to the mesh size. Finally, an extension to the method in 3D is delicate.

Methods based on linear fracture mechanics

The mechanics of linear fracture is based on a macroscopic and energetic description of the cracking process described by the pioneers' work Grith & Irvin. In its classical theoretical framework the solutions (obtained under restrictive hypotheses) are limited with the increase of the complexity of the problem (geometry, heterogeneities...), involving an increase in mathematical complexity when computed analytically. Combining such approaches with numerical methods like FEM allows dening accurate crack propagation criteria for more complex congurations.

FEM remeshing techniques

The classical use of FEM in the context of crack propagation is based on: (i) estimation of the stress state at the crack tip; (ii) dening a propagation condition and (iii) remeshing to update the nite element mesh used to again approximate the stress eld.

The estimation of the mechanical state of the material around the crack tip within the framework of the linear fracture mechanics can be quantied either by the stress level or energetically.

For the stress approach, the critical criterion is based on an asymptotic solution [Irwin, 1957] at the crack tip under the assumption of a material with linear elastic behavior. For example the stress solution at the crack tip is given in 2D under plane strain conditions by:

               σ rr K I 4 √ 2πr [5 cos(θ/2) -cos(3θ/2)] + K II 4 √ 2πr [-5 sin(θ/2) + 3 sin(3θ/2)] , σ θθ K I 4 √ 2πr [5 cos(θ/2) + cos(3θ/2)] + K II 4 √ 2πr [-3 sin(θ/2) -3 sin(3θ/2)] , σ rθ K I 4 √ 2πr [sin(θ/2) + sin(3θ/2)] + K II 4 √ 2πr [cos(θ/2) + 3 cos(3θ/2)] ,
(1.13) where (r, θ) is the local coordinate basis dened on the crack tip as illustrated in Fig. 1.11 and (K I , K II ) are stress intensity factors.

Figure 1.11 Local crack basis and J-integral domain.

These quantities can be evaluated numerically based on the FEM results in the mesh elements.

However, the singularity at the crack tip requires a very ne mesh and a volume average of stress intensity factors may result in a poor accuracy. The G-θ method consists in estimating the rate of energy restitution in the form of a contour path integral (J-integral [Rice, 1968] see Fig. 1.11) closed by the lip of the crack, thus avoiding the numerical issues at the crack tip according to

J = J ∂u ∂x • (σ • n) - 1 2 (σ : ε)n x dS , (1.14)
where the normal n or partial derivative operators have to be constructed and J calculated. The direction (angle θ c ) of propagation is found through verifying a given criterion, such as maximum radial stress of minimum energy. After each crack propagation, the mesh must be reconstructed from the crack and boundary, which might be costly and non-robust, especially for 3D geometries or changes of topology of the crack.

eXtended Finite Flement Method (X-FEM)

The eXtended Finite Element Method aims at avoiding the remeshing step in the above FEM crack propagation process within Finite Elements and linear fracture mechanics. The main idea is to enrich the FEM discretization to introduce both discontinuities and singularities related to the crack propagation directly within the FEM approximation through additional degrees of freedom at the nodes, which are associated to specic discontinuous and singular shape functions.

The method is actually a special case of the Partition of Unity Method [START_REF] Babu²ka | The partition of unity method[END_REF], which consists in improving a solution of a FEM problem by injecting, into the classical approximation, an additional function involving a basis that is constructed based on an expected solution (enrichment):

u(x) = i N i (x)u i + i,p N i (x)Φ i,p (x)a i,p , (1.15)
where a i,p denote additional unknowns associated to the enriched approximation. [START_REF] Belytschko | Elastic crack growth in nite elements with minimal remeshing[END_REF] applied this method to crack propagation by dening Φ i,p as basis functions 1.3. Computational Modeling of Crack Propagation describing the asymptotic solution in the vicinity of the crack:

u(x) = r→0                  1 2µ r 2π (K I cos( θ 2 ))(K -cos θ) + K II sin( θ 2 )(K + cos θ + 2) , 1 2µ r 2π (K I sin( θ 2 ))(K -cos θ) -K II cos( θ 2 )(K + cos θ -2) , 2 µ r 2π (K III sin( θ 2 )) ,                  (1.16)
where K is the Kolosov constant is dened for isotropic, plane stress (as K := 3 -4ν) and plane strain (as K := (3 -ν)/(1 + ν)), (r, θ) denes the geometry of the crack by a local coordinate system based on the crack tip, and (K I , K II , K III ) are the stress intensity factors dened for the basic fracture modes from linear elasticity [Irwin, 1957]. [START_REF] Moës | A nite element method for crack growth without remeshing[END_REF] enriched the basis with another additional discontinuous function H separating the crack lips and alleviating the remeshing step within the FEM procedure:

u(x) = u(x) = i N i (x)u i + i,p N i (x) Φi,p (x)a i , p + i H(x i )a i ,
(1.17)

with the discontinuous function H being dened by:

H(x) = + 1, if φ > 0 , -1, if φ < 0 , (1.18)
and φ is a local coordinate attached to the crack surface dened by a couple of level set functions (φ, ψ) in 3D [START_REF] Stolarska | Modelling crack growth by level sets in the extended nite element method[END_REF] (see Fig. The introduction of a damage model within a micromechanics framework was rst conducted in [Kachanov, 1958]. The model characterizes damage induced by the formation of micro-cracks.

Chapter 1. Introduction

In this context, an isotropic degradation function was introduced and involves a scalar damage parameter d, ranging from 0 to 1. The constitutive equation relating the stress eld σ and the strain eld ε of an isotropic damage model is written as:

σ = (1 -d)C : ε , (1.19)
where C is the stiness tensor of the elastic, healthy material. In the case of anisotropic damage eects, a damage tensor D must be introduced in lieu of d. In addition, an evolution law is necessary for d. This damage law may be chosen in order to reect the behavior of the considered material. For example, for quasi-brittle materials, the following model has been adopted in [START_REF] Peerlings | Gradient-enhanced damage modelling of concrete fracture[END_REF]]:

d = 0 if κ < κ 0 , 1 -κ 0 κ [(1 -α) + α exp -β(κ-κ 0 ) ] .
(1.20)

In (1.20), the scalar parameter β describes the softening behavior, α is a scalar which controls the residual state in the post peak stage, κ 0 is the threshold for the initiation of damage and κ is a history scalar parameter which takes the largest value of an equivalent strain ε function of ε (see below). Damage evolution is governed by the Kuhn-Tucker inequalities as follows:

κ ≥ 0 , f (ε, k) ≤ 0 , ḋf (ε, k) = 0 , (1.21)
where f (ε, k) = ε-κ is the loading function driving the evolution of damage. Early developments in the context of numerical methods can be found in [Krajcinovic, 1983,Chaboche, 1988,Lemaitre and Chaboche, 1994]. Various denitions for ε have been later proposed. For example, according to the [Mazars, 1984] criterion, cracks can only propagate due to tensile strains, according to:

ε(ε) = ε i : ε i , (1.22)
where ε i are principle strains and ε i = |ε i |+ε i 2 . For ductile fracture, the modied von Mises equivalent strain is usually dened as:

ε(ε) = k -1 2k(1 -2ν) I 1 (ε) + 1 2k (k -1) 2 (1 -2ν) 2 I 2 1 (ε) + 12k (1 -ν) 2 J 2 (ε) , (1.23)
where k is the tensile/compressive strength ratio which is adapted depending on the material, ν is the Poisson's ratio, I 1 and I 2 are the rst two invariants of the strain tensor. Another choice is the so-called smooth Rankine calibration [START_REF] Jirásek | Numerical aspects of the crack band approach[END_REF]:

ε(ε) = 1 E σ i : σ i , (1.24)
where σ i is the principle stress tensor, and E is the Young's modulus. Such local models induce numerical issues when implemented using the nite element method.

In particular, a non-convergence of the response with respect to mesh density is typically observed (see the left panel in Fig. 1.13), for the strain localizes in an innitely narrow band of elements as the mesh is rened, and crack trajectories (identied as localization bands of damage) turn out to be highly sensitive to mesh construction (regular or unstructured). To avoid these issues, dierent regularization techniques were introduced, including the use of: Figure 1.13 Local and regularized damage models: response of a concrete three-point bending beam [Jirásek, 2004].

A Cosserat continuum or a micropolar model (see, e.g., [START_REF] Lakes | Fracture mechanics of bone with short cracks[END_REF], Borst, 1991]) that involves a higher-order continuum model, and an associated internal length scale regularizing the energy of the system across the mesh.

An articial viscosity technique [START_REF] Etse | Failure analysis of elastoviscoplastic material models[END_REF].

Gradient Enhanced Damage (GED) models [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]].

An integral-type regularization [START_REF] Baºant | Nonlocal continuum damage, localization instability and convergence[END_REF], JiráSek, 2007, Baºant and Oh, 1983] where a convolution operator is used to regularize the strain eld.

Higher-order gradients of deformation [Needleman, 1988, Needleman, 1990, Bourdin et al., 2000, Karma et al., 2001, Miehe et al., 2010a].

An illustration of this regularization of damage model on the mechanical response compared to the classical damage model is depicted in the right panel in Fig. 1.13.

Pros and cons of these methods have been discussed in [START_REF] Borst | Fundamental issues in nite element analysis of localization of deformation[END_REF]. Among these techniques, the last two are the most used in computational analysis and are called regularization techniques. A very popular method in this context is the so-called phase eld method to fracture [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], Karma et al., 2001, Miehe et al., 2010a], which will be detailed in the sequel and used in this thesis.

Phase Field (PF) method

The so-called "phase eld method" presented in the following is a non-local damage method which has been initiated by dierent communities in both physics and mechanics. The name "phase eld" comes from the fact that in the physics community, this approach has been developed following the well-known material phase change models with smeared interfaces (see below)

to describe the transition between an undamaged phase and a damaged phase in the material, mimicking the transition between one material phase to another one. In the mechanics community, the technique has been developed from the concept of variational approach to fracture, as presented below. Both origins of the method are energetic minimization concepts. In what follows, we present the main ideas and ingredients of the method. Recent and complete reviews on phase eld methods can be found e.g. in [START_REF] Ambati | A review on phaseeld models of brittle fracture and a new fast hybrid formulation[END_REF] or [START_REF] Wu | Phase eld modeling of fracture[END_REF].

Chapter 1. Introduction

Variational approach to fracture

The variational approach to fracture has been originally proposed by [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]]. The method has been recast in a regularized form by [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and is nowadays called phase eld method by most authors [Kuhn andMüller, 2008, Miehe et al., 2010a]. While in its original form the method estimated the energy based on unknown sharp discontinuities, the regularized approach uses a continuous (called damage) eld to describe the discontinuities thanks to a Mumford -Shah functional as proposed in [Ambrosio and Tortorelli, 1990], which gratefully simplies the minimization process with respect to both displacements and damage elds. The regularization process involves a parameter , which denes an internal length variable. The obtained models are close to gradient-enhanced damage models [START_REF] Borst | Gradient damage vs phaseeld approaches for fracture: Similarities and dierences[END_REF] The solid may contain cracks denoted collectively as Γ. The total energy of the system is dened, in the absence of body forces, as:

E = Ω Ψ(ε, Γ)dΩ + g c Γ dΓ - ∂Ω F F * • udΓ , (1.25)
where Ψ(ε, Γ) is the elastic strain density function and g c is the critical energy release rate in the sense of Grith. A regularized form is given by [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], Miehe et al., 2010a]:

E = Ω Ψ(ε, d)dΩ + g c Ω γ(d, ∇d)dΩ - ∂Ω F F * • udΓ , (1.26)
where γ denotes the crack density function, g(d) is a degradation function such that g(0) = 1, g(1) = 0 and g (1) = 0 and Ψ is an elastic strain density function (see examples in the following).

Computational Modeling of Crack Propagation

The variational approach to fracture as proposed in Bourdin, Francfort and Marigo [Francfort and Marigo, 1998, Bourdin et al., 2000, Bourdin et al., 2008] and developed in a convenient algorithmic setting by Miehe [Miehe et al., 2010a], and is presented here. The phase eld formulation implies: (a) minimization of the total energy with respect to the displacement eld u and (b) minimization of the energy with respect to the scalar eld d describing the crack surface in a smooth manner. This second minimization is subjected to an inequality constraint: ḋ ≥ 0.

D δu L = 0 , (1.28) D δd L = 0, 0 ≤ d n ≤ d n+1 , (1.29)
where D δv f (u) is the Gateaux (directional) derivative, dened by:

D δv f (u) = f dα (u + α δv) α=0
.

(1.30)

The rst equation (1.28) denes the mechanical problem, whereas the second one (1.29) denes the so-called phase-eld problem.

Link with material phase eld

The phase eld approach has been derived independently within the Physics community, by considering the damage process as a phase transition between a sound phase and a damaged phase. To develop this idea, [START_REF] Aranson | Continuum eld description of crack propagation[END_REF] introduced the following free energy functional:

Ψ = Ω D δξ |∇ξ| 2 + φξ dΩ , (1.31)
where ξ = (1 -d), d ∈ [0, 1] describes the phase in the domain (d = 0) outside of the crack and (d = 1) within the crack, D . denoting the directional derivative and φ being a polynomial function. This formulation leads to a time crack evolution problem, dened as:

ḋ = -D δξ ∆φ + φ(d)[α 1 (1 + (tr(ε) -α 2 )φ -γ u)∇φ] ,
(1.32)

where α i is a model parameter and .

is the temporal derivative operator. This formulation of the free energy was enriched in [START_REF] Karma | Phase-eld model of mode iii dynamic fracture[END_REF], by integrating the balance momentum equation, by changing the stress/strain relationship and crack evolution, and by dening the form of the energy functional as

Ψ(u, d) = Ω [g(d)(ψ 0 (∇(u)) -ψ c ) + V (ξ) + 1 2 D δξ |∇ξ| 2 ]dΩ , (1.33)
where g is a function coupling elasticity and damage. Above, V is a double-well potential, ψ c is a critical strain energy, and ψ 0 the elasticity part subjected to some dierent model as:

   ψ 0 = ε : C : ε , [Karma, 2001] ψ 0 = 1(trε + )ψ 0 + 1(trε -)(ψ 0 - 1 2 αk(tr(ε)) 2 )
, [START_REF] Henry | Dynamic instabilities of fracture under biaxial strain using a phase eld model[END_REF] (1.34)

where k is the bulk modulus, α a scalar parameter (α > 1) guarantying a safe material under compression and 1(trε + ) = 1, if tr( ) > 0 (0 otherwise) respectively 1(trε -) = 1, if tr( ) < 0 (0 otherwise). The two models are subject to mode III fracture [Karma, 2001] and mode I & II fracture [START_REF] Henry | Dynamic instabilities of fracture under biaxial strain using a phase eld model[END_REF].

The phase eld method has several crucial advantages as compared to other numerical methods for crack modeling. Specically:

(i) The initiation of cracks from undamaged structures or materials can be handled.

(ii) Branching and merging of cracks are naturally taken into account.

(iii) Arbitrary geometrical congurations of crack networks can be treated (see, e.g., Fig. 1.15).

(iv) The phase eld problem dened with the rst-order damage gradient can be solved with classical nite elements, without modifying existing codes. One drawback which can be reported is the requirement of a ne mesh along the crack path.

The denition of a suitable crack density function, the treatment of self-contact within the crack and the denition of a proper degradation function constitute key ingredients of the phase eld formulation. These aspects are reviewed below.

Denition of a crack density function

Many choices are possible regarding the denition of the crack density function γ in (1.26).

In [Wu, 2017], a general denition was introduced as:

γ(d, ∇d) = 1 c 0 2 D(d) + 2 |∇d| 2 , (1.35)
where c 0 a scaling parameter and D is a function that characterizes the distribution of the smooth cracks, with c 0 = 4 1 0 (D(x)) 1/2 dx. The function D is expressed as

D(d) = ξd + (1 -ξ)d 2 ,
(1.36)

where ξ ∈ [0, 2] is a scalar parameter, and D(d) ∈ [0, 1] ∀d ∈ [0, 1]. In this context, the prole of the diused damaged normal to the crack is dened in 1D according to

d(x) = Argmin Ω * γ(d, ∇(d))dΩ , (1.37)
where d is dene on its domain Ω * . Several particular cases of this framework can be found in the literature.

The choice D(d) = d 2 and ξ = 0 is widely used in the literature [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] and

leads to γ(d, ∇d) = 1 2 1 d 2 + ∇(d) • ∇(d) .
(1.38)

The corresponding 1D prole is then given by d(x) = exp(-|x|/ ).

The denition D(d) = d and ξ = 1 considered 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] yields

γ(d, ∇d) = 8 3 1 d + |∇d| 2 .
(1.39)

This form induces a linear part of the response before failure.The 1D prole is given by

d(x) = (1 -(|x|/2 )) 2 ∀x ∈ Ω * = [-2 , 2 ] 3 .
The case D(d) = 2d -d 2 and ξ = 2 was proposed in [Wu, 2017] and gives

γ(d, ∇d) = 1 π 1 (2d -d 2 ) + |∇d| 2 .
(1.40)

The 1D prole of the damage is a sinus function d

(x) = 1 -sin(|x|/ ) ∀x ∈ Ω * = [-π /2, π /2] 3 .
An extension involving higher-order gradients was introduced in [START_REF] Borden | A higher-order phase-eld model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF] as

γ(d, ∇d) = 1 2 1 d 2 + 2 ∇(d) • ∇(d) + 3 16 ∆ 2 (d) .
(1.41)

This form leads to a more regular damage prole given in 1D by d(x) = exp(-|x|/2 )(1 + |x|/(2 )). The main drawback of this choice is that C 1 continuity is required for the nite element scheme, which may increase the computational complexity, especially in 3D.

These examples of regularization are illustrated in Fig. 1.16 (left). The reader is referred to [START_REF] Wu | Phase eld modeling of fracture[END_REF] for an extensive list of regularized functions. Note that each model of regularization introduces the regularization parameter . The value of depends on the critical stress value and the degradation function aecting the maximal force response (see the right panel in Fig. 1.16, based on [START_REF] Kuhn | On degradation functions in phase eld fracture models[END_REF]; see also Table 1 .3.1.2). A discussion on the choice of can be found in, e.g., [START_REF] Nguyen | On the choice of parameters in the phase eld method for simulating crack initiation with experimental validation[END_REF].

- Critical stress with respect to . Treatment of self-contact (strain split)

To describe the dierence between compressive and tensile damage and to model self-contact within the crack, the strain density function can be split into two parts as:

Ψ = Ψ + + Ψ -, (1.42)
where Ψ + and Ψ -will be dened momentarily. The idea is then to associate damage with the positive part only, according to:

Ψ = g(d)Ψ + 0 + Ψ -, (1.43)
where Ψ + 0 denotes the positive part of the strain density energy for the undamaged material, and g is a degradation function that depends on the damage parameter d. Several choices were proposed for such a decomposition in the case of isotropic damage: [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] proposed a decomposition based on the spherical and deviatoric parts of the strain tensor (denoted as ε H and ε D respectively). In this model the damage is generated by a positive spherical and deviatoric strain.

Ψ + (ε) = 1 2 k[tr(ε)] 2 + µε D : ε D , Ψ -(ε) = 1 2 k[-tr(ε)] 2 , (1.44)
where k and µ are the bulk and shear moduli, and tr(.) is the trace operator. This choice introduces damage for the positive hydrostatic strain (ε H := tr(ε)/n, n the spatial dimension) only, while the deviatoric part (ε D := εε H ) does not induce damage:

σ = g(d)1 (ε + H ) [kε H + 2µε D ] + 1 (ε - H ) [kε H ] ,
(1.45)

where 1(ε

+ H ) = 1 if ε + H > 0 (0 otherwise) and 1(ε - H ) = 1 if ε + H < 0 ( 0
otherwise) and λ and µ denote the elastic Lamé constants.

[ Miehe et al., 2010a] introduced another model, based on the spectral decomposition of the strain tensor:

Ψ ± (ε) = λ 2 T r (ε) ± 2 + µT r ε ± 2 , (1.46)
where λ and µ are the Lamé constants. In Eq. (1.46), ε + and ε -are such that ε = ε + +εand are dened by

ε ± = m i=1 < κ i > ± ϕ (i) ⊗ ϕ (i) , (1.47) in which {(κ i , ϕ (i) )} m
i=1 are the pairs of associated eigenvalues and eigenvectors of the strain tensor ε, and < • > ± is the operator given by

< z > ± = 1 2 (z ± |z|) , ∀z ∈ R .
(1.48)

Within this scheme, the stress tensor is derived according to:

σ = g(d)1(tr(ε + ))[λtr(ε) + 2µε + ] + 1(tr(ε -))[λtr(ε) + 2µε -] .
(1.49) [START_REF] He | Closed-form coordinate-free decompositions of the two-dimensional strain and stress for modeling tensioncompression dissymmetry[END_REF] new failure model is based on the decomposition of the strain tensor into two complementary parts, which are orthogonal in the sense of an inner product where both fourth-order elastic stiness and compliance tensors act as metric operators. Let C be the fourth-order stiness tensor, the traction (denoted by (.) + ) and compression parts (denoted by (.) -) of the strain energy density can be expressed by

ψ e± (ε) = 1 2 ε ± : C : ε ± .
(1.50)

This model can be applied to arbitrary initial anisotropic elastic behavior, in contrast to the above other models.

The orthogonality condition for the positive/negative parts ε ± can be dened as follows

ε + : C : ε -= C : ε + : ε -= 0 .
(1.51)

The requirement described in Eq. (1.51) can be ensured through a method based on elastic energy preserving transformation. Within this framework, the square root of the elastic stiness tensor is introduced

C 1/2 = i Λ 1/2 i ω i ⊗ ω i and C -1/2 = i Λ -1/2 i ω i ⊗ ω i , (1.52)
where Λ i are the eigenvalues of C, and ω i are second-order eigentensors associated to Λ i . Dening ε± as the positive and negative parts of C 1/2 : ε, we compute ε ± by ε ± = C -1/2 : ε± .

(1.53)

The derivatives of ε± with respect to the transformed strain tensor ε also dene two projection tensors P± ( ε) = ∂ ε [ ε± ( ε)], which can be determined following the approach proposed by Miehe [Miehe, 1998]. This implies the complete formulation for the proposed decomposition scheme as follows ε ± = C -1/2 : P± : C 1/2 : ε .

(1.54)

Compared to the scheme proposed by Miehe et al. [Miehe et al., 2010a], this model is computationally more ecient due to the very simple and analytical expressions of the dierent operators, which do not require numerical evaluation of eigenvalues and eigenvectors of the strain tensor.

In this framework, the Cauchy stress σ is obtained as

σ(ε, d) = g(d) ∂ψ e+ (ε) ∂ε + ∂ψ e-(ε) ∂ε = C(d) : ε .
(1.55)

From (1.54) and by introducing P ± = C -1/2 : P± : C 1/2 , the general form of the elastic tensor accounting for damage is dened by C(d) = g(d) P + : C : P + + P -: C : P -.

(1.56)

An illustrative example for a shear test is presented to demonstrate the impact of the aforementioned modeling choices on both the mechanical response (see Fig. 1.17) and crack propagation (see Fig. 1.18). This three force response tests have been performed until the crack aected the boundary of the domain for the Miehe's model.Concerning the shear tests, the strains (positive and negative) are localized axisymmetrically along the pre-crack axis. As the rst damage model [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]] aected all the strain energy (positive and negative), the crack propagates symmetrically and diagonally. For the [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF]

Denition of the degradation function

Several choices can also be made for the degradation function g. The latter is required to satisfy the following conditions:

g(0) = 1 and g(1) = 0 must hold to ensure that the material is initially undamaged and ultimately, fully damaged, respectively.

[ [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]] [Amor et al., 2009] [Miehe et al., 2010a] Figure 1.18 Damage elds associated with the shear test, obtained with three dierent formulations.

g (1) = 0 must hold to ensure a nite value of the stress at the crack tip.

g must be monotonically decreasing to ensure damage decrease away from the crack.

Many degradation functions satisfying these conditions have been proposed in the literature.

[ [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] have used a quadratic polynomial function borrowed from [Ambrosio and Tortorelli, 1990]. Higher-order (quadratic, cubic, and quartic) polynomial degradation functions were studied in [START_REF] Kuhn | On degradation functions in phase eld fracture models[END_REF], and the impact on the mechanical response is illustrated in Fig. 1.19. Note that plastic softening can be taken into account in the formulation proposed in [START_REF] Borden | A phase-eld formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality eects[END_REF]. Some examples of degradation functions are listed below. [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] g [START_REF] Borden | A phase-eld formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality eects[END_REF] g(d) = (1 -d) 2 (quadratic) [START_REF] Kuhn | On degradation functions in phase eld fracture models[END_REF] g

Degradation function Reference

g(d) = (1 -d) 2 [
(d) = (3 -s)(1 -d) 2 -(2 -s)(1 -d) 3 [
(d) = 3(1 -d) 2 -2(1 -d) 3 (cubic)
[ [START_REF] Kuhn | On degradation functions in phase eld fracture models[END_REF] g

(d) = 4(1 -d) 3 -3(1 -d) 4 (quartic)
[ [START_REF] Kuhn | On degradation functions in phase eld fracture models[END_REF] g(d) = 1-d 1-d+md (quasi-linear, m ≥ 1) [START_REF] Geelen | A phase-eld formulation for dynamic cohesive fracture[END_REF] g(d) =

(1-d) 2 (1-d) 2 +md(1+pd) (quasi-quadratic p ≥ 1) [Lorentz, 2017] [START_REF] Miehe | Thermodynamically consistent phase-eld models of fracture: Variational principles and multi-eld fe implementations[END_REF] shows the impact of the η parameter on the post-rupture mechanical response). Other solutions can be found for that purpose, such as in [Lorentz, 2017] and [START_REF] Geelen | A phase-eld formulation for dynamic cohesive fracture[END_REF].

2011], the development of approaches incorporating multiscale and probabilistic ingredients all

together is a natural path to extend the predictive capabilities in fracture simulations.

In this chapter, presented in [START_REF] Hun | Construction of simplied models for crack propagation in random heterogeneous media[END_REF], Hun et al., 2019a] and adapted from the published paper [START_REF] Hun | Stochastic multiscale modeling of crack propagation in random heterogeneous media[END_REF], we propose a stochastic, multiscale-informed phase-eld approach to model crack propagation in random quasi-brittle materials. The objective is to construct a simplied mesoscopic model which does not require an explicit description of ne scales heterogeneities, while retaining the stochastic features that allow reproducing the crack paths and response of the related structures. In the proposed framework, the parameters involved in the elasticity-phase-eld formulation are specically dened through multiscale analysis with non-separated scales. This particular setting ensures consistency with critical subscale information, and allows for the propagation of stochasticity at the macroscopic level. Similar ideas were pursed in the very recent work [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF], with a few noticeable dierences though.

First, the approach developed in the above reference is concerned with dynamical fracture, solved using an asynchronous space-time discontinuous Galerkin method, and it is focused on fracture strength random elds. A phase-eld approach to brittle fracture modeling is alternatively considered and extended here, in which validation is further assessed on a macroscopic quantity of interest. Second, and while both contributions invoke information theory as a rationale to dene probability measures, stochastic modeling aspects and related methodological issues are addressed more extensively hereinafter. Note also that crack paths are simulated in the sequel by propagating a pre-existing crack, whereas crack nucleation sites are identied, for each sample of the microstructure, as the weakest material points in [START_REF] Acton | Voronoi tessellation based statistical volume element characterization for use in fracture modeling[END_REF].

This chapter is organized as follows. The computational approach enabling the description of crack propagation at the microscopic scale is rst detailed in Section 2.2. The phase eld formulation with mesoscopic descriptors (elasticity and damage) is then introduced, and some results comparing the mechanical response and the crack paths at micro and meso-scale are presented in Section 2.3. Stochastic methodologies to represent and subsequently identify the mesoscopic descriptors are introduced and applied to the aforementioned microstructure. The relevance of the framework is nally assessed by comparing macroscopic predictions based on either the reference microscopic model or the proposed mesoscopic stochastic modeling.

Microscopic Crack Propagation

In this section, we dene the random heterogeneous microstructures studied in this work and describe the Monte Carlo approach used to generate realizations of crack paths at the microscale.

Governing equations for the crack propagation problem are introduced within the phase eld method, in the continuity of Section 1.3.4.2. These realizations will be used, in Section 2.4.1, to identify a mesoscale stochastic model constructed in Section 2.3.2.1.

Generation of heterogeneous random microstructures

In order to illustrate the methodology, a prototypical stationary, isotropic random microstructure made up of a homogeneous matrix and monodisperse spheres is selected hereinafter (N p ∈ {1, 2} is dened as a phase index). A two-dimensional square domain Ω = (]0, L[) 2 is considered, with L = 1 mm, and the radius of the inclusions is set to R = 0.04 × L. Plane strain conditions are assumed. A set of θ obs = 1000 realizations was generated by using the molecular-dynamicstype algorithm (based on event-driven molecular dynamics [START_REF] Donev | Neighbor list collisiondriven molecular dynamics simulation for nonspherical hard particles. i. algorithmic details[END_REF]) used in [START_REF] Skoge | Packing hyperspheres in high-dimensional euclidean spaces[END_REF] for hard-sphere packings (under periodic boundary conditions). Each periodized sample contains N inc = 50 non-overlapping heterogeneities. Four independent realizations of this microstructure are shown in Fig. 2.1.

Figure 2.1 Independent realizations of the periodized random microstructure.

Phase-eld formulation

Using the phase-eld framework to simulate crack propagation, variation of (1.26) with respect to u and d leads to the coupled equations:

   g c (d -2 ∆d) -2(1 -d)H(ε) = 0 , ∇ • σ(u, d) = 0 , (2.1)
Where ∇•(.) is the divergence operator. It should be noted that the elastic energy density Ψ uses the model (1.46), the crack description follows the model (1.38) and a displacement boundary problem is only considered (F * = 0). H is a strain density history function, used to prescribe damage irreversibility [Miehe et al., 2010a]:

H(x, t) = max τ ∈ [0,t] Ψ + (x, τ ) . (2.2)
For an isotropic medium, the stress tensor for the damaged material reads as

σ(u, d) = ∂Ψ ∂ε = g(d)1(tr(ε + ))[λT r(ε) + 2µε + ] + 1(tr(ε -))[λT r(ε) + 2µε -] .
The system of equations ( 2.1) is complemented by the following boundary conditions

     u = u D on ∂Ω u , σn = F * on ∂Ω F , ∇d • n = 0 on ∂Ω , (2.3)
where ∇(.) is the gradient operator, u D and t N are prescribed vector elds of displacements and tractions, and n is the outward-pointing normal vector on ∂Ω (see Fig. 1.14 (a)).

The classical weak form associated with Eq. ( 2.1) is given by

       Ω 2H + g c dδd + g c ∇d • ∇(δd) dΩ = Ω 2HδddΩ , Ω σ(u, d) : ε(δu)dΩ = ∂Ωt t N • δudS , (2.4)
where (δd, δu) are test functions belonging to appropriate functional spaces. It should be noticed that in the phase eld method, crack propagation is described through the evolution of the damage eld x → d(x), which is updated at each load step by solving the coupled equations (2.1) (under Eq. ( 2.3)). Algorithmic details on the method can be found in, e.g., [Miehe et al., 2010a[START_REF] Nguyen | On the choice of parameters in the phase eld method for simulating crack initiation with experimental validation[END_REF].

Notice that the denition of Ψ ± in the phase eld formulation (see (1.46)) encapsulates the stochastic aspect of the propagation, since the elasticity eld corresponds (at microscale) to a realization of the elasticity tensor random eld {C(x), x ∈ Ω} given by

C(x) = Np i=1 1 Ωp (x)C p (λ p , µ p ) , (2.5)
where {1 Ωp (x), x ∈ Ω} and C i are the indicator function and elasticity tensor of phase p (assumed to be isotropic here), occupying the domain Ω p . The couple of Lamé coecients (λ p , µ p ) is constant in each phase and follows the random spatial distribution of inclusions. Respectively we dene the toughness random eld as:

g c (x) = Np i=1
1 Ωp (x)g c,p .

(2.6)

For random microstructures, the indicator functions are, indeed, non-Gaussian random elds:

in practice, it is thus required to proceed to Monte Carlo simulations of these elds, and to solve the coupled elasticity-phase-eld problem for each realization of Ω. This strategy allows the variability in crack paths (and consequently, in the nonlinear part of the macroscopic response)

to be simulated at microscale, as illustrated in the next section.

Statistical analysis of crack trajectories at microscale

In the simulations presented throughout this chapter, the constitutive materials are assumed isotropic, and the bulk and shear moduli are denoted as (k m , µ m ) and (k i , µ i ) for the matrix phase (N p = 2) and inclusions (N p = 1), respectively. Accordingly, g c,m and g c,i denote the toughness of the matrix and inclusions. These properties are chosen such that the mechanical

contrast α satises α = k i /k m = µ i /µ m = g c,i /g c,m
, where the properties of the matrix are taken as

k m = 175 [GPa], µ m = 81 [GPa] and g c,m = 2.7 × 10 -3 [kN.mm -1 ]. Dirichlet boundary conditions are applied in the form u D (x 1 , 0) = 0 and u D (x 1 , L) = u D e (1) , for 0 x 1 L, in which e (1) = (1, 0) is the rst vector of the canonical basis in R 2 (note that ∂Ω N = ∅) and u D ∈ [0, u D ] (the nal displacement is taken for u D = 2 × 10 -2 [mm]); see Fig. 2.2. The incremental displacement value ∆u D = 2×10 -5 [
mm] is selected, and an initial crack is positioned as described in Fig. 2.2. These boundary conditions correspond to a pure shear loading (see Fig. 2.2), and the evolution of the damage eld and displacement-force curve (associated with the microstructural sample shown in Fig. 2.2) can be seen in Figs. 2.3 and 2.4, for = 0.0075 [mm] (in the two phases) and α = 10. The nite element mesh adapted to this value of contains about 150, 000 (linear triangular) elements, with a mesh size comprised between 3 × 10 -4 and 1.5 × 10 -2 [mm]. As expected given the selected contrast in toughness, the crack exclusively propagates within the matrix phase.

Since the crack path Γ is, by denition, identied as the collection of points x Γ for which d ≈ 1, and upon restricting the analysis to congurations containing a single crack, the variability in the crack propagation generated by the underlying microstructural randomness can be observed by considering the stochastic process {x In these gures,

Γ 1 (x 2 ), x 2 ∈ X 2 }, in which X 2 ⊆ [0, 0.5]. Likewise,
u D → F (u D ) = E{F (u D )} and u D → σ F (u D ) = E{F (u D ) 2 } -E{F (u D )} 2
represent the mean and variance functions for the macroscopic force, and the statistical estimators and envelopes are obtained using 100 independent realizations of the microstructure. Similar notations are used for studying the second-order properties of the process {x Γ 1 (x 2 ), x 2 ∈ X 2 }.

Mesoscopic Modeling of Crack Propagation

In this section, we propose a simplied model of crack propagation at a mesoscopic scale where the piecewise constant elds of elastic properties dening the detailed microstructure are replaced by smooth approximations.

Phase-eld formulation

In order to characterize crack propagation using coarse descriptors, we specically introduce an upscaled version of the elasticity-phase-eld problem as follows:

   ∇ • σ(u, d) = 0 , g c (d -2 ∆d) -2(1 -d) Ψ + (ε) = 0 , (2.7)
It should be noted that the choice to use the same formulation as in micro scale has been made.

Only the parameters of the problem change. This study is part of an empirical, exploratory approach which aims to become an operational tool (insofar as there is no homogenization theory in fracture).

where the stress tensor

σ(u, d) = g(d) ∂ Ψ + (ε(u)) ∂ε(u) + ∂ Ψ -(ε(u)) ∂ε(u) (2.8)
is here expressed as a function of a new, mesoscopic stored energy function

Ψ ± (ε) = ε ± : C : ε ± , (2.9)
in which x → C(x) is a mesoscopic elasticity tensor eld, here chosen as isotropic (so that the negative/positive decomposition of Miehe can be used), g c represents an equivalent toughness for the mesoscale medium and is the characteristic length associated with the regularized description at the mesoscale. Note that the restriction related to isotropy could be alleviated by using the recent framework proposed in [START_REF] He | Closed-form coordinate-free decompositions of the two-dimensional strain and stress for modeling tensioncompression dissymmetry[END_REF] (see Section 1.3.4.2): this decomposition was not explored hereafter, due to time constraints and to the fact that the assumption of isotropy is, indeed, reasonably accurate for the system under investigation (see Section 2.3.2.3).

Note that (i ) this mesoscale model diers from a fully homogenized one, as it maintains some statistical uctuations raised by microstructural randomness; (ii ) the formulation remains predictive to study crack propagation, in contrast with an approach that would describe cracks through rst-order, averaged characteristics (such as crack density).

In what follows, the denition of the mesoscopic elasticity and equivalent toughness is investigated through a two-step methodology:

First, the denition of the elasticity eld x → C(x) is achieved using a moving-window upscaling approach under dierent types of boundary conditions. This point is discussed in Section 2.3.2.1.

Second, the denition of the toughness g c is addressed in Section 2.3.3, where a statistical inverse problem involving the peak force at the macroscopic scale is introduced.

The results from the rst step will be used, in Section 2.4.1, to construct a stochastic surrogate for the elasticity eld. This model will enable us to draw additional samples of the mesoscopic elasticity eld without having recourse to the homogenization solver.

Construction of mesoscopic elasticity

In this section, we dene the technique used to construct a (smooth) mesoscopic denition of the heterogeneous elastic medium from fully detailed realizations of microstructures. Such an approach has been extensively discussed over the past two decades. Within a multiscale setting, this can be achieved by using a local homogenization (see [Ostoja-Starzewski, 2008] for a survey) or a ltering approach [START_REF] Yvonnet | A consistent nonlocal scheme based on lters for the homogenization of heterogeneous linear materials with non-separated scales[END_REF], Bignonnet et al., 2014, Tran et al., 2016]. For averaging-type upscaling, kinematic and static uniform boundary conditions (which are denoted by KUBC and SUBC hereinafter) can be considered [Ostoja-Starzewski, 1998] and provide bounds for the eld of apparent tensors [Huet, 1990, Hazanov andHuet, 1994]. Alternatively, periodic boundary conditions (PBC) can be invoked, especially when a fast convergence toward the eective properties is sought; PBC were employed in the so-called moving-window approach [START_REF] Graham | Non-gaussian simulation of local material properties based on a moving-window technique[END_REF], for instance. In the sequel, kinematic uniform boundary conditions (KUBC) and static uniform boundary conditions (SUBC) are selected.

While square-shaped domains are typically used in the literature of homogenization, a circular moving window is considered to prevent the generation of spurious anisotropic features at the mesoscale created by the corners of a square window (see, e.g., [START_REF] Bignonnet | Macroscopically consistent non-local modeling of heterogeneous media[END_REF]), see also [START_REF] Salmi | Apparent and eective mechanical properties of linear matrix-inclusion random composites: Improved bounds for the eective behavior[END_REF] for the consideration of, and comparison with, alternative boundary conditions based on Voronoi cells. A method for obtaining a smoothed, equivalent eld with uctuations is described in the following.

Methodology

Let Ω x be a circular domain of radius R, centered at x ∈ Ω, with boundary ∂ Ω x . For one realization of the microstructure, { C(x), x ∈ Ω} is obtained by performing a homogenization locally in Ω x , for both KUBC and SUBC (see Fig. 2.6). As x moves within Ω, we obtain a smooth, equivalent medium characterized by wavelengths associated with the radius R (see the illustration of this process in Fig. 2.7 for dierent radii R). We recall that KUBC correspond to the following boundary conditions:

u (ij) D (z) = [E (ij) ]z , ∀z ∈ ∂ Ω x , (2.10)
where in the present 2D plane strain context, the indices i and j run over {1, 2}, leading to:

[E (11) ] = 1 0 0 0 , [E (22) ] = 0 0 0 1 , [E (12) ] = 1 2 0 1 1 0 . (2.11)
The strategy is schematically depicted in Fig. 2.6. For SUBC, traction vectors are applied in the Figure 2.6 Dening mesoscale elds of stiness through homogenization (case of KUBC).

form

t N (z) = [Σ (ij) ] n(z) , ∀z ∈ ∂ Ω x ,
(2.12)

where n(z) is the outward-pointing normal vector at point z ∈ ∂ Ω x . Combinations of indices similar to those introduced for KUBC are considered, with

[Σ (11) ] = 1 0 0 0 , [Σ (22) ] = 0 0 0 1 , [Σ (12) ] = 1 2 0 1 1 0 . (2.13) Let { C KU BC (x), x ∈ Ω} and { C SU BC (x),
x ∈ Ω} be the random elds of mesoscopic elasticity tensors dened under the aforementioned boundary conditions. Note that when x approaches the boundary ∂Ω, the realization of the microstructure is virtually replicated, by periodicity, and the mesoscopic tensor is still well dened.

Denoting by R the radius of the inclusions, the ratio R/R plays an important role in dening a continuous transition from the microscale ( R/R → 0 + ) to the macroscale ( R/R → +∞).

Additionally, it species the level of anisotropy exhibited by the local apparent elasticity tensor, ranging from microscopic isotropy to mesoscopic anisotropy, and then to macroscopic isotropy (in the present case). This aspect turns out to be critical for the phase-eld approach, since the constitutive model used for the damaged material is based on the isotropy of the background media (see Eq. 2.9). In particular, the local isotropy of the underlying mesostructure is implicitly assumed in the presented formulation. The denition of the isotropic approximation for the mesoscale elasticity tensor eld is addressed in the following.

Construction of the mesoscopic isotropic approximations

Let { C iso BC (x), x ∈ Ω} denote the isotropic approximation of the eld { C BC (x), x ∈ Ω}, where the subscript BC refers to the type of boundary conditions under consideration (KUBC or SUBC).

In two-dimensional elasticity, and assuming plane strain conditions, the Voigt-type matrix representation {[ C iso BC (x)], x ∈ Ω} of the aforementioned eld is given by

[ C iso BC (x)] =     k BC (x) + 4 3 µ BC (x) k BC (x) -2 3 µ BC (x) 0 k BC (x) -2 3 µ BC (x) k BC (x) + 4 3 µ BC (x) 0 0 0 µ BC (x)     , ∀ x ∈ Ω , (2.14)
where { k BC (x), x ∈ Ω} and { µ BC (x), x ∈ Ω} are the random elds of three-dimensional bulk and shear moduli dening the isotropic approximation, in plane strain elasticity, of the actual

(anisotropic) elastic tensor [ C BC (x)].
To obtain these coecients, we minimize the distance (in the sense of the metric dened below) between [ C BC (x)] and [ C iso BC (x)] (see [Guilleminot andSoize, 2012b, Tran et al., 2016] and the references therein):

( k BC (x), µ BC (x)) = argmin k>0, µ>0 [ C BC (x)] -[ C iso BC (x)] 2 F , (2.15) in which [ C BC ] =    C 11 C 12 C 13 C 12 C 22 C 23 C 13 C 23 C 33    , (2.16)
where the Voigt's notation has been used for the dierent components of the tensor and • F is the Frobenius norm.

The optimization problem dened by Eq. ( 2.15) can be solved by a direct dierentiation of the cost function, and the mesoscopic moduli of the isotropic approximation are found as

k(x) = 1 60 (11 C 11 (x) + 11 C 22 (x) -4 C 33 (x) + 38 C 12 (x)) (2.17) and µ(x) = 1 5 ( C 11 (x) + C 22 (x) + C 33 (x) -2 C 12 (x)) , (2.18)
where the subscript BC has been dropped for notational convenience (this convention will be used in the sequel when no confusion is possible).

Identication of mesoscale elastic properties

In this section, we rst analyze the error generated by the isotropic approximation at mesoscale.

A few fundamental properties of the eld thus constructed are then investigated.

Isotropic mesoscopic approximation

The eect of the mesoscopic resolution on the isotropic approximation (for KUBC) is qualitatively

shown on a single realization in Fig. 2.7. As expected, the eld becomes more homogeneous as the ratio R/R increases. In order to further assess the relevance of the isotropic approximation,

A BC (x) = [ C BC (x)] -[ C iso BC (x)] F [ C BC (x)] F , ∀ x ∈ Ω (2.19)
The graphs of the elds of mean and standard deviation evaluated for a coarse mesoscopic grid and their isotropic approximations remain small in mean and variance, and that the error is larger in the case of KUBC. The approximation in the set of isotropic tensors is satisfactory and allows the phase-eld approach to be applied with an isotropic background medium. In this rest of this chapter, we will consider the characterization and simulation of the non-Gaussian elds

{ k BC (x), x ∈ Ω} and { µ BC (x), x ∈ Ω} for R/R = 3.

Statistical analysis on mesoscopic elasticity

The graphs of the rst-order marginal probability density functions for the bulk and shear moduli are shown in Fig. 2.10. The well-known ordering with respect to boundary conditions is observed almost surely (that is, k SU BC k KU BC and µ SU BC µ KU BC for each microstructural sample), and it is seen that the level of statistical uctuations associated with KUBC is larger than for SUBC.

The estimated normalized correlation functions along e (1) and e (2) for the elds of bulk and shear moduli (for SUBC) are shown in Fig. 2.11. In these gures, the notation τ i → R data k (τ i ) indicates that the correlation function of the bulk modulus random eld is evaluated along the unit vector e (i) (a similar notation is used for the shear modulus). It is seen that the dierences between the correlation functions for the two random elds are almost indistinguishable, due to the very strong cross-correlation between the two properties. Moreover, it can be observed that the correlation rst decreases over the range [0, L/2] (with L = 1) and then starts increasing on [L/2, L], in accordance with the periodicity of the underlying background medium. This information will be used to select an appropriate form of the correlation functions for the random eld models, constructed in Section 2.4.1.1.

Mesoscopic Toughness

This part is concerned with the identication to the mesoscopic toughness, solving an inverse problem. A validation study involving both the mechanical response (under some given, macroscopic loading) and crack path variability is then undertaken. ) and e (2) , estimated from the simulated data, for the random elds of bulk (left) and shear (right) moduli.

Inverse problem strategy

Let us now turn to the identication of the fracture and phase-eld parameters g c and at the mesoscale (see Section 2.3.1). In this work, the characteristic length of the mesoscopic medium is set to be equal to the characteristic length at the microscale, previously denoted by , and g c is assumed constant. These choices are supported by a set of parametric studies, unreported hereinafter for the sake of conciseness, and by the numerical results presented in [START_REF] Nguyen | Identication of fracture models based on phase eld for crack propagation in heterogeneous lattices in a context of non-separated scales[END_REF] for similar propagation regimes. It is worth noticing that while the length scales are taken similar in both the microscopic and mesoscopic phase eld equations (given by Eq. (2.7)), the mesh involved in the mesoscale description only needs to be rened in the vicinity of the crack.

An adaptive meshing strategy can then be deployed to substantially reduce the computational cost (as opposed to the microscopic description in which the mesh must be uniformly ne over the entire microstructure). The mesoscopic toughness parameter g c is next identied by solving an inverse problem involving the peak force F max = max u D F (u D ). More precisely, g c is calibrated by imposing a match between the mean value F max of F max , estimated with 500 independent microstructural samples and the ne-scale elasticity-phase-eld simulations (detailed in Section 2.2.2), and the mean value determined with the mesoscopic description introduced in Section 2.3.1, denoted by F max . In the latter description, samples of the elasticity random eld are estimated through the mowing-window homogenization procedure, and the mesoscopic toughness then appears as the unique unknown parameter. Let F max ( g * c ) be the mean peak force associated with the candidate value g * c for the mesoscopic toughness (that is, by substituting g * c for g c in Eq. (2.7)). An optimal value can thus be dened by minimizing the relative error function

J( g * c ) = |F max -F max ( g * c )| F max (2.20)
over the admissible set [g c,m , g c,i ]:

g c = argmin g * c ∈ [gc,m, g c,i ] J( g * c ) .
(2.21)

Identication results and validation study

Since the mesoscopic elasticity eld depends on the boundary conditions applied, the optimization problem dened in the previous section must be solved independently for KUBC and SUBC. The The optimal values are obtained as

g c = 1.022 × g c,m ≈ 2.75 × 10 -3 [kN.mm -1 ] (2.22)
for KUBC and is shown in the right panel in Fig. 2.12; see also Fig. 2.13 for a comparison on crack path (for a given microstructural samples). This gure shows that the mesoscale formulation identied under SUBC provides a fairly accurate estimate of the mean macroscopic response. This conclusion similarly holds for KUBC, although this type of boundary conditions leads to a stiening of the response, in accordance with the fact that apparent tensors obtained under KUBC constitute upper bounds for the mesoscopic elasticity (it becomes surrealistic in particular when R/R → 1).

In both cases, the variability is underestimated due to ltered elasticity uctuations. These eects can clearly observed in Fig. 2.14, where the mean and standard deviation on crack paths are reported for the two types of boundary conditions. It is seen that the cracks paths obtained with the mesoscopic formulations are localized near the mean crack path at microscale, for both KUBC and SUBC. Not surprisingly, the mesoscopic-based crack paths exhibit a variability that is much smaller than the one obtained at microscale (where the crack trajectory is constrained by the radius of the inclusions). While these results support the relevance of the formulation, the latter necessitates solving a very large number of homogenization problems to represent the uctuations of the elasticity eld at the mesoscopic scale. In the next section, we address the construction of a stochastic model that enables the elasticity tensor random eld to be sampled in a robust manner.

Γ micro Γ ups KU BC Γ micro Γ ups SU BC

Stochastic Modeling of Mesoscopic Crack Propagation

In this section, we rst address the construction and identication of a stochastic representation for the mesoscopic elasticity tensor random eld. We subsequently discuss the validation of the proposed model using direct and indirect data.

Stochastic modeling

The construction of stochastic models for random elds of elasticity tensors exhibiting arbitrary material symmetries has been investigated in [START_REF] Guilleminot | Stochastic model and generator for random elds with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF], Guilleminot and Soize, 2013a, Staber and Guilleminot, 2017] using an information-theoretic formulation [Soize, 2006] (see [START_REF] Guilleminot | Non-Gaussian Random Fields in Multiscale Mechanics of Heterogeneous Materials[END_REF] for a survey), and in [START_REF] Malyarenko | A random eld formulation of hooke's law in all elasticity classes[END_REF]] using a spectral expansion. In what follows, we seek to dene a representation that is consistent, on the one hand, with the observations drawn from the analysis of the samples obtained through the moving-window upscaling and, on the other hand, with the theoretical results derived in the aforementioned references.

Construction of random eld models

For a given type of boundary conditions, let { k BC (x), x ∈ Ω} and { µ BC (x), x ∈ Ω} be the random elds of three-dimensional bulk and shear moduli as dened in Eq. (2.14).

Let {Ξ(x) = (Ξ 1 (x), Ξ 2 (x)), x ∈ Ω} be a bivariate Gaussian eld with statistically independent, normalized components. These components are dened by the correlation functions (x, y) → R Ξ 1 (x, y; α (1) ) and (x, y) → R Ξ 2 (x, y; α (2) ) (α (i) is a model parameter dened for Ξ i ). The non-Gaussian random elds of elastic moduli are then dened through the nonlinear transformations

k BC (x) = F -1 G(p k ,q k ) F N (0,1) (Ξ 1 (x)) (2.24) and µ BC (x) = F -1 G(p µ ,q µ ) F N (0,1) (ρ Ξ 1 (x) + 1 -ρ 2 Ξ 2 (x)) , (2.25)
where F -1 G(p,q) is the inverse cumulative distribution function of the Gamma law, ensuring the pos- itivity of the modulus and with (shape and scale) parameters p and q (note that the dependence of these parameters on the KUBC and SUBC boundary conditions of homogenization procedure is not reported for notational convenience), F N (0,1) is the cumulative distribution function of the standard Gaussian law and ρ controls the dependance between k BC (x) and µ BC (x).

Given the stationarity and the form of the correlation functions estimated for the random elds of elastic moduli (see Fig. 2.11), the following separable form is retained:

R Ξ i (x, y; α (i) ) = r(τ 1 ; α (i) 1 ) × r(τ 2 ; α (i) 2 ) , ∀ τ ∈ ([0, L]) 2 , i ∈ {1, 2} , (2.26)
where τ j = |x j -y j | is the lag distance along e (j) , j ∈ {1, 2}, and the one-dimensional normalized correlation function τ → r(τ ; α) is dened as r(τ ; α) = exp -2 α 2 sin 2 πτ L .

(2.27)

In Eq. (2.27), α is a model parameter related to the internal length

L = L/2 0 |r(τ ; α)| dτ , (2.28)
which is interpreted, in the periodic setting under consideration, as the spatial correlation length of the Gaussian random eld along the associated basis vector (e.g., along e (1) if the function

τ 1 → r(τ 1 ; α (i)
1 ) is considered). It can be shown that the correlation length L < L/2 reads as

L = L 2 exp{-α -2 }I 0 (α -2 ) , (2.29)
where I 0 denotes the zero-order modied Bessel function. It should be noted that by construction, one has < L/2. The graph of τ → r(τ ; α) is shown in Fig. 2.15 for dierent values of α.

The following properties can easily be deduced.

The rst-order marginal probability measure is a bivariate Gamma law [Moran, 1969, Arnst andPonthot, 2014], which is consistent with previous results derived within the framework of information theory (see [START_REF] Guilleminot | Non-Gaussian Random Fields in Multiscale Mechanics of Heterogeneous Materials[END_REF] and the references therein). The mean values k BC = E{ k BC (x)} and µ BC = E{ µ BC (x)} read as

k BC = p k × q k , µ BC = p µ × q µ , (2.30)
and the coecients of variation are given as

δ k BC = 1 √ p k , δ µ BC = 1 √ p µ .
(2.31)

These properties are chosen independent of location x, owing to the stationarity of the random elds.

The elds of stiness and compliance tensors are of second-order:

E{ [ C iso BC (x)] 2 F } < +∞ , E{ [ C iso BC (x)] -1 2 F } < +∞ , ∀ x ∈ Ω .
(2.32) hence ensuring that the stochastic linear elastic boundary value problem is well posed [Soize, 2006].

The random elds { k BC (x), x ∈ Ω} and { µ BC (x), x ∈ Ω} are mean-square continuous and mean-square dierentiable.

From a computational standpoint, the underlying Gaussian elds are sampled using a truncated Karhunen-Loève expansion. In order to reduce the associated computation time, the random elds are sampled on a grid that is coarser than the one used to solve the elasticityphase-eld problem at the mesoscale. Realizations of elds are then obtained by interpolating on the ne mesoscopic grid. To that end, the coarse mesoscopic grid is specically dened so that the correlation structure is properly discretized. In the results presented hereinafter, the coarse mesh includes six Gauss points per correlation length, along each direction.

Identication of the elasticity random eld

The probabilistic model involves two sets of parameters controlling (i) the joint probability density function of the elastic moduli at a given location, and (ii) the correlation structure of the underlying Gaussian elds. The rst set of parameters gathers (p k , q k ) and (p µ , q µ ) (or equivalently ( k BC , δ k BC ) and ( µ BC , δ µ BC ), in view of Eq. ( 2.30) and (2.31)), as well as the the coecient of correlation ρ. These hyperparameters can be estimated from the database using standard statistical estimators, here with 500 sample realizations (which ensures the convergence of the estimators) and for the same resolution R/R = 3:

For KUBC, we found ρ = 0. Since the transformations given by Eqs. (2.24) and (2.25) are nonlinear, the correlation functions associated with the random elds of elastic moduli cannot be inferred explicitly. In this case, the identication of the vector-valued hyperparameters α (1) and α (2) is performed through the following two-step procedure. Let the correlation functions of { k BC (x), x ∈ Ω} and { µ BC (x), x ∈ Ω} be written as τ → R model k 2) )), respectively: this notation emphasizes the underlying dependence on the parameters of the Gaussian elds (see Eqs. (2.24) and (2.25)). By a slight abuse of notation, these correlation functions will also be denoted as

(τ ; α (1) ) and τ → R model µ (τ ; (α (1) , α (
τ j → R model k (τ j ; α (1) j ) and τ j → R model µ (τ j ; (α (1) j , α (2)
j )) when evaluated along e (j) , j ∈ {1, 2}. In a rst step, the components of α (1) are identied, for a given type of boundary conditions, as

α (1) 1 = argmin α > 0 J (1) 1 (α) , α (1)
2 = argmin α > 0 J (1) 2 (α) , (2.37)
where the cost functions are given by J

(1)

1 (α) = L/2 0 (R data k (τ 1 ) -R model k (τ 1 ; α)) 2 dτ 1 1/2
(2.38) and

J

(1)

2 (α) = L/2 0 (R data k (τ 2 ) -R model k (τ 2 ; α)) 2 dτ 2 1/2
.

(2.39)

In a second step, the hyperparameters controlling the correlation structure of the Gaussian random eld {Ξ 2 (x), x ∈ Ω} are identied as

α (2) 1 = argmin α > 0 J (2) 1 (α) , α (2)
2 = argmin α > 0 J (2) 2 (α) , (2.40) in which J (2) 1 (α) = L/2 0 (R data µ (τ 1 ) -R model µ (τ 1 ; (α (1) 1 , α)) 2 dτ 1 1/2
(2.41) and J

(2)

2 (α) = L/2 0 (R data µ (τ 2 ) -R model µ (τ 2 ; (α (1) 2 , α)) 2 dτ 2 1/2
.

(2.42)

In Eqs. (2.41) and (2.42), the values of α

(1)

1 and α

(1) 2 are those obtained within the rst step of the methodology (see Eq. (2.37)). Above, the estimations of the correlation functions associated with the stochastic model are obtained as follows. For given values of the hyperparameters, a set of 500 independent realizations of the random elds is rst generated on a coarse grid with equidistant points. Usual statistical estimators are then used to estimate a set of correlation functions, indexed by the reference point (which is the point with respect to which the lag vector τ is dened). Spatial averaging over properly selected reference points is nally applied to improve the quality of the estimations. The optimal values are found as α (1) = (0.4624, 0.4574) , α (2) = (0.4043, 0.4014) (2.43) for the elasticity random elds identied under KUBC, and α (1) = (0.4654, 0.4604) , α (2) = (0.4694, 0.4654) It is seen that the calibrated model allows the decays of the correlation functions to be accurately reproduced, which is key to mimicking the mesoscopic elasticity (and in particular, the frequency of sample path oscillations that has a substantial impact on the crack paths in the phase-eld formulation at mesoscale).

(

Comparison of stochastic microscopic and mesoscopic modeling

The comparison is then carried out regarding the crack path and the force response, to determine whether the statistical structure of the fracture model and generated elasticity elds is sucient to represent the crack test.

Results on the crack propagation

The variability in crack paths can be observed in Fig. 2.20 for the reference microscale-based computations and the mesoscopic formulation dened with either locally-homogenized microstructural samples (following the approach detailed in Section 2.3.2.1) or the elasticity eld stochastic model. It can be observed that the crack paths corresponding to a description at microscale the error generated by the proposed model-based, mesoscopic formulation can be characterized as

Γ micro Γ ups KU BC Γ mod KU BC Γ micro Γ ups SU BC Γ mod SU BC
ε = |E{ u D 0 F (u D ) 2 du D } 1/2 -E{ u D 0 F (u D ) 2 du D } 1/2 | E{ u D 0 F (u D ) 2 du D } 1/2 , (2.45)
where the stochastic process { F (u D ), u D ∈ [0, u D ]} implicitly depends on the boundary conditions applied at mesoscale. The right-hand side term in Eq. 2.45 can be estimated through Monte Carlo simulations, and the relative error remains small for the two types of boundary conditions, with ε ≈ 1.6% for KUBC and ε ≈ 2.2% for SUBC. The error for the prediction of the mean peak force can be characterized by

ε max = |F max -F max |/F max . (2.46)
The error measure is given by ε max ≈ 0.084% for KUBC and ε max ≈ 0.037% for SUBC, showing that an accurate prediction of the mean peak force can be obtained with the two types of boundary conditions. Finally, the fact that the mesoscopic approach underestimates the variability in the macroscopic response is expected, given the nature of the propagation at microscale. In this context, the deviation from the mean crack path is more contained than in the microscopic simulations where the crack essentially propagates around the heterogeneities. Depending on the application of interest, one possible way to compensate for this intrinsic eect could be to adopt a goal-oriented strategy where (ctitious) anisotropic uctuations are incorporated into the mesoscopic elasticity eld (by means of a generalized stochastic model; see [Guilleminot andSoize, 2012a, Guilleminot and[START_REF] Guilleminot | Stochastic model and generator for random elds with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF]), and where the hyperparameters of the stochastic model are calibrated by solving a statistical inverse problem on the macroscopic response (using an appropriate identication metric that is sensitive to both the mean and variance along the macroscopic loading path).

Conclusion

A stochastic approach to model crack propagation in random media has been proposed in this chapter. The formulation relies on a phase-eld formulation where material coecients are dened and identied through multiscale computations. Monte-Carlo simulations were rst performed using a description at the microscopic scale. These computations enable the characterization of subscale-induced randomness on the macroscopic response of the domain and were subsequently used as reference results to assess the relevance of the framework. The denition of the mesoscopic parameters was then addressed. The elasticity eld at mesoscale was specically dened as the isotropic approximation of spatially dependent homogenized tensors, obtained by means of a moving-window upscaling approach (under kinematically and statically uniform boundary conditions). The (deterministic) mesoscopic toughness was identied by solving an inverse problem related to the mean peak force. It is shown that the formulation under statically uniform boundary conditions allows for an accurate prediction of the mean elastic response and mean peak force. In contrast, kinematically uniform boundary conditions generate a stiening of mesoscale elasticity, in accordance with theoretical results derived elsewhere. An informationtheoretical probabilistic model for the elasticity random eld was then constructed and allows for a fast, robust sampling of mesoscopic elasticity. The results obtained by feeding this stochastic surrogate model into the phase-eld formulation were nally compared with those corresponding to the full-scale, microscopic model. It is shown, in particular, that the model-based, mesoscopic elasticity-phase-eld formulation associated with statically uniform boundary conditions allows for an accurate prediction of both the mean elastic response and mean peak force. Extensions of the present framework are further discussed in the perspective chapter of this thesis (see Chapter 5.5) .

Chapter 3

Phase-eld formulation for nite strains and shrinkage

Introduction

In this chapter, a phase-eld model for nite strains is presented. This formulation is motivated by the large deformations that are experimentally observed during drying tests on clay samples, as discussed in chapter 4. Several works have been devoted to the extension of the phase-eld method (PFM) to nonlinear behaviors or nite strains in the literature. In [Miehe andSchänzel, 2014, Miehe et al., 2015], the PFM has been adapted to model polymers with thermoplastic properties. The energy was specically described through a Neo-Hookean model coupled with a damage term. In [START_REF] Ambati | A phase-eld model for ductile fracture at nite strains and its experimental verication[END_REF], the authors also used a Neo-Hookean model which diers from the aforementioned one by decomposing the potential into deviatoric and volumetric parts.

The damage aects this energy according to the decomposition model proposed by [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] in the case of small deformations. Theirs works was carried out on the quasi-static ductile fracture model and compared with experimental data, from the literature, for steel and aluminum samples. Good predictions on the force response in the plastic regime, and also on the crack path, have been presented. In [Raina and Miehe, 2016, Gültekin et al., 2018], another model was developed for soft tissue materials by incorporating in the energy an additional potential taking into account the preferential brous orientation in the material.

In this chapter, we present our phase-eld formulation for nite strains, reviewing both the formulation and its numerical implementation. These aspects are mostly borrowed from the literature. In addition, an extension to drying shrinkage is proposed, together with an analysis on parameter sensitivity, solving strategies, and a comparison with the linear case. For convenience, the nite strains and the shrinkage extension are studied independently. The chapter 5 will deal with the two models together.

Finite strain kinematics and mechanical modeling

Background in nite elasticity

Let X be the position of a material point in the reference conguration Ω, and let x = φ(X) denote the position of this material point in the actual conguration Ω t , with φ is the deformation map. The second-order deformation gradient tensor F is dened by:

F = ∇ X (φ(X)) , (3.1)
where ∇ X (.) denotes the gradient operator with respect to the reference conguration Ω. The right Cauchy-Green deformation tensor is dened as

C = F T F , (3.2)
and the Green-Lagrange strain tensor is given by

E = 1 2 (C -I) , (3.3)
where I is the second-order identity tensor. In the present work, we adopt the linear Saint-Venant-Kirchho material constitutive law dened by the strain density function:

Ψ(E) = 1 2 E : C : E . (3.4)
For this choice of strain density function, the associated second Piola-Kirchho stress tensor reads as

S =

∂Ψ(E) ∂E = C : E . where P = F S is the rst Piola-Kirchho stress, and N is the outward-pointing unit vector on boundary ∂Ω N in reference conguration. The total energy of the system is then given by:

E(u, d) = Ω Ψ(E(u), d) dΩ + g c Ω γ(d, ∇(d)) dΩ - ∂Ω N t • u dΓ , (3.8)
where Ψ is the degraded strain energy function and γ denotes the crack density function. In what follows, we only consider the case of local traction strain elds and do not take into account closure of cracks, for the sake of simplicity. Formulations considering the asymmetry of fracture (by spectrally decomposing in this context the local deformation operator F ) can be found elsewhere in the literature; see, e.g., [START_REF] Hesch | Thermodynamically consistent algorithms for a nite-deformation phase-eld approach to fracture[END_REF]. We thus dene

Ψ(E, d) = g(d)Ψ(E) , (3.9)
where g(d) is the degradation function described in section 1.3.4.2 (see [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]).

The second Piola-Kirchho stress tensor is then given by

S(E, d) = ∂ Ψ(E, d) ∂E = g(d) ∂Ψ(E) ∂E = g(d)(C : E) .
The function γ takes the form

γ(d, ∇d) = 1 2 1 d 2 + ∇(d) • ∇(d) .
(3.10)

The coupled problem is nonlinear due to geometric nonlinearities. In this work, the solution is obtained by using a standard staggered scheme recalled in section 3.3.1.

Computational aspects

Staggered resolution strategy

As previously mentioned, a staggered solution strategy is adopted in which the damage and mechanical problems are alternatively solved for each quasi-static load increment. For a given state of damage, the mechanical problem u = argmin

u * ∈ Su E(u * , d)) , (3.11)
where S u is a set of kinematically admissible elds, is solved by using the Newton-Raphson method. The associated weak form can be expressed as :

Ω (P (u), ∇δu)dΩ = ∂Ω N (t(u), δu)dS ,
(3.12)

with P the rst Piola-Kirchho stress. It can also be written as :

Ω S(E(u)), δE(δu, u) dΩ = ∂Ω N (t(u * ), δu)dS , (3.13)
with S the second Piola-Kirchho stress express in Eq. (3.64) and E the Green-Lagrange strain tensor. Regarding the phase-eld problem, a weak form similar to the one obtained at small strains is derived :

Ω g c + 2H (d, δd) + g c (∇d, ∇δd) dΩ = Ω (2H, δd)dΩ , (3.14)
where the history function is given by

H(x, t) = max τ ∈ [0,t] {Ψ (E(x, τ ))} , (3.15)
describing the irreversible process of damage in time and Ψ is dened in (3.4).

The algorithmic structure for this procedure is described in Fig 3 .1.

Figure 3.1 Algorithm used to solve the phase-eld formulation at nite strains.

Linearization of the mechanical problem

The residual associated with the weak form of the mechanical problem in the initial conguration is dened as :

R(u, δu) = Ω S(E(u), d) : δE(u, δu) dΩ - ∂Ω N t(u) • δu dΓ , (3.16)
where d is temporarily frozen (staggered scheme) and δE is the variation of E obtained through

the directional derivative. Let R(u, δu) = R int (u, δu) -R ext (u, δu) = 0 , (3.17)
where R int and R ext are dened as

R int (u, δu) = Ω S(E(u), d) : δE(u, δu) dΩ, R ext (u, δu) = ∂Ω N t • δu dΓ . (3.18)
Here an implicit assumption is made on R ext that there is no tracking force, such as pressure.

By applying the Newton-Raphson method, the linearization of the residual is written

R(u k+1 , δu) ≈ R(u k , δu) + D ∆u R(u k , δu) , (3.19) must be evaluated, where D ∆u R(u k , δu) = D ∆u R int (u k , δu)
∂E ij ∂E kl = 1 2 (δ ik δ jl + δ il δ jk ) and ∂E pp δ ij ∂E kl = δ ij δ pk δ pl = δ ij δ kl , (3.23)
where I is the fourth-order (symmetric) identity tensor and ⊗ denotes the tensor product, the tangent material tensor C t M can be dened as:

C t M (E(u), d) = g(d)[λI ⊗ I + 2µI] . (3.24) Combining D ∆u R(u k , δu) = -R(u k , δu) , (3.25)
the nite element discretization then leads to a linear system

[K u (u)]{∆u} = {F u (u)} (3.26)
to be solved for the increment ∆u. The nite element discretization is presented in the next section.

Finite element discretization

This section is devoted to the nite element discretization of the linearized problem.

Mechanical problem

La solution en déplacement u, la fonction de test δu et le déplacement incrémental ∆u sont interpolés sur chaque élément comme :

u(x) = [N u (x)]{u i e } , (3.27) δu(x) = [N u (x)]{δu i e } , (3.28) ∆u(x) = [N u (x)]{∆u i e } , (3.29)
where [N u (x)] is the matrix of shape functions and the subscript e indicates vectors of i-nodal values. The gradient ∇ X u is dened according to the convention

∇ X u = (u 1,1 , u 1,2 , u 2,1 , u 2,2 ) T (3.30) in 2D and ∇ X u = (u 1,1 , u 1,2 , u 1,3 , u 2,1 , u 2,2 , u 2,3 , u 3,1 , u 3,2 , u 3,3 ) T (3.31)
in 3D, and is computed through (3.37)

∇ X u(x) = [G(x)]{u i e } , ( 3
Expressions for these matrices can be found in Appendix. (A.2.2) for the 2D and 3D cases.

The matrix form of the fourth-order elasticity tensor at small strains is written as

[C] =   λ + 2µ λ 0 λ λ + 2µ 0 0 0 2µ   (3.38)
in 2D plane strain, and as

[C] =         λ + 2µ λ λ 0 0 0 λ λ + 2µ λ 0 0 0 λ λ λ + 2µ 0 0 0 0 0 0 2µ 0 0 0 0 0 0 2µ 0 0 0 0 0 0 2µ         (3.39) in 3D.
Introducing the above discretization in Eq. (3.25), the problem can be written in each element (3.46)

The displacement increment ∆u is thus obtained by solving the linear problem

[K u (u)]{∆u} = {F u (u)} , (3.47)
where [K u (u)] and {F u (u)} are obtained by assembling elementary matrices:

[K u (u)] = Ne A e=1 [k u (u; e)], [F u (u)] = Ne A e=1 [f u (u; e)] ,
(3.48)

where A symbolically denotes here the assembly operator for displacement eld problem. Note that the construction of all operators is based on the notation presented in [START_REF] Bonnet | Analyse des solides déformables par la méthode des éléments nis[END_REF].

History functional

The computation of the history functional (see Eq. (3.68)) requires the Green-Lagrange matrix evaluated with displacement u:

E = 1 2 (∇ X u + ∇ X u T + ∇ X u T ∇ X u) .
(3.49)

The function is then computed on the integration points according to:

H(x, t) = max τ ∈ [0,t] λ 2 (tr (E(x, τ ))) 2 + µtr (E(x, τ )) 2 .
(3.50)

Damage problem

We proceed similarly for the damage problem dened identically as in small strain. So we recall the damage scalar eld d and the associated test function δd are interpolated as The above integrals are computed using a standard quadrature rule (see Eq. (A.9)).

d(x) = [N d (x)]{d i e } (
[k d (e)] = Ωe g c + 2H [N d (x)] T [N d (x)] + g c [B d (x)] T [B d (x)] dΩ ,
Finally, the unknown damage eld vector d is obtained by solving the global problem:

[K d ]{d} = {F d } , (3.57)
where [K d ] and {F d } are the assembled matrix and assembled force vector of the global damage problem

[K d ] = Ne A e=1 [k d (e)] , [F d ] = Ne A e=1 [f d (e)] ,

Modeling of hydric shrinkage

In this section, we introduce the large deformation shrinkage model applied to 2D plane strain and 3D cases in order to model the experiments that will be presented in chapter 4, assuming either uniform strain in the thickness of clay the samples, as a rst modeling attempt, or taking into account the full 3D complexity of the experiments and in particular heterogeneous shinkages along the thickness of the samples. Note here that the coupling with the phase eld has been formulated in the retraction model to gather all the equations to deal with the complete problem: nite strain and shrinkage.

Phase-eld formulation including hydric strains

Models coupling nite strains with thermal expansions can be found in, e.g., [Lu andPister, 1975, Erbts andDüSter, 2012]. Here, we pursue a similar approach, substituting hydric parameters for thermal coecients. The total deformation gradient is thus dened through the multiplicative decomposition: (3.59) where F e and F h are associated with the elastic and hydric parts of the strain; see Fig. 3.2.

F = F e F h ,
In the present work, the hydric strain is modeled as a purely volumetric contribution and is dened as

F h = F(h)I , (3.60)
where F is a nonlinear function of the hydric evolution parameter, denoted by h. This function can be identied based on an experimental study; see chapter 5. Below, we simply exemplify the procedure and choose for the sake of illustration. It is known that the Lagrangian strain E corresponding to F is given by

E = E h + F T h E e F h , (3.62)
where E h and E e are the Green-Lagrange strain tensors associated with F h and F e , respectively.

By combining Eq. (3.61) and Eq. ( 3.62), it follows that (3.63) where E h = α h 1 with α h = 1 2 (h 2 -1). Using the constitutive relationship dened in Eq. (3.5), we deduce

E = E h + h 2 E e ,
S(E(u), E h , d) = g(d) h 2 [λtr(E) + 2µE] -[λtr(E h ) + 2µE h ] .
(3.64)

History functional

To take into account the damage only within traction, the strain density function is modied according to

Ψ(E, E h ) = Ψ + (E, E h ) + Ψ -(E, E h ) , (3.65) with Ψ ± (E e ) = λ 2 tr (E e ) ± 2 + µtr E ± e 2 .
(3.66)

The spectral decomposition is applied in this case to E e tensor as:

E ± e = m i=1 < κ i > ± ϕ (i) ⊗ ϕ (i) , (3.67) in which {(κ i , ϕ (i) )} m
i=1 are the pairs of associated eigenvalues and eigenvectors of the strain tensor E e , and < • > ± is the same operator dened in (1.48).

The history function H is actualized in nite strain, on the positive part(see e.g. [START_REF] Areias | Phase-eld analysis of nite-strain plates and shells including element subdivision[END_REF]) and caused only by the total strain tensor E according to

H(x, t) = max τ ∈ [0,t] Ψ + (E (x, τ )) . (3.68)
The equations of the model with hydric shrinkage are summarized as follows.

Strong forms

The cracking problem is dened as two incremental problems (mechanical and damage) solved with a history function updated at each time step τ ∈ [0, t] ; the strong forms associated are dened as:

Mechanical problem:            Div X (P ) = 0 on Ω S(E(u), E h , d) = g(d) h 2 C(E -E h ) E h = α h I u = u D on ∂Ω D P • N = (t) on ∂Ω N (3.69)
where Div X (.) denotes divergence operator with respect to the reference conguration and assuming t well known (no tracking forces).

Damage problem:

   2(1 -d)H(τ ) -gc l {d -∆ X d} = 0 on Ω d(x) = 1 on Γ ∇ X d(x) • n = 0 on ∂Ω (3.70)
where ∆ X (.) denotes the Laplacian operator with respect to the reference conguration. (3.74)

Weak forms

The stiness matrix [k u (u, e)] is decomposed according to:

[k u (u; e)] = [k 1 u (u; e)] + [k 2 u (u; e)] ,

2.7e

-3

[mm] 7.5

-3 7.5 -3 ∆u [mm] 1.0e -5 1.0e -5 h [mm] 1.0e -3 -3.0e -2 1.0e -3 -3.0e -2
Table 3.1 M I and M II problem parameters for mesh convergence.

We use here 6 dierent mesh sizes to study the mesh convergence h ∈ [1.0e -3 : 3.0e -2 ] mm, rened within the areas where the crack propagates. We can observe in Fig. 3.5 the force response of the system for the M I and M II problems for the dierent mesh sizes h, which allows verifying the convergence with respect to the mesh discretization. These results show the convergence with respect to the mesh size. It has been noticed in the literature [Nguyen et al., 2016a] that the convergence of the phase-eld simulation is achieved when the relation between mesh size and the length of regularization respect: h ≤ 2 . In the studied examples, convergence according to the regularization criterion is estimated at L/h ≥ 300. For the next following study, we have chosen the mesh size h = 2 × 10 -3 corresponding to L/h = 500.

Sensitivity with respect to load increments

In the previous tests, a constant load step ∆u = 1.0e -5 has been used (see Table 3.1). In the present example, (see Fig. 3.7), dierent load steps have been used. We can observe the inuence of the load step on the mechanical response in the nonlinear phase after the maximal force F c has been reached, as well as on the speed of crack propagation.

2.7e

-3

[mm] 7.5

-3 7.5 -3 ∆u (d<0.4)) [mm] 1.0e -4 1.0e -4 ∆u (d≥0.4)) [mm] 1.0e -5 -1.0e -6 1.0e -5 -1.0e -6 h [mm] 2.0e -3 2.0e -3
Table 3.2 M I and M II problem parameters. In Table 3.2 a coarse rst step ∆u (d< d)) = 1.0e -4 is chosen until the system reaches a threshold of damage state d (here d = 0.4 ). Then, the load displacement step is changed to ∆u (d≥0.4)) (here ∆u (d≥0.4)) = 1.0e -5 or 1.0e -6 ). In the present case of nonlinear solving procedure, the load step has been observed to yield instabilities with divergences of the Newton-Raphson for large loading steps and especially when d locally reaches 1. In the next numerical test, the ne load displacement step will be xed for ∆u (d≥0.4)) = 1.0e -6 .

Comparison between the linear and nonlinear formulations

The aim of this section is to qualitatively compare the responses obtained with the linear (Lin) and nonlinear (NL) phase-eld formulations. For this rst comparative study, we use the same denitions for problems M I and M II tests as previously. The crack propagation and the me- chanical response are observed. The idea of the test is to change synthetically the g c parameter controlling the crack apparition and to estimate the error between the classical linear mechanical problem and this non-linear formulation. Then, a parameter is dened as :

r gc = g c ḡc , (3.84)
where ḡc is the Grith energy parameter dened in Table 3.1 and commonly used in the literature as concrete property.

We rst compare the critical force response (see Fig. 3.8). We can observe that for dierent values of r gc , the dierence of the maximal force value between Lin and NL formulation is important when increasing the critical displacement u c controlled by an increase of the g c pa-

rameter. Physically, we can interpret this as follows: when the crack nucleates for a prescribed displacement, the linear formulation underestimates the force response and overestimates the critical displacement as compared to the geometrical non-linear formulation. More specically, by using the two comparison parameters (related to the critical force F c or critical displacement u c ), we dene the error between the linear and non-linear responses as:

err(X ) = 100 × 1 - (X) Lin (X) NL . (3.85)
Based on the previous tests, we can note that a signicant dierence between the Lin and NL formulation is achieved:

For the M I test :

The error err(F c ) reaches ≈ 16% when the crack starts to propagate around u c < 6%L

(r gc = 100).

The error on the critical displacement is stable and around err(u c ) = 3 -4%. For the M II test :

The error err(F c ) reaches ≈ 100% when the crack starts to propagate around u c = 12%L (r gc = 100)

The error on the critical displacement is signically larger than for tension (always above 6

These errors have a consequence on the crack initiation, as observed in Fig. 3.11 for the M I test and in Fig. 3.12 for the M II test. We observe that the crack propagation in the NL context appears before the classical Lin formulation.

Based on these rst results we conclude that even considering moderate amplitudes of loads, taking into account the geometrical nonlinear behavior has a signicant inuence on both the overall response of the sample as well as on the crack path.

3-point bending test

The three-point bending (3-PB) test is now investigated. The boundary conditions allow to reach the nonlinear regime quickly. The interest is therefore to observe with this formulation in large displacement how the structure behaves as compared to the small strain formulation, without changing the damage model. The geometry of the system is dened in Fig. 3.13 and the parameters are dened in Table 3.3.

Figure 3.13 3-PB geometry and boundary conditions.

In the 3-PB test, the stress concentration at the point of application of the u d displacement induces premature damage at this point. In order to avoid this, the evolution of the history function is deliberately locked at zero, over the entire (e = 0.2) top thickness of the sample.

As in the previous tests the mechanical responses are summarized in Fig. 3.14. Linear and non-linear formulations are compared directly with the elasticity without damage (∀t ∈ [0, t n [ d = 0). We notice once again that there is a signicant dierence between the Lin and NL responses. However, it is shown here that the Lin model overestimates the mechanical response as compared to the NL model.

Parameter [Unit]

Value

λ [GPa]
12.00

µ [GPa]
8.00

g c [kN/mm] 2.5e
-4

[mm] 3.0e

-3 ∆u (d<0.40)) [mm] 1.0e -4 ∆u (d≥0.40)) [mm] 1.0e -6 h [mm] 1.0e -3
Table 3.3 3-PB problem parameters.

These dierent tests (M I , M II , and 3-PB) allow to highlight, both quantitatively and qualitatively, the dierences observed on the maximum force response and crack initiation between the models derived at small and large strains. It is noted that signicant, problem-dependent dierences are observed.

Hydric shrinkage of a homogeneous sample

In this numerical test, we consider a cylindrical sample as described in Figs. 3.17 and 3.18. The parameter h of retraction is related to the hydric strain through:

E h = α h I , (3.86)
with α h = 1 2 (h 2 -1). The hydric shrinkage is modeled as a compressive strain that increases linearly with time, according to :

α h (t) = α h0 t , ∀t ∈ [0, 1] , (3.87)
where the strain coecient α h0 is arbitrarily set to (0, 0.3).

The parameters for the simulations are dened in Table 3.4). The geometries for the 2D and 3D cases are as follows: in 2D, the sample is dened in an initial circular domain of r 0 and the shrinkage parameter α h acts on plan. In 3D we use the same radius and a height h 0 = 0.1 (see Fig. 3.18), α h aects also the height.

Parameter [Unit] Value E [MPa] 1.00 ν [] 0.30 r 0 [mm] 1.0 h [mm] 4.0e -2
Table 3.4 Shrinkage problem parameters.

Note that no cracking parameters are used because we study in this part only the hydric shrinkage.

The evolution of the radius with respect to the retraction coecient α h is shown in Fig. 3.16, for both the small strain and nite strain formulations. As expected, we can note from Fig. 3.16 that the two models signicantly dier from one another for large values of the retraction coecient.

From a qualitative standpoint, the dierence becomes noticeable at α h = 30%, as seen in Fig. 3.17 for 2D plane strain and in Fig. 3.18 fo 3D, where contours of deformed congurations are displayed.

Figure 3.17 Dierent shrinkage radius in 2D for initial (black) and after drying process with linear (blue) and non-linear (red) formulation.

Figure 3.18 Dierent shrinkage sample and the projection of their radius in 3D for initial (black)

and after drying process with linear (blue) and non-linear (red) formulation.

The purpose of this test is to evaluate the evolution of the radius with the function of the water retraction parameter α h and then identify it experimentally.

Hydric shrinkage induced crack propagation

In this next example, we apply the above framework including nite strains, shrinkage strains and crack propagation. We consider a 2D circular domain as described in Fig. 3.19 (a)) containing a square sti inclusion in its center. The radius is r 0 = 50 mm and the length of the square is L i = 30 mm. The geometry was inspired from an experimental result of clay desiccation with such a geometry provided in [Barnier, 2015]. An illustration of the cracks induced in such test are provided in Fig. 3.20). The boundary conditions are as follows: the displacements on the external boundary ∂Ω (see Fig.3.19 (a) ) are set to zero in the beginning of the simulation and the boundary of the square inclusion are stress-free. When cracks reach the external boundary, the boundary conditions are changed to stress free conditions. This choice is a priori opposed to experience. It is rather the contrary for an ideal drying experience, i.e. without friction, the displacement would be blocked around the rigid square and free edges at the outer contour. Thus the clay would want to shrink, but the square prevents it from moving. Finally, the conditions adopted here are closer to a sticking of the sample on the edges of the cup. However, the experiment is probably imperfect and would contain surface forces (below and/or around) that would change the loading conditions. The sample would contract but outwardly, possibly due to increased friction (such as contact surfaces) from the middle, which would cause a crack near the centre. These conditions, which we are going to use all the same, represent an advantage in terms of simplicity of numerical implementation.

Figure 3.19 (a) Geometry and boundaries conditions, (b) Mesh, (c) experimental result [Barnier, 2015]. We use a linear evolution of α h (t) with respect to time in Eq. (3.86) in the form

Parameter [Unit] Value E [MPa] 1.00 ν [] 0.30 g c [N/m] 1.15 [mm] 2 r 0 [mm] 50.0 L i [mm] 30.0 ∆t(d < 0.1)) [] 1.0e -2 ∆t (d≥0.1)) [] 1.0e -4 h [mm] 0.8
α h (t) = A.t , (3.88)
with A a constant A = -0.1, for a range time of t ∈ [0, ], using rst, a time step ∆t = 10 -2 s and then ∆t = 10 -4 s when d(x) > 0.1. All parameters of the simulation are listed in Table 3.5

Even though this example is used for illustration only of the theoretical model framework and does not intend to reproduce the experiments, we can appreciate that qualitatively, the cracks start from the boundaries of the inclusion and then propagate to the external boundary. After the propagation, the remaining parts continue to shrink as in the experiments.

Conclusion

In this chapter, we have presented a phase-eld formulation for nite strains that accounts for compressive shrinkage strains. Borrowing most of the formulation from the literature, we rst focused on parametric analyses, with the aim of characterizing the sensitivity to model parameters. In particular, we have investigated the eects of geometrical nonlinearities, even at moderate loading, and quantied the eects of the parameters on the convergence and stability for the macroscopic response. We then introduced a simple formulation where shrinkage is taken into account. A simple application example has been presented involving a heterogeneous sample, where both shrinkage and cracking were involved. In Chapter 5, the proposed model will be used for comparisons with the experimental results presented in the next chapter.

Introduction

This chapter presents an experimental study of the desiccation processes of a model heterogenous clayey medium. First, the clay material properties, the setup and the experimental conditions are presented. This allows observing and quantifying the shrinkage, the initiation and the prop-

Experimental setup 4.2.1 Material

The material used in this experimental study is Romainville clay, found in the East Paris Basin.

Its ability to shrink, swell and crack causes a lot of damage to buildings, especially in periods of severe drought. That is why several study campaigns have been conducted [START_REF] Audiguier | Caractérisation au laboratoire de la sensibilité au retrait-gonement des sols argileux[END_REF], Zemenu et al., 2009], to understand and analyze these phenomena. Geologically this clay is a composition of dierent sediments: Illite and Scmectite, carbonate, quartz and feldspath.

The physical properties of this clay are referenced in Tab. 4.1.

Soil properties Values

Density of solid phase [Mg.m

-3 ] 4.1 Physical properties of Romainville clay [START_REF] Tang | Experimental characterization of shrinkage and desiccation cracking in thin clay layer[END_REF].

(1) Liquid Limit is the water content at which soil changes from a plastic to a liquid state, the indication is related to the mass proportion of water in the mixture.

(2) Plastic Limit is the water content at the change from a plastic to a semisolid state.

(3) The plasticity index of a soil is the dierence between its liquid and plastic limits.

4.2. Experimental setup

Experimental method

For the preparation of the experiment, the clay paste was conditioned as in the studies of [START_REF] Tang | Experimental characterization of shrinkage and desiccation cracking in thin clay layer[END_REF], Tang et al., 2012] and many others. In this work the same clay was used, and experimental investigations were carried out on cracking under hydro-mechanical eect, cyclic loading, for dierent temperatures, and other physical phenomena. For the sample preparation :

the clay was taken directly from the site (Romainville, Paris-Est, France) by block;

the material was cut into small 2 cm pieces and immersed in distilled water to liquify it for a period of t = 24 hours; then the whole set was passed through a s = 2 mm sieve, allowing us to lter out the larger particles;

the resulting mixed sludge contains a water content of about wc = 170% at this stage, and was poured into D = 116 mm diameter petri dishes. Rigid inclusions (cylinders with circular cross section made of PMMA) were randomly placed in number (n = 0, 1, 3 or 6)

as illustrated in Fig. 4.1. Note that the Petri dishes in which the clay is poured, is covered with a Teon lm in order to limit the eects of friction on the surface (both on the lateral walls and the bottom) as much as possible. A diagram of the preparation is shown in the last step is the removal of air bubbles from the clay slurry: the samples are placed under vacuum for a period of 2 hours and left covered for sedimentation for 72 hours; before starting the analysis, the clay mud has evacuated the supernatant water lm on the surface; the water content at t= 0 s of the test is wr ∼ 115% and the thickness is H = 8 mm.

The clay sample is ready to be analyzed, and is placed in a wooden box of dimension V = 50×50×50 cm, preserves the sample from rapid variations of thermal or humidity variations or the surrounding the ambiant air. This box is however not tight, and its internal hygrothermal,conditions evolve with current laboratory conditions. Several measuring tools are placed in this box, namely: a scale, a thermometer and an ambient humidity probe, giving access to the water loss of the sample (and thus after processing of the full data the water content), the temperature and the humidity of the environment throughout the drying test.

After 72 hours, the sample has suciently stiened and has been painted with black speckles on its surface, allowing to follow the retraction kinematics and to access the local deformations using a digital camera (Canon EOS TTL camera, equipped with a 18-55mm zoom lens, providing 8bit color images with 5184x3456 pixels). The image acquisition is done at time intervals of ∆t = 10 min for dierent samples with dierent numbers of inclusions, as schematized in Fig. These images are then analyzed using the digital image correlation technique (DIC) [START_REF] Bornert | Mesures de champs et identication en mécanique des solides, chapter corrélation d'images[END_REF], Wang et al., 2018], and provide access to the 2D deformation eld at the upper surface of the sample. The experiment thus provides many information during the drying process such as: temperature, humidity, water volume fraction, crack evolutions, and also the deformation eld. It is noted that many simulation models in literature only provide crack pattern images.

Here the DIC provides a new, more accurate and relevant level of understanding of drying kinematics and crack dynamics, which will help enrich our model.

Results

The collection of images, is then analyzed by DIC, using an image correlation software developed (CMV) at Ecole des ponts ParisTech, France. The complete description of the calculation of the local deformations is described in the work of [START_REF] Wang | Nucleation and propagation mechanisms of soil desiccation cracks[END_REF] using the same software. The processed images provide the local displacement eld and nally a map of the local deformations including two quantities of interest that will be exploited: the spherical (or hydrostatic) strain S part and the deviatoric part D dened as S = 1 2 tr( ) , and

D = 2 3 (¯ 2 -¯ 1 ) , (4.1)
¯ 1 and ¯ 2 are the eigenvalues of the 2D Green-Lagrange strain tensor, and tr( ) = 1 + 2 measures the variation of surface. It is emphasized that s is dierent from a volume variation which cannot be measured locally by purely surface investigations. These quantities are evaluated for a gauge length limited by the typical length scale of the speckle painting and the (adapted) optical magnication of the camera. With about 2700 pixels along the diameter of the Petri dish, there are about 50 independent local evaluations of the 2D strain tensor along such a distance, so that the spatial resolution is of the order of 2mm. Accuracy of the measurement is limited by image quality and noise and other artifacts not detailed here, but is better than 1%, and thus sucient to quantify the heterogeneity of the investigated strain eld.

Results

Experimental conditions

The global study has been performed on dierent samples :

S 0 : a clay sample without inclusion, S 1 , S 3 , S 6 : clay samples containing n = 1, 3, 6 randomly-distributed rigid inclusions, S 1c : sample with a unique inclusion placed at the center of the circular geometry.

The rst results, shown in Fig. 4.4 rst indicate the conditions of the experiment, with an ambient relative humidity of RH = 48% ± 6, a temperature of T = 20 ± 2 o C, and the water content (ratio of weight of water and clay) of wc = 115-10%. Then the relatively small variations of these parameters from one sample to another show that the experiment was carried out under similar measured conditions and shows the repeatability of the environment for a drying study.

Even if the conditions might be slightly dierent form one sample to the other, it turns out the drying rate is very slow (tests last about one week) with respect to the typical time of the moisture transfer within the clay matrix, so that the water content in the clay material can be considered uniform and is directly given by the overall mass loss quantied by the scale. The latter is almost linear with time, as observed in the last graph of Fig. 4.4. Moreover, as viscous eects are likely to be negligeable in such experiments, the current loading state of a sample is directly provided by this water content which apperas thus as the main loading parameter, to be used in particular for comparisons with numerical simulations.

Crack process

Figs. 4.6,4.8,4.7,4.9 and 4.10 show the states of the material at dierent times t = 4000, 6000, 8000 min each representing an initial state without cracks, the rst propagations, and the nal state.

The crack paths obtained on each specimen show some similar features of the propagation pattern. The cracks initiate near the rigid inclusions and propagate towards to the outer edge of the specimen. Once all cracks have propagated, the material shrinks over its entire range, dividing the material into small independent subassemblies. It should be noted that the crack paths over the entire test set consists of 3 or 4 branches with angles betwen these branches θ c such that: π/2 θ c π/3.

Strain analysis

In this part the deformation maps from the DIC technique are presented for the three steps in the cracking process: shrinkage without cracking (t =4000), crack initiation (t =6000), then propagation (t =8000). S B is a ring surface dened as near-edge surface ring , S + is another ring surface, the closest to S B , S i is the rest of the surface introduced as an interior surface.

S 0 test

The spherical and deviatoric deformation show a quite remarkable radial deformation pattern (see Fig. 4.6):

First, at t = 4000 min, the intensity of the spherical deformation is rst negative (as expected) on S i , of low intensity S ∼ -2%, homogeneous on the surface except on the surface close to the edge S B ∪ S + where another regime settles. We observe a surface in a positive deformation crown (S + ) then another negative deformation crown (S B ) more intense than in the center of the sample. Concerning the deviatoric strain, it appears to be almost zero in the central area S i , indicating a purely isotropic surface deformation as expected from a homogeneous shrinkage. However the deviatoric strain exhibits a small value D ∼ 2% in the surrounding area S B ∪S + where the surface strain was heterogeneous.

It appears thus that the sample shrinks almost homogeneously and isotropically (in 2D) in its central part, as expected, but that it also undergoes some more complex deformation at its outer edge, due to some boundary eect. A possible explanation might be associtated with some adherence of the sample on the walls of the dish.

Secondly at t = 6000 min, a wall debonding appears at the top left of the sample. The spherical and deviatoric strain maps are modied in consequence. On the rst crown (S i ) the spherical strain remains at S ∼ -2%, while the second S + crown splits into two distinct parts: a rst part which is not aected by the wall detachment as seen at t = t 1 , then the second part close to the detachment where there is no longer an extension but now a strong compression S ∼ -14%. These observations conrm the adhesion of the sample to the external walls. The last crown on the periphery of the border, which remains negative, is of strong intensity. Concerning the deviatoric strain, the obverall features observed for the previous step are still present. But the debonding induces a more intense deviatoric strain in its vicinity. In addition, some heterogeneities in the central part S i are now much more visible.

In the last step at t = 8000 min, the sample is fully detached form the dish walls and the strain eld exhibits again an essentially axisymetric distribution. The contrasts in strains become more accentuated, with two distinct zones: the external zone (S B ) of strong intensity S ∼ -40% and S ∼ 15%, and the central zone S + ∪ S i with a less intense but nevertheless important shrinkage of S ∼ 15% associated with a moderate deviatoric strain D ∼ 5%. Some heterogeneities are also observed in areas where they had already been seen in the earlier steps.

S 1c test

The S 1c test is comparable to the test without inclusion (S 0 test), with additional features in the vicinity of the inclusion, see Fig. 4.7:

At step t = t 1 , an additional ring surface surrounds the inclusion (S B ). This area has a strong retraction S = -10%, equal to the deformation on the outer edge. Moving away from the crown close to the inclusion, we have, as in the previous case (S 0 test), a zone S i of negative strain of weak intensity S ∼ -2%. Then a last zone (S + ), in weak extension ( D ∼ +4%). For the deviatoric part, the intensity is strong ( D ∼ 5%) near the edges S B , then weak elsewhere S + ∪ S i (with D ∼ 0%). These observations show that the perturbations induced by the presence of the inclusion are concentrated in area with a diameter of the order of three times the diameter of the inclusion. Because of the central position of the inclusion, these perturbation do no interact with the phenomena observed at the periphery of the sample. There is an intermediate zone which behaves almost as in the S 0 test.

At t = t 2 a debounding of the clay on the wall appears at the bottom of the sample. The strain maps are modied as in the case without inclusion. The rst crown S B linked to the inclusion is still in maximum shrinkage as well as the outer edge ( S ∼ -14%). Then the surface S i keeps a weak retraction ( S ∼ -2%). Then the last crown S + splits into two distinct parts: a rst part which is not aected by the wall detachment, identically to the previous observation (t = 4000 min). For the deviatoric part, the overall intensity increases with a greater concentration on the outer edges S B and the inclusion ( D ∼ 10%) for a relatively low intensity ( D < 4%) on the remaining surfaces S + ∪ S i . The debonding at the lower part induces also a local increase of the deviatoric strain.

Finally, at t = t 3 , the clay separates from the edge over 40% of the outer edge. In contrast to the previous case (S 0 test), the upper half of the crown S + remains in extension ( S ∼ +5%). For the other part, as in the previous commentary, there is a strong shrinking zone ( S ∼ -20%) on the borders S B and a lower one ( S ∼ -5%) on the last remaining area S i . for the deviatoric part, intense strains are again observed in the edges S B (outer part and near the inclusion). It is almost null in areas far away from the borders, the inclusion, especially on the non-debonded half part of the sample.

Note also that the very intens value of the deviatoric strains at the borders of the cracks is a post-processing artefact and should not be interpreted as the local deformation of the clay. Indeed the displacement jump over the crack is integrated in the local evaluation of the strain which is based on some nite dierences (see [START_REF] Allais | Experimental characterization of the local strain eld in a heterogeneous elastoplastic material[END_REF] for the details about the calculation of the strain).

S 1 test

In this test, a single rigid inclusion is inserted within the sample, but is here o-centered to break the symmetry of revolution.

At t = t 1 , as in the previous cases, the high intensity areas are located on the edges S B of the sample exhibiting a large contraction, of intensity S ∼ -10%. One may also notice that the eld surrounding the inclusion is also more heterogeneous than in the S 1 c sample.

In this test the crown S + is not apparent, the remaining surface S i being a mixture of low intensity positive and negative surface stress | S | < 5%. For the deviatoric part, the intensity is relatively, as in the previous samples, strong D ∼ 5%, on the outer areas S B , and then weaker in the inner area S + ∪ S i , but with some noticeable heterogeneities.

At t = t 2 , the intensities on the outer edges of S B increase with S ∼ -14%. The S + zone does not seem to develop in this case. The inner zone S i remains of relatively weak intensity S ∼ 2% with nevertheless some scattered zones with a stronger intensity up to S ∼ +10%. Regarding the deviatoric part, the intensity continues to increase on the edges S B and remains small on the interior surface S + ∪ S i .

At the last step t = t 3 , the specimen maintains and conrms a trend of strong retraction S ∼ 20% on the outer edges S B . The inner zone S i has extensional strains S ∼ [5%, 10%], larger than in the previous S 0 and S 1c tests but remaining relatively low. For the deviatoric part, the trend seems to be conrmed with an intensication D ∼ [7%, 15%] on the outer edges S B and weaker values in the inner areas S i .

In this sample with an excentred inclusion, debonding at the outer edge was limited. Adhesion induces very strong contractive strains in areas very close to these edges, together with extensive strains of lower intensity in larger zones and far away from the edges. So that the average surface change is zero, as required from the absence of debonding. Some uctuations are observed in these area in extension but are not associated with strong deviatoric strains. The presence of the inclusion in the lower part generates rst local deviatoric strains and then cracks. This allows the surrounding clay to shrink, rst in the near-neighborhood of the inclusion and then in the larger areas around the cracks. A general observation is that the overall features of the strain eld around the inclusions in S 3 are similar as those observed in S 1 and S 1 c. This is also true for the crack geometries:cracks emitted from an inclusion do not interact with the others. This is no longer true for the S 6 sample in which positive S strains are no longer observed between the closer inclusions. In addition cracks have dierent features: some are generated in areas with positive S , and one clearly connects two inclusions (the lower left ones).

Discussion

Main mechanisms

A clear pattern of deformation and cracks is observed for all the tests carried out. At the beginning of the desiccation process, the material generates a high-intensity shrinkage deformation close to all edges and dened by the surface S B . The rest of the surface is shared by (i) a low intensity, and overall shrinking portion, (ii) or possibly an extensional portion dened as the higher intensity S + surface. This extended area appears relatively far from the edges, and may partially disappear as the material debounds from the wall, or as a crack passes through it. Eventually, a large number of cracks have been generated, and the sample globally shrinks, facilitating and enlarging the cracks.

Stochastic aspects

The deformation and cracking processes share some common deterministic features as discussed above. But their are also characterised by some strong stochastic aspects, regarding in particular the detailed geometry of the cracks. This is in particular illustrated in Fig. We provide here some observations which might be taken into account for future stochastic modeling studies of the shrinkage process:

1. The ground friction has been limited by the application of a non-stick Teon layer. However a quantication of this surface interaction remains to be identied. Numerical studies on an experimental basis [START_REF] Sima | Numerical simulation of desiccation cracking in a thin clay layer using 3d discrete element modeling[END_REF], Amarasiri et al., 2011] have studied the inuence of this interaction on the cracking process.

2. The detachment at the outer walls inuences the overall kinematics of retraction and has been found in many similar experimental works. The construction of a model specically accounting for this phenomenon could be proposed and inferred.

3. This debonding could be modeled through interface models, possible accounting for the variability. Interface models using the phase eld modeling has been carried out in [Verhoosel and de Borst, 2013[START_REF] Nguyen | A phase-eld method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography[END_REF].

Introduction

In this chapter, the numerical modeling tools developed in chapter 3 are used to simulate crack initiation and propagation during hydric shrinkage of clay samples. The experimental results presented in chapter 4 are then used to identify and carry out comparisons with the numerical simulations. This work and results has been presented in [START_REF] Hun | Computational modeling of crack propagation in a heterogeneous medium under drying conditions[END_REF], Hun et al., 2019a] and adapted from the scientic paper [START_REF] Hun | Crack propagation in a heterogeneous clay material under drying conditions, experimental study and numerical simulation[END_REF].

Description of the problem

We consider the numerical model schematically described in Fig. 5.1. For a given test (that is, for one sample and a given number and distribution of obstacles), a digital replica is dened using the experimental data described in Chapter 4. Based on measurements, the initial radius of the cylindrical cup is set to r 0 = 58 mm, the radius of inclusions to r inc = 8 mm, and the height to h 0 = 8 mm. The clay material is contained within the domain Ω (see Fig. 5.1 (a) for an example with 3 inclusions).

In order to account for the damage mechanisms observed experimentally, boundary conditions are dened as follows:

Zero-displacement Dirichlet boundary conditions are initially applied on ∂Ω, while stressfree conditions are considered at the boundaries of the inclusions, denoted collectively as ∂Ω . Shrinkage is then induced by increasing the eigenstrain E h = α h I (see Eq. (3.69) 3 ), where the parameter α h is to be dened. The problem (3.69) is then solved at each time step, and fracture is observed if localized traction area are created within the sample.

When the condition d = 1 is met on a subset ∂Ω deb of ∂Ω (note that ∂Ω deb is not connected in general), boundary conditions on ∂Ω deb are switched to free tractions, with the aim of modeling the sudden debounding observed during the experiments.

Concerning the bottom of the sample, denoted by ∂Ω z=0 , the boundary conditions are vertically xed to 0. The mechanical parameters used in the model described in section 3.3.3 are extracted from the literature. The Young's modulus for the clay material is taken as E = 1 MPa, according to the study in [START_REF] Vo | Numerical modelling of desiccation cracking of clayey soil using a cohesive fracture method[END_REF], where the simulation of crack propagation in the 2D plane with a cohesive zone model on clayey materials was conducted. This modulus was dened experimentally in [El Mountassir et al., 2014] where a variable Young's modulus was dened as a function of the compaction of the material. In our clayey material, we assume to be in a similar situation when the material solidies according to the decrease of the water content. In [START_REF] Vo | Numerical modelling of desiccation cracking of clayey soil using a cohesive fracture method[END_REF], this modulus was averaged over the range of compaction. This hypothesis will be retained for the numerical study. In other works, [START_REF] Cajuhi | Phase-eld modeling of fracture in variably saturated porous media[END_REF] proposed a variation of this modulus according to the water content in the experimental work of [Lakshmikantha, 2009] on the Barcelona soil, and in [START_REF] Peron | Desiccation cracking of soils[END_REF] where Bioley clayey silt was investigated. The Poisson's ratio will be chosen as ν = 0.3 [START_REF] Vo | Numerical modelling of desiccation cracking of clayey soil using a cohesive fracture method[END_REF].

Fracture parameters

The damage parameters used in the following simulations require the estimation of the toughness g c and of the regularization length , as this parameter is here interpreted as a material parameter (see a discussion in [START_REF] Nguyen | On the choice of parameters in the phase eld method for simulating crack initiation with experimental validation[END_REF]). In the literature, many authors have studied the impact of water content on crack resistance. [START_REF] Kodikara | Desiccation cracking in clayey soils: mechanisms and modelling[END_REF] has provided a comprehensive review of environmental factors in the cracking process through multiple experiments.

More specically, numerous studies have allowed identifying parameters of linear mechanical fracture models through experiments (see [START_REF] Wang | Experimental study on fracture toughness and tensile strength of a clay[END_REF], Lakshmikantha et al., 2009, Peron et al., 2009]). The parameters used are often the stress intensity factor K I or the tensile limit stress σ t , which can then be theoretically related to the pair of parameters (g c , ) of the phase eld method (see [START_REF] Nguyen | On the choice of parameters in the phase eld method for simulating crack initiation with experimental validation[END_REF]). Following these works, a numerical study on clay desiccation using the phase eld method was reported in [START_REF] Cajuhi | Phase-eld modeling of fracture in variably saturated porous media[END_REF]. The parameters used are a cracking energy depending on the water content (g c = 1.12(wr = 0%) -1.71(wr = 100%) N/m) and a constant regularization length = 2e -3 m. For our study, we will consider the toughness as a constant, with g c = 1.15 N/m, and set the regularization length to = 2 mm.

Hydric model

In this section, we describe the procedure to identify the function α h (t) of the coupled hydric shrinkage model (see Eq. (3.69)) from the experiments. The experimental test S 0 described in section 4.6.1 has been used, where a homogeneous circular sample is considered. In this test, the sample shrinks without cracks and remains circular during the drying process. Denoting the experimental radius by R exp (t), it is possible to record the evolution of this radius as a function of time. This evolution is depicted in Fig. 4.6 and the radius is reported in Fig. 5.2 (Left).

The experimental evolution of the radius (red dots) has been tter by the R T function dened as: R T (t) = A(1 + tanh(B(t -φ)) + R ∞ , (5.1) where A is the amplitude dened by

A = 1 2 (R ∞ -R 0 ) , (5.2)
with R 0 and R ∞ denote the initial and nal radius, B denes the transition time between radius R 0 and R ∞ according to Then, the parameter α h was identied using the numerical application in section 3.5.5, and is directly related to the radius retraction measurement. Thus we obtain the time function of α h , illustrated in Fig. 5.2 (right). In the next simulations, the identied function α h (t) is used, even the ones involving cracks, keeping in mind and being aware that this is a strong simplication of the model. The others parameters are summarized in (with 1TB of RAM and 2×24 cores rated at 2.7 GHz) ranged from 0.7 to 3 seconds for the 2D case, and from 55 to 150 seconds in the 3D case. Total computational times based on 1000 time steps were 2.5 hours and 42 hours for each 2D or 3D case, respectively.

B = 1 ∆T , ( 5

Simulation results

Crack patterns

In order to obtain a meaningful comparison between the experimental tests and the numerical simulations, we dene the following characteristic times: t = t 1 is associated with damage initiation within the (solid) sample; t = t 2 corresponds to the time at which the rst crack starts propagating; t = t 3 is the time at which all cracks have propagated and do not evolve anymore.

The above times were estimated for all experiments, and are given by: t 1 = 4000 min, t 2 = 6000 min, and t 3 = 8000 min. In what follows, we compare the crack patterns observed experimentally with the proles predicted by the numerical model.

S 1C conguration

First, we consider the specimen S 1c , which contains a single, centered inclusion; see Fig. 5.3. It is observed that the specimen exhibits 3 crack branches in the physical experiment. The 2D simulation also gives rise to 3 branches, whereas only 2 primary branches (that subsequently split into 2 secondary branches each) are obtained in the 3D case. However, it should be noted that the solution is not unique, and that the propagation only depends on the numerical perturbations related to mesh construction in the present case. The discrepancy between the 2D and 3D cases may also be explained by the fact that the 3D simulation includes through-thickness shrinkage that the 2D model does not.

S 1 conguration

We now consider a second example (specimen S 1 ) where a non-centered inclusion breaks the symmetry of revolution; see Fig. Based on these results, we note that the crack patterns predicted by the 2D and 3D computational models present good qualitative agreement with the experimental observations. In particular, it is noticed that crack propagation always initiates in the vicinity of an inclusion and remains mostly radial (until the boundary of the sample is reached, at which point debounding occurs). Coalescence between cracks having initiated at adjacent inclusions is also observed in the simulations and on certain experiments (as S 6 ).

Spherical (hydrostatic) strain distribution

S 1C conguration

In the case of the centered inclusion (see Fig. 5.7), we compare the distribution of the hydrostatic eld s between the experiment and the 2D/3D simulations. The same conclusions can be drawn for the conguration with a single centered inclusion; see Fig. 5.12.

Figure 5.12 Deviatoric strain map in S 1 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right). Here, a better agreement on the strain elds is obtained. In the experiments, areas of higher intensity around the inclusions and the outer edges appear during the damage initiation. These eects are well captured in the simulation, but in the vicinity of the inclusions only. In the rest of the sample, the deformation values remain quite low and in agreement with the experiment overall. Finally, the locations of inclusions seem to play a more predominant role in the simulations where interaction eects may give rise to privileged crack paths.

Discussion

In the above comparisons, we have observed some encouraging qualitative agreement on some points between experiments and simulations, while some others remaining to be claried. These points can be summarized as follows.

1. The zero-displacement boundary condition applied at the outer edge is a rst approximation of the frictional forces that prevent the clay matrix from retracting freely, hence creating the conditions for cracking. However, it has been observed that the clay matrix slightly lifts o the bottom at the edges, and that the amplitude of this phenomenon increases with time. This results in a soil-clay interaction zone oset. The modeling of frictional forces is therefore a natural extension to the present model. To illustrate this point, results are presented in Fig. 5.15 where the sample is not retained on its outer edges, but in a zone of small thickness close to the edges and in the clay. In this case, the edge eects observed experimentally are qualitatively recovered.

2. A localized, compressive strain zone is experimentally observed close to the outer edge and at the interface between the inclusions and the matrix. This complex interface phenomenon requires further studies and is presently not captured by our model. 3. The parameters which have been used in the simulations have been chosen as constant.

However, some experimental results available in the literature show that these parameters may be spatially varying and depend on the state of water content. In this context, a more advanced experimental campaign involving variable water content could enhance the predictive capability of the computational model.

4.

A large discrepancy is also observed in terms of crack propagation velocity (quantied by crack extension in a quasi-static regime). While crack propagation is typically observed over long time ranges in the experiments (about 3000 to 4000 min), crack propagation is about 10 to 20 times faster in the corresponding simulations.

5. Due to time constraints, it was not possible to include and model stochastic aspects such as those developed in Chapter 2. It is however reasonable to assume that random subscale details within the clay induce substantial variability in the elastic and fracture properties.

Other sources of uncertainties related to the experiments have also been listed at the end of Chapter 4. One priority for future works will be to include all these uncertainties within the modeling framework.

Conclusion

In this chapter, we have compared results obtained through the drying experiments detailed in Chapter 4 with predictions delivered by 2D and 3D numerical models. The experimental results were notably used to identify an appropriate form for the hydric shrinkage strain. An extensive simulation campaign was then carried out to investigate the impact of inclusions on crack propagation, for various congurations.

From a qualitative standpoint, predictions were shown to be in reasonable agreement with the experimental results. The strain elds and the hydrostatic strain are mostly negative, and localize on the edges of the inclusions. As expected, localized strain elds are also observed in the vicinity of crack paths.

Many issues remain, however, to obtain better quantitative agreement. The present model lacks in representing some features observed during the experiments, such as local compressive strain at the matrix-inclusion interface and near the external boundary. This strain induces a large discrepancy in the numerical values of local strains and crack paths.

We have suggested several improvements to enhance the present models, including a deeper analysis of the cohesion at the boundaries and at the matrix-inclusion interface, the introduction of friction at the bottom surface, a more realistic phase change model for the clay (which passes from liquid to solid state), and the introduction of uncertainties within the model. identication of the hydric strain function based on the variation of the sample geometry during the shrinkage has been carried out. We then addressed validation with dierent congurations of obstacles, using both 2D and 3D numerical models. This validation step constitutes a rst step towards the construction of such models, and while some qualitative agreement was found overall, some major discrepancies were also identied. This suggests that the present model could be enhanced in many ways and possible strategies to improve predictions, reported as perspectives for future work below, were proposed.

Perspectives

There are many perspectives for this work, which are listed below.

The stochastic model presented in Chapter 2 could be extended by generating random elds of damage parameters. In the present work, we have used a random elastic medium, but the eective stochastic medium was limited to the elastic aspects and the eective fracture was considered as constant. Extension to random elds of fracture models constitutes a direct perspective for this work.

This stochastic fracture model could be applied to a large class of other materials such as porous media, lattice models or polycrystalline materials. A study using this identication model applied to other geometries or stresses would be another perspective.

As mentioned above, the constructed phase eld model for the drying-induced model still Finally, a complete description of uncertainties of the process and its introduction of the model will be an important improvement. A rst class of uncertainties are related to the experimental conditions (decohesion of the boundaries, friction on the bottom, etc).

Another class of uncertainties are related to the material, where the heterogeneities at the microscale induce random eects during the drying and inuence the crack paths. This is obvious in homogeneous samples, where the crack patterns are highly random, and more deterministic when obstacles induce predetermined localization paths during the process.

An analysis on how these macro heterogeneities (obstacles) change the stochastic nature of the crack path will be a challenging and exciting perspective for this work.

n : the number of nodes in one element Ω E , p : the dimension of the space, x p (ξ) : the p-coordinates of x, ξ p : the p-coordinates of ξ dene on the isoparametric basis Ē, X p : the p-coordinates of X dene on the physical basis Ē, (i) p : the p-coordinates of the (i)-node, N i (ξ) the shape functions dene on isoparametric basis.

X
Jacobian matrix of the transformation T :

[J] ij (ξ) := ∂x j (ξ) ∂ξ i = x j ,ξ i , A.4) [J(ξ)] (pxp)

(
p=2 := x 1,ξ 1 x 2,ξ 1 x 1,ξ 2 x 2,ξ 2 (2x2) (2D case) p=3 :=  
x 1 (ξ) ,ξ 1 x 2 (ξ) ,ξ 1 x 3 (ξ) ,ξ 1 x 1 (ξ) ,ξ 2 x 2 (ξ) ,ξ 2 x 3 (ξ) ,ξ 2x 1 (ξ) ,ξ 3 x 2 (ξ) ,ξ 3 x 3 (ξ) ,ξ 3   (3x3)

(for 3D case), using the Eq. A.3 interpolation the jacobian can be written as :

[J] (A.6) it can be deduced that : ∇ x (f (x)) = [J] -1 ∇ ξ (f (x)), (A.7) in particular f can be interpolated using in the same way the Eq. A.3 :

f (x) = n i=1 N i (ξ)f (i) , (A.8)
Here f (i) is the nodal value of f function estimated on the i ∈ 1, ..., n node . computational evaluation of integral : The integral is estimated by the classical Gauss quadrature for example here :

Ω E f (x)dΩ E = Ω Ē f (
:=     N u1,1 (x) 0 N u2,1 (x) 0 N u3,1 (x) 0 N u1,2 (x) 0 N u2,2 (x) 0 N u3,2 (x) 0 0 N u1,1 (x) 0 N u2,1 (x) 0 N u3,1 (x) 0 N u1,2 (x) 0 N u2,2 (x) 0 N u3,2 (x)     {u i e },
or for 3D (T4 element) as :

{∇u} (3D) :=               N u1,1 0 0 N u2,1 0 0 N u3,1 0 0 N u4,1 0 0 N u1,2 0 0 N u2,2 0 0 N u3,2 0 0 N u4,2 0 0 N u1,3 0 0 N u2,3 0 0 N u3,3 0 0 N u4,3 0 0 0 N u1,1 0 0 N u2,1 0 0 N u3,1 0 0 N u4,1 0 0 N u1,2 0 0 N u2,2 0 0 N u3,2 0 0 N u4,2 0 0 N u1,3 0 0 N u2,3 0 0 N u3,3 0 0 N u4,3 0 0 0 N u1,1 0 0 N u2,1 0 0 N u3,1 0 0 N u4,1 0 0 N u1,2 0 0 N u2,2 0 0 N u3,2 0

Figure

  Figure 1.4 DEM Procedure.

  Figure 1.12 Level set functions (φ, ψ) dening a 3D crack.

F

  Figure 1.14 (a) sharp description of a cracked solid; (b), (c), (d): smeared description within the phase eld framework (damage proles d ∈ (0, 1) for dierent regularization lengths ).

Figure 1 .

 1 Figure 1.15 Multiple, complex cracks networks in a model concrete sample simulated by the phase eld[Nguyen et al., 2016a]).

Figure 1 .

 1 Figure 1.16 Inuence of the crack density function and of parameter .

Figure 1 .

 1 Figure 1.19 Inuence of the degradation function g on the material response U * → F * (U * ), where F * is the adimensional force response and U * denotes displacement.

Figure 2 . 2

 22 Figure 2.2 Boundary conditions applied to the domain and initial crack.

Figure 2 . 3

 23 Figure 2.3 Evolution of the damage eld x → d(x) for the microstructure shown in Fig. 2.2.

Figure 2 . 4

 24 Figure 2.4 Simulated displacement-force response for the microstructure shown in the left panel in Fig. 2.2. the stochasticity induced on the macroscopic response can be characterized by computing the horizontal force on the top edge (x 2 = L) of the samples, denoted by F . The mean and standard deviations for these quantities of interest, together with sample-based envelopes are shown in Fig. 2.5.

F

  (uD) F (uD) + σF (uD) Envelope over sample

Figure 2 . 5

 25 Figure 2.5 Envelope, mean and standard deviation for the displacement-force curve (left) and the crack path (right).

  Figure 2.7 One realization of the random eld { [ C iso KU BC (x)] F , x ∈ Ω} for R/R ∈ {0.2, 0.6, 0.8, 1, 2, 4} (from left to right). the following random eld {A BC (x), x ∈ Ω} is introduced [Guilleminot and Soize, 2012b, Tran et al., 2016]:

(Figure 2 . 8

 28 Figure 2.8 Graphs of the mean functions for the random elds {A KU BC (x), x ∈ Ω} (left) and {A SU BC (x), x ∈ Ω} (right) for R/R = 3.

Figure 2 . 9 Figure 2 .

 292 Figure 2.9 Graphs of the elds of standard deviation for the random elds {A KU BC (x), x ∈ Ω} (left) and {A SU BC (x), x ∈ Ω} (right) for R/R = 3.

Figure 2 .

 2 Figure 2.11 Graph of the correlation function along e (1) and e(2) , estimated from the simulated

Figure 2 .

 2 Figure 2.12 Graphs of the cost function for KUBC and SUBC (left), and envelopes of the macroscopic response for the reference microscale model and the mesoscopic formulations (right).

  g c = 1.255 × g c,m ≈ 3.375 × 10 -3 [kN.mm -1 ] (2.23) for SUBC. A comparison of the macroscopic responses obtained with the mesoscopic formulations (for KUBC and SUBC), parametrized with the identied values, and the reference computations

Figure 2 .

 2 Figure 2.13 Illustration of dierent crack path for the same microstructure: microscale description (left), KUBC-based description (middle), SUBC-based description (right).

Figure 2 .

 2 Figure 2.14 Comparison of crack paths obtained with the microscopic description (black dashed/solid lines) and the mesoscopic formulation. Left panel: case of KUBC. Right panel: case of SUBC.

Figure 2 .

 2 Figure 2.15 Plot of the correlation function τ → r(τ ; α) for dierent values of α.

Figure 2 .

 2 Figure 2.16 Kernel density estimates for the probability density function of the bulk (left) and shear (right) moduli, for the two types of boundary conditions KUBC (blue) and SUBC (red).

Figure 2 .

 2 Figure 2.17 Kernel density estimates for the joint probability density function of the bulk and shear moduli, for homogenization-based (left) and model-based (right) samples (for SUBC).

Figure 2 .

 2 Figure 2.18 Normalized correlation function of the bulk modulus random eld along e (1) (left) and e (2) (right), estimated from the multiscale data (red line) and the calibrated stochastic model (black line).

Figure 2 .

 2 Figure 2.19 Normalized correlation function of the shear modulus random eld along e (1) (left) and e (2) (right), estimated from the multiscale data (red line) and the calibrated stochastic model (black line).

Figure 2 .Figure 2 .

 22 Figure 2.20 Comparison of crack paths obtained with the microscopic description (black dashed/solid lines), the mesoscopic formulation where the elasticity is obtained from microstructural samples (blue lines), and the mesoscopic formulation involving the stochastic model for elasticity tensors (red lines). Left panel: case of KUBC. Right panel: case of SUBC.

  Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Finite strain kinematics and mechanical modeling . . . . . . . . . . . 3.2.1 Background in nite elasticity . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Phase-eld formulation for nite strains . . . . . . . . . . . . . . . . . . 3.3 Computational aspects . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Staggered resolution strategy . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Linearization of the mechanical problem . . . . . . . . . . . . . . . . . . 3.3.3 Finite element discretization . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Modeling of hydric shrinkage . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Phase-eld formulation including hydric strains . . . . . . . . . . . . . . 3.5 Numerical applications . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Convergence with respect to mesh size . . . . . . . . . . . . . . . . . . . 3.5.2 Sensitivity with respect to load increments . . . . . . . . . . . . . . . . 3.5.3 Comparison between the linear and nonlinear formulations . . . . . . . 3.5.4 3-point bending test . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.5 Hydric shrinkage of a homogeneous sample . . . . . . . . . . . . . . . . 3.5.6 Hydric shrinkage induced crack propagation . . . . . . . . . . . . . . . . 3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  -eld formulation for nite strains Consider the body Ω, subjected to mixed boundary conditions P • N = t on ∂Ω N (3.6) and u = u on ∂Ω D ,

  is the directional derivative of the residual at u k in the direction of ∆u, with u k+1 = u k + ∆u. This derivative readsD ∆u R(u, δu) = D ∆u Ω S(E(u), d) : δE(u, δu) dΩ ,(3.20) where the superscript indicating the iteration is omitted for notational convenience. Proceeding with calculus leads to D ∆u R(u, δu) = Ω D ∆u S(E(u), d) : δE(u, δu) dΩ + Ω S(E(u, d)) : D ∆u δE(u, δu) dΩ = Ω ∆E(u, ∆u) : ∂ ∂E S(E(u), d) : δE(u, δu) dΩ + Ω S(E(u), d) : ∆δE(∆u, δu) dΩ = Ω ∆E(u, ∆u) : C t M (E(u), d) : δE(u, δu) dΩ + Ω S(E(u), d) : ∆δE(∆u, δu) dΩ,where ∆δE = δF T ∆F + ∆F δF T , with ∆F = ∇(∆u), and C t M is the fourth-order tensor for isotropic material and dened by

2 u

 2 Ω e as ∀Ω e , [k u (u; e)]{∆u i e } = {f u (u; e)} , (3.40) where damage variable d is omitted to simplify notation, and the stiness matrix [k u (u, e)] is decomposed as [k u (u; e)] = [k 1 u (u; e)] + [k 2 u (u; e)] , are interpreted here as the elementary elastic stiness and the geometric stiness, dened as [k 1 u (u; e)] = Ωe g(d) [B t (x, u; e)] T λ(e)[I ⊗ I] + 2µ(e)[I] [B t (x, u; e)] dΩ x; e)] T [B p (S(u), x; e)]dΩ .

  rst and second members in the left-hand side in Eq. (3.25). In the equations above, brackets indicate matrix forms of tensor-valued quantities, and [B p ] is dened in Appendix. A.2.3. Concerning the force vector in Eq. (3.40), we write {f u (x, u; e)} = {f ext u (x, u; e)} -{f int u (x, u; e)} , (3.44) where the rst term, associated with the external force contribution, is expressed as {f ext u (x, u; e)} = Ωe [N u (x, u; e)] T {t(e)}dΩ , (3.45) with {t(e)} the force vector dened by the Neumann boundary condition and applied on element Ω e . The assumption of a well known force ( without tracking force) has been made. The second term related to the internal force is dened as: {f int u (u; e)} = Ωe [B t (x, u; e)] T {S(u; e)}dΩ .

∇

  X d(x) = [B d (x)]{d i e } (3.53) and ∇ X δd(x) = [B d (x)]{δd i e } .

  N d (x)] and[B d (x)] are the matrix of shape functions and their derivatives, respectively. The damage stiness matrix [k d (e)] in element Ω e is dened as

  vector {f d (x)} reads as {f d (e)} = Ωe 2[N d (x)] T HdΩ .

  symbolically denotes here the assembly operator for damage scalar problem relating the elementary operator ([k d (e)], {f d (e)} ) to the global system matrix ([K d ], {F d }).

F

  Figure 3.2 Decomposition of the total deformation gradient.

  u), E h , d) : δE(u, δu) dΩ = 2H (d, δd) + g c (∇d, ∇δd) dΩ = Ω (2H, δd)dΩ .

  (3.71) and (3.64), problem (3.69) can be written as:Ω g(d) h 2 (C : (E(u) -E h )) : δE(u, δu) dΩ = ∂Ω N t • δu dΓ .

  linearization procedure described in Eq (3.3.2) the problem can be formulated as ∀Ω e , [k u (u; e)]{∆u i e } = {f u (u; e)} .

Figure 3 . 4

 34 Figure 3.3 Denition of boundary conditions for M I (left) and M II (right) tests.

Figure 3 . 5 Figure 3 . 6

 3536 Figure 3.5 Material response for M I (left) and M II (right) tests for dierent mesh size h.

Figure 3 . 7

 37 Figure 3.7 Inuence of loading step with respect to the mechanical response.

  Figure 3.8 Mechanical response for linear (Lin) and non-linear (NL) models for dierent values of fracture energy r gc for (left) M I and (right) M II tests.

Figure 3 . 9 Figure 3 .

 393 Figure 3.9 Error on the maximum force response for several fracture energies r gc for (left) M I and (right) M II tests.

Figure 3 .Figure 3

 33 Figure 3.11 Damage evolution M I test (from left to right) for u D = {1.10e -1 , 1.14e -1 , 1.18e -1 } mm, with (top) linear formulation and (bottom) non-linear formulation.

Figure 3 Figure 3

 33 Figure 3.14 Comparison between linear and non-linear 3-PB response with and without damage. The points on the curves relates to the crack path corresponding to the same displacement ūD = {3.00e -2 , 3.27e -2 , 3.33e -2 , 3.50e -2 } mm (see Fig. 3.15).

Figure 3 .

 3 Figure 3.16 Radius evolution with respect to hydric parameter α h for the linear and nonlinear formulations for 2D and 3D.

Figure 3

 3 Figure 3.20 Dierent shrinkage states α h = 0.01% -1% -10%.

  agation of cracks in the clay. The other objective is to identify the mechanisms involved during the shrinkage. An image correlation analysis on the deformation kinematics is also presented and adds valuable additional information to the observation. The presented experiments have been conducted at laboratoire Navier in the context of the post-doctoral project of Abdellali Dadda, also funded by Labex MMCD. The experimental program and the associated numerical simulations presented in chapter 5 were regularly discussed within the context of this collaborative project.

Figure 4 . 1

 41 Figure 4.1 Experimental sample : S 0 , S 1 , S 1c , S 3 and S 6 .

Fig. 4 . 2 ;Figure 4 . 2

 4242 Fig. 4.2;

Figure 4 . 3

 43 Figure 4.3 Experimental setup photos.

  These maps are dened on the set of samples (S 0 , S 1 , S 1c , S 3 , S 6 ), and show the spherical and deviatoric deformation of the samples. The spherical and deviatoric deformations are dened by the Eq.4.1 in percentage. The intensity of these strains are indicated by a color, giving an indication of the kinematics of the drying cracking process in a qualitative and quantitative way.

Figure 4 . 5

 45 Figure 4.5 Denition of certain areas of interest.

Figure 4 .

 4 Figure 4.6 S 0 test with shrinkage, spherical and deviatoric strain for t = 4000, 6000, and 8000 min.

Figure 4 .

 4 Figure 4.7 S 1c test with shrinkage, spherical and deviatoric strain for t = 4000, 6000, and 8000 min.

Figure 4 .

 4 Figure 4.8 S 1 test with shrinkage, spherical and deviatoric strain for t = 4000, 6000, and 8000 min.

Figure 4 .

 4 Figure 4.9 S 3 test with shrinkage, spherical and deviatoric strain for t = 4000, 6000, and 8000 min.

Figure 4 .

 4 Figure 4.10 S 6 test with shrinkage, spherical and deviatoric strain for t = 4000, 6000, and 8000 min.

  4.11, showing three tests with similar positions of the inclusions but very dierent crack geometries. These stochastic behaviours are of course related to the detailed position of the inclusions, but also to some other features of the experiments which are much more dicult to control. In particular, the friction betwwen the lower surface of the sample and the teon substrate involves local uctuations, whichare hard to control. The same holds for the adhesion or debonding of the lateral edge of the sample to the teon wall. The three examples of Fig.4.11 illustrate the potentiel inuence of this interfacial behaviour on the cracking perocess: for the test at the left, the debonding was total, as in S 0 ; in the central text, it was partial as for most other earlier examples; it the right example, there was almost no debonding as in S 6 .

4 .

 4 Elastic and fracture properties are spatially-varying, time-dependent model parameters, due to the heterogeneous nature of the mixture, as well as to the possible complexity of the drying process, when such heterogeneity is taken into consideration.In particular, the local uctuations observed in the central part of sample S 1 might be linked to such constitutive heterogeneities.

Figure 4 .

 4 Figure 4.11 Dierent crack patterns are observed on similar congurations.4.8 ConclusionIn this chapter, we have presented an experimental campaign focusing on the characterization of crack propagation in clay samples lled with rigid inclusions. The analysis relies on the use of a digital image correlation technique, which is used to monitor displacement and deformation elds and crack patterns during the drying process.The observations made on the dierent samples showed several interesting mechanisms. The rst concerns the initiation and propagation of cracks in the presence of heterogeneities. Indeed, a clear pattern was observed for all samples in which damage rst initiates near the inclusions and crack propagation (driven by shrinkage) then occurs. The deformation maps thus obtained also give us additional information about local deformations. In particular, a concentration of hydrostatic and deviatoric strains close to the inclusions is observed, leading to the initiation of damage. Potential sources of uncertainties are nally listed with the aim of explaining the variability observed in the experiments. In particular the adhesion or debonding the the clay at the outer boundary of the sample, and possibly also the friction properties at the interface of sample and substrate might have a prominent role of the deformation and cracking patterns.

Figure 5 . 1

 51 Figure 5.1 Problem denition for a specic conguration: (a) geometry and boundary conditions for the 2D model; (b) 2D mesh; (c) geometry and boundary conditions for the 3D model; (d) 3D mesh.

  .3) with ∆T the transition time, and φ represents the shift of the hyperbolic tangent function. The other parameters are obtained for φ = 10750 and B = 7.0.10 -4 corresponding to ∆T = 1.4.10 3 . The Fig. 5.2 shows both experimental and tted radius function (in black line).

Figure 5 . 2 (

 52 Figure 5.2 (Left) Identied R T function, (right) the identied shrinkage parameter α h .

From a computational

  standpoint, both 2D plane strain and 3D simulations have been carried out. The geometry was discretized with T3 elements in the 2D case, resulting in a mesh with N d dof = 15 × 10 3 degrees of freedom for the damage problem, and N e dof = 2N d dof degrees of freedom for the elastic problem; T4 elements for the 3D case. The mesh contains N e = 400 × 10 3 elements, corresponding to N d dof = 80 × 10 3 degrees of freedom for the damage problem, and to N e dof = 3 × N d dof for the elastic problem. An in-house implementation combining C++ and Matlab algorithms was used. Computational times per time step on a computer equipped with an Intel(R) Xeon (R) platinum 8168 CPU

Figure 5 . 3

 53 Figure 5.3 Crack pattern in S 1c : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

5 . 4 .

 54 In this case, 4 cracks are observed experimentally, while both 2D and 3D simulation results only present 3 cracks. It is noticed that despite the loss of symmetry, the angle between the directions of propagation for two adjacent cracks remains almost constant over crack pairs and is roughly equal to 2π/3.

Figure 5 . 4

 54 Figure 5.4 Crack pattern in S 1 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Figure 5 . 5

 55 Figure 5.5 Crack pattern in S 3 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

  be drawn from the last example (specimen S 6 ), for which 6 inclusions are considered; see Fig.5.6.

Figure 5 . 6

 56 Figure 5.6 Crack pattern S 6 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Figure 5 . 7

 57 Figure 5.7 Spherical strain map in S 1C : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

  Fig. 5.10.

Figure 5 .

 5 Figure 5.10 Spherical strain map in S 6 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Figure 5 .

 5 Figure 5.11 Deviatoric strain map in S 1C : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

S 3

 3 and S 6 congurations Interaction eects are more pronounced for the congurations with three and six inclusions, illustrated in Fig.5.13 and Fig.5.14 respectively. Here, cracks initiate near the inclusions and merge in between, in both the experiments and the simulations.

Figure 5 .

 5 Figure 5.13 Deviatoric strain map in S 3 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Figure 5 .

 5 Figure 5.14 Deviatoric strain map in S 6 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Figure 5 .

 5 Figure 5.15 Deviatoric (left) and spherical (right) strain maps in S 6 : experimental (top) and simulation results (bottom).

  requires improvement to provide satisfactory agreement with the experiments. The possible improvements of the model include: (a) a better description of the experimental boundary conditions, like the friction at the bottom of the specimen and the role of decohesion on the external boundaries; (b) a better description of interfacial conditions around the obstacles within the sample; (c) elucidating the localization of compressive strain near the matrix/interface inclusions; (d) a nonlinear model for model parameters, including a dependence on the water content, and a full phase change description of the clay through its liquid-solid state during the drying.

  ξ) ,ξ 1 ... N i (ξ) ,ξ 1 ... N n (ξ) ξ) ,ξp ... N i (ξ) ,ξp ... N n (ξ) the gradient of f (x) function with respect to ξ ∇ ξ (f (x)) := f (x), ξ 1 f (x), ), x 1 f (x), x 2 .

  ξ)det([J(ξ)])dΩ Ē N g p=1 ω p f (ξ p )det(J(ξ p )),(A.9) p = {1, 2, ..., N g} is the Gauss quadrature points, N g the maximal number of Gauss quadrature, and ω p the gauss weight associated to the N g quadrature chosen.

  e ] and [B u ] operators Those operators can be explicitly describe for 2D (T3 element) as : }G(1, :) + {∇u(4)}G(4, :) + {∇u(7)}G(7, :) {∇u(2)}G(2, :) + {∇u(5)}G(5, :) + {∇u(8)}G(8, :) {∇u(3)}G(3, :) + {∇u(6)}G(6, :) + {∇u(9)}G(9, :) (1/ √ 2)({∇u(3)}G(2, :) + {∇u(6)}G(5, :) + {∇u(9)}G(8, :) + {∇u(2)}G(3, :) + {∇u(5)}G(6, :) + {∇u(8)}G(9, :)) (1/ √ 2)({∇u(3)}G(1, :) + {∇u(6)}G(4, :) + {∇u(9)}G(7, :) + {∇u(1)}G(3, :) + {∇u(4)}G(6, :) + {∇u(7)}G(9, :)) (1/ √ 2)({∇u(2)}G(1, :) + {∇u(5)}G(4, :) + {∇u(8)}G(7, :) + {∇u(1)}G(2, :) + {∇u(4)}G(5, :) + {∇u(7)}G(8, :)) B p ] operator [B p ] operator is dened for 2D case (T3 element) as: S(5)G(3, :) + S(6)G(2, :)] S(2)G(2, :) + (1/ √ 2)[S(4)G(2, :) + S(6)G(1, :)] S(3)G(3, :) + (1/ √ 2)[S(4)G(2, :) + S(5)G(1, :)] p ] operator with shrinkage [B p ] operator is dened for 2D case (T3 element) with : S = C(E -E h ) where E h = α[1, 1, 0] as : (T4 element) with the same way and where E h = α[1, G(7, :) + (1/ √ 2)[S(5)G(9, :) + S(6)G(8, :)] S(2)G(8, :) + (1/ √ 2)[S(4)G(9, :) + S(6)G(7, :)] S(3)G(9, :) + (1/ √ 2)[S(4)G(8, :) + S(5)G(7, :)]

Table 1

 1

	.1 Degradation function examples.
	Note that the degradation function may include a small regularizing parameter 0 < η	1,
	introduced to ensure the well-posedness of the boundary value problem. In this case, a modied
	function g * such that g

* (d) = g(d) + η is considered ( [

  2.3. Mesoscopic Modeling of Crack Propagationgraph of the cost function obtained for a resolution parameter R/R = 3 is shown in the left panel in Fig.2.12.

		0.15									0.9
			KUBC SUBC							0.8
		0.1									0.7
											0.6
	J (γ(g c,m ))	0.05								F max [kN]	0.4 0.5
											0.3
		0									0.2						micro
											0.1						KUBC SUBC
		-0.05									0
		1	1.05	1.1	1.15	1.2	1.25	1.3	1.35	1.4	0	0.002	0.004	0.006	0.008	0.01	0.012
						γ								u [mm]

  9775, with k KU BC = 258.7 [GPa] , δ k KU BC = 10.25% , KU BC = 137.7 [GPa] , δ µ KU BC = 13.30% . SU BC = 238.6 [GPa] , δ k SU BC = 9.3% ,

															(2.33)
	and
				µ (2.34)
	For SUBC, ρ = 0.995, with
				k (2.35)
	1.8	10 -11					3.5	10 -11
	1.6						3
	1.4
							2.5
	1.2
	1						2
	0.8						1.5
	0.6
							1
	0.4
	0.2						0.5
	0						0
	1.5	2	2.5	3	3.5	4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	2.2
						10 11									10 11

and µ SU BC = 112.2

[GPa]

, δ µ SU BC = 10% .

(2.36)

The kernel density estimations of the rst-order marginal and joint distribution for the bulk and shear moduli obtained with the data and the model-based samples are shown in Figs. 2.16 and 2.17 (recall that the resolution is xed by R/R =3 here). A very good agreement is observed

  2.44) for the case of SUBC. The graphs of the normalized correlation functions estimated with the data and with the stochastic model thus identied (under SUBC) are shown in Figs.2.18 and 2.19 for the bulk and shear moduli random elds, respectively.

	1										1
	0.8					R data k R model (τ1) k (τ1)				0.8						R data k R model (τ2) k (τ2)
	0.6										0.6
	0.4										0.4
	0.2										0.2
	0										0
	-0.2										-0.2
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
					Lag distance									Lag distance

Table 3 .

 3

5 Shrinkage induced crack problem parameters.

  Figure 4.4 humidity, temperature and water content experimental conditions.

	humidity [%]
		S 0	S 1	S 1c	S 3		S 6
	1000	2000	3000	4000	5000	6000	7000	8000
				t [min]
	C]
	temperature [ o
		S 0	S 1	S 1c	S 3		S 6
	1000	2000	3000	4000	5000	6000	7000	8000
				t [min]
	100
	water content [%]
		S 0	S 1	S 1c	S 3		S 6
	0
	1000	2000	3000	4000	5000	6000	7000	8000
				t [min]

Table 5

 5

	.1.

Table 5 .

 5 1 Identied parameters for the shrinkage model.

Acknowledgements

Chapter 3. Phase-eld formulation for nite strains and shrinkage The modied operator is described in Appendix A. 2.4. [k 2 u (u; e)] = Ωe g(d) h 2 [G(x; e)] T [B p (S(E(u) -E h ), x; e)]dΩ .

(3.77)

The force vector in Eq. (3.74) becomes:

{f u (x, u; e)} = {f ext u (x, u; e)} -{f int u (x, u; e)} , (3.78) where the rst term, associated with the external force contribution, is still expressed as {f ext u (x, u; e)} = Ωe [N u (x, u; e)] T {t(e)}dΩ .

(3.79)

where {t(e)} is the force vector dened by the assumed xed Neumann boundary condition and applied on element Ω e . The second term related to the internal force is dened as:

(3.80)

or explicitly written as:

The displacement increment ∆u is thus obtained with shrinkage contribution by solving the linear problem

where [K u (u)] and {F u (u)} are obtained by assembling elementary matrices:

(3.83)

where A symbolically denotes the assembly operator for displacement problem. The phase eld discretization does not change and is formulated in section 3.3.3.

Numerical applications

Some computational tests will be presented, including a sensitivity analysis of the mesh, load increment, followed by a comparison test between linear and nonlinear modeling of the cracking model on sandard tests such as traction and pure shear. Finally, the last tests using the drying model will be presented.

Convergence with respect to mesh size

We rstly investigate a S = 1 × 1 mm 2 square plate with a half notch placed on the left side of the specimen. The setup of the problem is shown in Fig. 3.3 for two boundary conditions. The rst problem is denoted by "M I " and is used to study the full traction problem attached to the mode-I fracture. The second example, referred to as "M II " is related to a mode-II fracture (see the crack paths on Fig. 3.4).

Parameters of the system and boundary conditions are summarized in Table 3.1.

S 1 conguration

Similar results are shown in Fig. 5.8 for the conguration with the o-centered inclusion.

Figure 5.8 Spherical strain map in S 1 : experimental results (top), numerical results in 2D (middle), and numerical results in 3D (bottom) at times t 1 , t 2 , and t 3 (from left to right).

Here, we observe the same characteristics regarding the amplitudes as in the previous case, both globally and at the edges, with some local dierences due to the lack of the geometry of revolution though. Indeed, we notice that the zone at the edges of the inclusion and the nearest outer edge zone are associated with localized extensive strains when damage initiates (t = t 1 ).

This zone does not seem to aect or to predetermine a path, moreover with this geometry, the 3D case is quantitatively closer to 2D, regarding t=t1.

S 3 conguration

For the 3-inclusion geometry presented in Fig. 5.9, interaction eects are observed in both 2D and 3D simulations at t = t 1 , with specic high-intensity areas located at the edges of the inclusions and facing each other. The crack paths follow these zones at t 2 and t 3 . Branching seems to occur at the center of the sample for t = t 2 in the 2D case, with the creation of a positive deformation zone (rouge). these areas of strong extension occurs exactly in the damaged areas (d ≥ 0.8) which stretches at t = (t 2 , t 3 ). It should be highligted that the positive deformations are here concentrated around the inclusions, with a loss of symmetry of revolution, whereas they are distributed in the experiment and preserve this symmetry locally around the inclusions. A complex phenomenon in the experiment causes a mechanical deformations distributed and which is opposed to the simulation. An explanation could come from a non-linear elastic behaviour law that would disadvantage high deformations (as in elastomer materials) or from an unknown complex friction eect.

Conclusion and Perspectives

General conclusion

In this thesis, we have investigated the fracture of heterogeneous quasi-brittle materials through both computational analysis and physical experiments, putting emphasis on random media and fracture induced by desiccation. The contributions of the present work are summarized as follows.

First, a stochastic approach to model crack propagation in random media was proposed in Chapter 2. The formulation relies on a phase-eld model involving mesoscopic material coecients. Monte-Carlo simulations were rst performed at microscale to characterize stochasticity in the macroscopic response. The denition of the mesoscopic parameters was then addressed.

The elasticity eld at mesoscale was dened as the isotropic approximation of spatially varying apparent tensors, obtained under kinematically and statically uniform boundary conditions. The mesoscopic toughness was identied solving an inverse problem related to the mean peak force.

A probabilistic model for the elasticity random eld was then constructed and allows for a fast, robust sampling of mesoscopic elasticity. The results obtained by combining this surrogate with the phase-eld model were nally compared with the reference, microscopic model. It was shown, in particular, that the model-based, mesoscopic phase-eld formulation associated with statically uniform boundary conditions allows for an accurate prediction of both the mean elastic response and mean peak force.

In Chapter 3, we reviewed a nite-strain phase-eld formulation, with the aim of capturing the eects of the large deformations occurring during the drying of clay samples. A parametric analysis was then conducted in order to quantify the impact of the numerical parameters in terms of convergence and accuracy of the solution. We then extended the model to account for the hydric eects arising during the drying tests. A generic, low-dimensional representation was retained to parameterize the driving term, hence allowing the choice of a calibration from the available experimental data.

Chapter 4 is concerned with the description of the physical experiments (conducted by Dr.

A. Dadda, Navier Laboratory, Ecole des Ponts ParisTech). In these experiments, cylindrical clay samples are dried and cracks start propagating due to shrinkage. In order to introduce and characterize the impact of potential interaction eects in a controlled manner, xed cylindrical obstacles were inserted within the samples. In addition to crack patterns, strain elds were also obtained as a function of time, using a digital image correlation technique. The combination of these results constitutes an original contribution enabling us to discuss some correlation eects between the nature (and evolution) of the strain elds and the creation of crack paths.

The experimental results were nally compared with numerical predictions in Chapter 5. The bijective geometrical transformation from the reference element to physical element have some properties as : local coordinate form:

x = [x 1 (ξ), x 2 (ξ)] T = [x 1 (X), x 2 (X)] T , in 2D case, (A.1)

x = [x 1 (ξ), x 2 (ξ), x 3 (ξ)] T = [x 1 (X), x 2 (X), x 3 (X)] T , in 3D case.

(A.2) which the following interpolation:

x p (ξ) = n i=1 N i (ξ)X (i) p , (A.3)