Ecole Sigle

Doctorale

Renée El

Melhem Bat

Blaise Pascal

Insa : R Gourdon

M Stéphane

M Gérard Scorletti
email: gerard.scorletti@ec-lyon.fr

Isabelle Mme

Von Bueltzingloewen

Keywords: composites, MultiGrid, laminated materials, brittle crack, phase field matériaux composites, Multigrille, stratifié, rupture quasi-fragile, méthode de champ de phase

With the development of imaging techniques like X-Ray tomography in recent years, it is now possible to take into account the microscopic details in composite material simulations. However, the composites' complex nature such as inclined and broken fibers, voids, requires rich data to describe these details and thus brings challenging problems in terms of computational time and memory when using traditional simulation methods like the Finite Element Method. These problems become even more severe in simulating failure processes like crack propagation. Hence, it is necessary to investigate more efficient numerical methods for this kind of large scale problems.

The MultiGrid (MG) method is such an efficient method, as its computational cost is proportional to the number of unknowns. In this thesis, an efficient MG solver is developed for these problems. The MG method is applied to solve the static elasticity problem based on the Lamé's equation and the crack propagation problem based on a phase field method. The accuracy of the MG solutions is validated with Eshelby's classic analytic solution. Then the MG solver is developed to investigate the composite homogenization process and its solutions are compared with existing solutions in the literature. After that, the MG solver is applied to simulate the free-edge effect in laminated composites. A real laminated structure using X-Ray tomography is first simulated. At last, the MG solver is further developed, combined with a phase field method, to simulate the brittle crack propagation. The MG method demonstrates its efficiency both in time and memory dimensions for solving the above problems.

List of Figures

) i, j,k (white boxes stand for the location of Q), black circles stand for the location of u associated and white circles stand for the location of left u). . . . . . 2.5 Finite difference discretization for ∂ ∂y (Q ∂u ∂x ) i, j,k (white boxes stand for the location of Q), black circles stand for the location of u associated and white circles stand for the location of left u). . . . . . a.error of initial guess (scale: 10 1 ), b.error after 10 relaxations on finest grid (scale: 10 0 ), c.error after the interpolation (scale: 10 -1 ), d.error after the 10 post relaxations on finest grid (scale: 10 -2 ). . .

Restriction operator I H

h and its weighting [START_REF] Boffy | Techniques multigrilles et raffinement pour un modèle 3D efficace de milieux hétérogènes sous sollicitations de contact[END_REF]

Background and motivation

Composite materials consist of two or more constituent materials with significantly different properties and exhibit far more improved qualities that neither of the constituent materials possesses [START_REF] Robert | Mechanics of composite materials[END_REF]. The advantages of using composites include high strength, light weight, improved corrosion resistance, etc. Nowadays composites like fibrous composite and laminate composite have been widely used in aircraft structures for the past decades. For example, more than 50% of the overall airframe in the Boeing 787 Dreamliner is made of composite materials, such as fuselage and wings, as shown in Fig 1 .1. The usage of the composite materials significantly reduces the aircraft's weight and hence the fuel consumption. Furthermore, the aircraft's strength and its resistance to damage and corrosion are improved due to the optimized properties, which further reduces the operating cost and improves the efficiency in the long term.

However, with the increased percentage of usage, the safety concerns with respect to composite materials in the airframe are rising as well, which stem mainly from the lack of information on the behavior of composite structures since they are a relatively new material compared to conventional metallic materials [START_REF] Baker | Composite materials for aircraft structures[END_REF]. The behavior of a large composite component undergoing complex external loading is difficult to predict and designers have to rely on their limited experience and take large safety factors. For example, the wing box in the Boeing 787, which is made of composite materials, had initially insufficient stiffness. Boeing had to add additional brackets to reinforce the structure of the wing box for those already built. For ones that are yet to be built, Boeing has to modify the design and delay the delivery.

Hence, numerous efforts are paid by researchers to investigate the performance prediction of the composites both in analytical and numerical ways. The difficulties in composite performance prediction stem mainly from the highly anisotropic microstructures. This anisotropy makes the classic mechanics with respect to homogeneous materials like metallic alloys improper any more for the composites. From an analytical point of view, it is possible to develop the prediction models at a microscopic scale directly [START_REF] Hashin | Analysis of composite materials: a survey[END_REF][START_REF] Mura | Micromechanics of defects in solids[END_REF][START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]. However, it is often far too complex to consider all the microstructures in a large composite component at the same time. An alternative way is to regard a composite structure as a "homogeneous" one at the macroscopic scale and obtain its global effective properties. This kind of approximation is commonly referred to as the homogenization methods [START_REF] Hashin | Analysis of properties of fiber composites with anisotropic constituents[END_REF].

The Homogenization method The homogenization method refers to the ap-proximation process that the global behavior of a composite structure can be obtained through its microstructural stress-stain field. Eshelby proposed a terminology "eigenstrain" that describes a kind of strain produced without external forces in an elastic, linear, homogeneous and infinite solid [START_REF] Jd Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Jd Eshelby | The elastic field outside an elastic inclusion[END_REF][START_REF] Jd Eshelby | Elastic inclusions and inhomogeneities[END_REF]. According to this idea, a so-called "Eshelby's equivalent inclusion method" was developed to determine the analytical solution of the stress field in a heterogeneous material containing of simple shaped inhomogeneities (ellipsoid, cuboid, ...) [START_REF] Mura | Micromechanics of defects in solids[END_REF]. Based on this, various rigorous homogenization methods for predicting the overall material properties have been developed. Hashin and Shtrikman [HS62a,[START_REF] Zshtr Hashin | A variational approach to the theory of the elastic behaviour of polycrystals[END_REF][START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] proposed an upper and lower bound for the effective moduli of composites by assuming that the inhomogeneity is spherical and is bonded perfectly to the matrix. The Mori-Tanaka method, which was proposed by Mori and Tanaka in 1973 [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF] and reformulated by Benveniste in 1987 [START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF], estimates the effective stiffness tensor by considering cylindrical, ellipsoidal or plane fibers or fractures embedded in an isotropic matrix transversely isotropic or orthotropic. Another method called self-consistent method [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF][START_REF] Budiansky | On the elastic moduli of some heterogeneous materials[END_REF][START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] takes account of the interaction between inhomogeneities by replacing the matrix domain with an effective medium. These methods are reviewed in detail in Chapter 3. More information about the homogenization methods in composites can be found in [START_REF] Oller | Numerical simulation of mechanical behavior of composite materials[END_REF].

From a numerical point of view, with the development of the computational ability and simulation methods in the last two decades, the numerical simulation is playing a more and more important role in the composite structure design since the analytical methods are limited to the simple geometry and the number of components in composites. Usually, an engineering product before commercial use goes through three main steps: design, prototyping and physical testing. In the aircraft industry, the steps of prototyping and physical tests are expensive and sometimes the testing time may not be compatible with human dimensions. With the help of the numerical simulations, designers are able to predict the performance, analyze reliability and potential failures, optimize construction, and export accurate information to manufacturing, all before a physical prototype with expensive composite materials is built [ANS].

Hence, a variety of numerical methods have been developed to simulate the composite materials. Among them, the most commonly mentioned methods are the finite element method and the finite difference method. Many mechanical problems in composite performance prediction can be described in terms of partial differen- tial equations (PDEs). The core of these numerical methods is to solve the PDEs approximately on a discretized domain.

The finite element method (FEM) The FEM was first developed in the 1950s in the aerospace industry. It is a numerical approximation method that subdivides the problem domain into smaller parts, called finite elements, where the simpler equations are developed and then assembled into a larger system of equations that model the entire problem [START_REF] Fish | A first course in finite elements[END_REF]. Fig 1 .2 (a) shows an example of a plate with a circular inclusion, discretized into smaller triangular elements. FEM is based on the energy principle such as the virtual work principle or the total potential energy principle. The system description of this method in solid mechanics is referred in O.C Zienkiewicz and R.L Taylor's book "The Finite Element Method: Solid Mechanics" [START_REF] Olek | The finite element method for solid and structural mechanics[END_REF]. The advantages of using FEM in composite structure analysis include:

• The FEM can readily handle the complicated geometries because of the abundant element types. Theoretically, any geometry can be approximated by triangular elements.

• Boundary conditions are easy to impose in the FEM.

• The FEM are used by millions of engineers and scientists, hence there are al- However, the computational time and memory cost are the problems in the FEM for simulating the composite materials with complex structures since the microstructure needs to be modeled accurately. For example, the fiber reinforced laminate, which is also the application in this work, consists of layers with different fiber orientations, where in each layer fibers are aligned along one direction (see Fig 1 .3). The coupling of different orientation layers results in the specific material properties, that differ in terms of the location. The number of layers in a single laminate structure can be dozens, hundreds, or even thousands. The memory cost in the FEM is proportional to the N 2 , where N is the number of degree of freedom (DoF). Furthermore in the step of solving system equations, the cost of computational time can be unacceptable for a large N, e.g N = 10 9 , since the complexity for obtaining the inverse matrix increases dramatically with the increase of N. Therefore, it can be a big challenge to use the FEM in composite simulations.

The finite difference method (FDM) The FDM is another method for solving PDEs, which began to be applied in numerical simulations since the early 1950s. The principle of the FDM is to use differential quotients upon the Taylor expansion to approximate the derivatives in PDEs. The approximate solution is solved using iterative methods like Jacobi method or Gauss-Seidel method on the discretized domain. as in the FEM example. The computation of solutions is carried out on the vertexes of rectangular meshes. A full description of the FDM can be found in [START_REF] Gordon | Numerical solution of partial differential equations: finite difference methods[END_REF].

The FDM can be more efficient than the FEM if the problem can be solved with a uniform structured grid. Compared to the FEM, the memory required in the FDM is less since it is proportional to the number of DoF N, which provides an opportunity to use dense meshes to depict the structural details in composite simulations. Nevertheless, after a few number of initial iterations, the convergence rate in the FDM slows down and deteriorates even worse when mesh size decreases. The reason is that the errors eliminated by the standard iterative methods (Jacobi, Gauss-Seidel, SOR) mainly belong to the short-wavelength components and then become smooth after the first few iterations. The standard iterations are inefficient to eliminate the long-wavelength (smooth) errors. Therefore, a naturally view is to treat these different scale errors differently, which drives an efficient numerical method-the MultiGrid method [START_REF] Yavneh | Why multigrid methods are so efficient[END_REF].

The MultiGrid method (MG) The MG method is an iterative acceleration method based on the FDM by solving the problem on different mesh sizes. It results in a hierarchy of meshes and the short-and long-wavelength errors are eliminated separately on the fine grid and coarse grid (see Fig 1 .4, same example as in the FEM and FDM illustration parts). Therefore, the MG method not only has the memory efficiency as in the FDM but also exhibits the computing time efficiency. According to the different coarsening processes, the MG method can be distinguished between algebraic MG (AMG) [START_REF] Jw Ruge | Algebraic multigrid[END_REF] and geometric MG. If there is no special indication, the term MG in the following sections refers to the geometric MG method.

The MG method is a relatively new method compared to the FEM and FDM. The first papers about the MG [START_REF] Petrovich | A relaxation method for solving elliptic difference equations[END_REF][START_REF] Petrovich | The speed of convergence of one iterative process[END_REF] were published in 1961 and 1964 by Fedorenko, who formulated a MG algorithm for solving the Poisson equation. After then, it was in the mid-seventies that Achi Brandt [START_REF] Brandt | Multi-level adaptive technique (mlat) for fast numerical solution to boundary value problems[END_REF][START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF][START_REF] Brandt | Multi-level adaptive techniques (mlat) for partial differential equations[END_REF] clearly pointed out the efficiency of the MG and proposed the main principles of the MG. The pioneering work by Brandt made the MG into practice and attracted wide attention to the MG research. His outstanding contributions include the introduction of a nonlinear multigrid method (FAS) and of adaptive techniques (MLAT), the discussion of general domains, the systematic application of the nested iteration idea (FMG) and, last but not least, the provision of the local Fourier analysis tool for theoretical investigation and method design [START_REF] Trottenberg | Multigrid. Academic press[END_REF]. An introduction to the MG method can be found in [BM + 00, Wes95]. More thorough books for the MG method are referred to see Brandt's guide book "Multigrid techniques: 1984 guide with applications to fluid dynamics" [BL11] and Trottenberg's comprehensive textbook "MULTIGRID" [START_REF] Trottenberg | Multigrid. Academic press[END_REF]. For more specific applications using the MG method, e.g the ElastoHydrodynamic Lubrication (EHL) problem, readers should refer to [START_REF] Venner | Multi-level methods in lubrication[END_REF]. A great number of MG source code and information in different fields can be found on the website MGNET: http://www.mgnet.org.

The application of the MG method in solid mechanics started from solving the membrane problems by G. Brand and D. Braess in the later 1980s [START_REF] Brand | Multigrid methods in finite element analysis of plane elastic structures[END_REF][START_REF] Braess | On the combination of the multigrid method and conjugate gradients[END_REF][START_REF] Braess | A multigrid method for the membrane problem[END_REF]. Parsons et al. [START_REF] Id Parsons | The multigrid method in solid mechanics: part iâ ȂŤ algorithm description and behaviour[END_REF][START_REF] Id Parsons | The multigrid method in solid mechanics: part iiâ ȂŤ practical applications[END_REF] investigated the performance of the MG method in 3D dimensional solid mechanics for homogeneous structures and concluded that the convergence rate is independent of the problem size but the MG performance is sensitive to the features of the problem, such as the degree of mesh non-uniformity. Fish et al. [START_REF] Fish | Multigrid method for periodic heterogeneous media part 1: Convergence studies for one-dimensional case[END_REF][START_REF] Fish | Multi-grid method for periodic heterogeneous media part 2: Multiscale modeling and quality control in multidimensional case[END_REF] first applied the MG method in numerical composite homogenization with periodic structures based on the finite element meshes and revealed that the convergence rate increases with the increase in material heterogeneity. M. F Adams [Ada99, Ada02, AT00, Ada04] focused on the algebraic unstructured grids which can be easily used in FE meshes and developed parallel algorithms to improve the efficiency. Recently, Watremetz et al. [START_REF] Watremetz | 2d thermo-mechanical contact simulations in a functionally graded material: A multigridbased approach[END_REF] applied the MG solutions to 2D graded materials and Boffy et al. developed an efficient MG solver for contact problems in 3D heterogeneous materials [START_REF] Boffy | An efficient 3d model of heterogeneous materials for elastic contact applications using multigrid methods[END_REF][START_REF] Boffy | Multigrid numerical simulation of contact mechanics of elastic materials with 3d heterogeneous subsurface topology[END_REF], considering multiple moving heat sources [START_REF] Boffy | Detailed modelling of a moving heat source using multigrid methods[END_REF] and for strongly heterogeneous materials [START_REF] Boffy | Multigrid solution of the 3d stress field in strongly heterogeneous materials[END_REF].

Objectives

The motivation of this thesis comes from the huge demand of large scale computing in composite simulations. The complex nature of the composites requires large data points to describe the micro-structures. And the development of imaging techniques like X-Ray tomography in recent years allows to consider the real structure in numerical simulations. However, as point in the previous parts, the developments of simulation methods fall behind this demand.

The efficiency of the MG method both in memory and time dimensions drives us to develop an efficient MG model for the composite material simulations. However, the related literature is limited compared to the one on FEM. To the author's knowledge, the performance of the MG method in composite simulations in detailed fields, such as the accuracy of the MG solution, the numerical homogenization, the free edge effects in laminated structure and the damage model, has not been investigated explicitly. Hence, the purpose of this thesis is to investigate the MG method in composite simulations and show the possibility of solving composite problems. The objective can be divided into the following sub-objectives:

(1) Develop an efficient MG solver for simulating the 3D heterogeneous materials in an elastic solid mechanics field. This work is partly based on H. Boffy's thesis [START_REF] Boffy | Techniques multigrilles et raffinement pour un modèle 3D efficace de milieux hétérogènes sous sollicitations de contact[END_REF].

(2) Once the efficient MG solver developed, the accuracy of the MG solutions should be investigated. The influence of the mesh size, the ratio of the 1.3. Outline material properties and the convergence rate should be discussed in detail.

(3) Develop numerical homogenization models using the MG method. The computational cost in traditional numerical homogenization models using the FEM method is expensive due to the complex nature in composites.

The advantages of the MG method shows the potential in this application. The boundary conditions in the MG model are discussed as well.

(4) Investigate the laminate structure and the related free edge effect using the MG method. A real laminate structure obtained from X-ray tomography is modeled using this MG solver. The influences of fiber orientation, interface layer thickness, fiber layer thickness as well as the different ratios of the Young's modulus between fibers and matrix are studied.

(5) Investigate damage models in composite materials and concentrate on the phase field methods. Develop a MG model combining with the phase field methods to simulate the crack propagation in heterogeneous materials. Apply the failure model to the laminate structure and investigate the delamination process.

Outline

Following the objectives outlined in previous sections, the structure of the thesis is elaborated as follows:

Chapter 1: In this chapter, the background, motivation and the objectives of the thesis are specified. Section 1.1 explains the background and the reason why select the MG method as the numerical tool in this thesis. The advantages of the MG method over the FEM and the FDM in composite simulations: time and memory, are illustrated. Detailed objectives of this thesis are proposed in Section 1.2.

Chapter 2: This chapter first reviews the basics of composite mechanics, the finite difference method and the basic notations of the MG method in solid mechanics. The linear elastic equation: Lamé equation that describes approximately the microstructure is presented and a proper finite difference scheme is derived in Section 2.2. Then in Section 2.3, some basic notations in MG: smoothing solver, transfer operators and coarse grid operator are explained. An example of a twogrid cycle is presented to show how it works for the Lamé equation to smooth the errors. After that, boundary conditions encountered in this thesis are classified in Section 2.4 and the related treatment in the MG routine is explained as well. Then in Section 2.5, a powerful analytical tool for analyzing the performance of the MG method is reviewed and the practical analysis for the Lamé equation is made. This shows an ideal convergence rate for the Lamé equation. In Section 2.6, an advanced MG technique: local grid refinement (LGR) is presented, which can help to decrease the computational time and memory consumption in MG. Finally, Section 2.7 reviews Alcouffe's techniques [START_REF] Re Alcouffe | The multigrid method for the diffusion equation with strongly discontinuous coefficients[END_REF] for restoring the good efficiency when facing the strong discontinuous problems in composite structures.

Chapter 3: In this chapter, the emphasis is put on the investigation of accuracy of the MG solutions and on the performance of the MG homogenization models. Section 3.2 presents an MG model that combines all the techniques mentioned in Chapter 2 to solve a simple but fundamental problem in composite mechanics: an inhomogeneity within a half-space. The MG solutions are compared with Eshelby's analytical solutions. Then the MG model is developed to be used in the homogenization application. Several analytical homogenization theories available in the literature are reviewed in detail in Section 3.3. The MG homogenization results are compred with these analytical methoods. The influence of the mesh size, the ratio of the material properties and the convergence rate are discussed in detail in Section 3.4. Finally, the MG homogenization model is used to investigate a more complex but common structure: woven material in Section 3.5.

Chapter 4: This chapter focuses on the application in laminate composites, which is the main form in composite materials. The famous stress concentration effect: free edge effect is considered in this simulation. Firstly, Section 4.2 reviews the free edge effect in literature and a comparison is made between the finite element homogenized model and the MG microstructural model. Then in Section 4.3, a real laminate structure obtained through X-ray tomography is simulated, which takes over 200 million DoFs. At last, the parameters that influence this phenomenon including fiber orientation, interface layer thickness, fiber layer thickness as well as the different ratios of the Young's modulus between fibers and matrix, are investigated in Section 4.4.

Chapter 5: In this chapter, the failure model in composite materials is investi-gated. The MG solver is further developed, combined with a phase field method, to simulate the crack propagation in brittle materials. First Section 5.2 reviews the basic notations in phase field method and a phase field method used in this work is illustrated in detail. Then in Section 5.3, the governing equations for both stress equilibrium and phase field is summarized. An overall algorithm suitable in the framework of MG is built. After that, some treatment in MG to efficiently solve the phase field equation are presented in Section 5.4. At last, three numerical crack propagation examples are presented in Section 5.5. The efficiency of the MG solver for simulating the brittle crack propagation using the phase field is pointed out as well.

Chapter 6: The last chapter summarizes the main results obtained in this thesis and suggests some directions for further work.

Introduction

In this chapter, the basics of linear elasticity theory and applied numerical simulation method -MultiGrid are reviewed. We present the basic notations in elasticity in Section 2.2 and emphasize the so called Lamé equation. A proper finite difference discretization scheme for the Lamé equation is derived. Then in Section 2.3, we review some fundamental notions in MG: smoothing solver, transfer operators and coarse grid operator. A two-grid cycle is presented and shows how it smooths errors efficiently for the Lamé equation. In Section 2.4, we present several relevant boundary conditions encountered in this thesis and explain how to deal with them in MG. Section 2.5 describes a powerful tool -local mode analysis for the quantitative analysis and the design of efficient MG methods with respect to the Lamé equation. Local grid refinement is presented in Section 2.6 to decrease the computational time and memory consumption. At last in Section 2.7 we focus on one of the main topics of this thesis -strong discontinuities and special treatments in MG are discussed to restore the good efficiency.

Linear elasticity theory and discretization

In linear elastic solid mechanics, an elastic body is one that deforms under external forces but can revert to the initial shape after removing the external forces [START_REF] Timoshenko | Theory of elasticity[END_REF]. Linearity refers to the relation between stress and strain, which is approximated linearly when deformation is small [START_REF] Ronald | The physics of rubber elasticity[END_REF] and usually called Hooke's law. In the following part of this section, the notations of strain, stress and Hooke's law are rapidly reviewed in order to derive a proper difference scheme at the end of this

Strain

Strain in elasticity is described as the ratio of deformation to the original dimension of the body. It can be decomposed into normal strain and shear strain (see Fig 2 .2). The normal strain, which is also called the extensional strain, is expressed as the change in length compared to the original length. The shear strain is defined as the tangent of the angle changed with respect to two specific directions. According to this definition, in the 3D case, it has three normal strain components (ε xx , ε yy and ε zz ) and six shear strain components (γ xy , γ xz , γ yx , γ yz , γ zx , γ zy ). Due to the symmetry, it has: γ xy = γ yx , γ yz = γ zy and γ xz = γ zx . So the strain tensor in 3D is:

ε =             ε xx γ xy γ xz γ xy ε yy γ yz γ xz γ yz ε zz             (2.1)
             (2.2)
or:

ε i j = 1 2 (u i, j + u j,i ). (2.3)

Stress

In mechanics, stress is defined as the force acting on a unit area. As in strain, stress can be divided into normal stress and shear stress (see Fig 2 .3). The normal stress is the force per unit area acting perpendicularly to a selected area, while a shear stress is a force per unit area acting transversely to the area. In 3D case, it has also three normal stress components (σ xx , σ yy and σ zz ) and three symmetric shear stress components (τ xy , τ yz τ xz ,). The stress tensor can be expressed in a matrix way:

σ =             σ xx τ xy τ xz τ xy σ yy τ yz τ xz τ yz σ zz             (2.4)
If we consider a static body without any body force, then the surface traction over the surface S is 0, which gives:

S σndS = 0. (2.5)
where n is the normal vector to the surface S . According to the Gauss divergence theorem, the surface integral can be converted to a volume integral:

S σndS = V ∇σdV = 0. (2.6)
Since the volume V is arbitrary, this requires that the integral be zero:

∇σ = 0. ( 2
= 0              (2.8)
or σ i j, j = 0 i, j = 1, 2, 3.

(2.9)

Hooke's law

The relations between stress and strain expressed in a general way are:

σ(ε) :                                      σ xx = f 1 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ) σ yy = f 2 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ) σ zz = f 3 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ) τ xy = f 4 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ) τ yz = f 5 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ) τ xz = f 6 (ε xx , ε xx , ε yy , ε zz , γ xy , γ yz , γ xz ).
(2.10)

where f i (i = 1, 2, ..., 6) depends on the inherent elastic body properties. With the small deformation assumption, equations 2.10 can be approximated by means of Taylor expansions. If we eliminate the second order and higher terms, for example, the first one in equations 2.10 can be written as follows:

σ xx = f 1 | 0 + ∂ f 1 ∂ε xx | 0 ε xx + ∂ f 1 ∂ε yy | 0 ε yy + ∂ f 1 ∂ε zz | 0 ε xx + ∂ f 1 ∂γ xy | 0 γ xy + ∂ f 1 ∂γ yz | 0 γ yz + ∂ f 1 ∂γ xz | 0 γ xz .
(2.11)

where f 1 | 0 denotes the value of f 1 when ε is 0 and stands for the initial stress in fact. Therefore, if there is no initial stress, Equations 2.10 can be approximated as below:

                              σ xx σ yy σ zz τ xy τ yz τ xz                               =                               C 1111 C 1122 C 1133 C 1112 C 1123 C 1113 C 2211 C 2222 C 2233 C 2212 C 2223 C 2213 C 3311 C 3322 C 3333 C 3312 C 3323 C 3313 C 1211 C 1222 C 1233 C 1212 C 1223 C 1213 C 2311 C 2322 C 2333 C 2312 C 2323 C 2313 C 1311 C 1322 C 1333 C 1312 C 1323 C 1313                                                             ε xx ε yy ε zz γ xy γ yz τ xz                              
(2.12) Equation 2.12 is called the generalized Hooke's law, where the coefficients C i jkl (i, j, k, l = 1, 2, 3) are elastic constants. Hooke's law can be written in a tensor form as: σ = Cε or σ i j = C i jkl ε kl i, j, k, l = 1, 2, 3.

(2.13)

C is called the stiffness tensor. For anisotropic materials, a maximum of only 21 of 36 elastic constants are independent because of symmetry, which means C i jkl = C kli j . While for isotropic materials, only 2 elastic constants are independent, which leads to:

                              σ xx σ yy σ zz τ xy τ yz τ xz                               =                               C 1111 C 1122 C 1122 0 0 0 C 1122 C 1111 C 1122 0 0 0 C 1122 C 1122 C 1111 0 0 0 0 0 0 (C 1111 -C 1122 ) 2 0 0 0 0 0 0 (C 1111 -C 1122 ) 2 0 0 0 0 0 0 (C 1111 -C 1122 ) 2                                                             ε xx ε yy ε zz γ xy γ yz τ xz                               (2.14)
If we define: C 1122 = λ, C 1111 -C 1122 = 2µ, then 2.14 can be simplified to:

σ i j = λε kk δ i j + 2µε i j i, j, k = 1, 2, 3. (2.15)
where δ i j is the Kronecker delta, which is:

δ i j =          1 i = j 0 i j.
(2.16)

λ and µ are called Lamé coefficients, they are a function of the Young's modulus E and the Poisson ratio ν.

λ(x, y, z) = E(x,y,z)ν(x,y,z) (1+ν(x,y,z))(1-2ν(x,y,z)) µ(x, y, z) = E(x,y,z) 2(1+ν(x,y,z)) .
(2.17)

Finite difference discretization

From Equations 2.9 and 2.15, we obtain:

σ i j, j = (λε kk δ i j + 2µε i j ) , j . (2.18)
With Equation 2.3, Equation 2.18 can be written using displacements:

(λu j, j ) ,i + (µu i, j ) , j + (µu j,i ) , j = 0.

(2.19)

In detail:

∂ ∂x ((λ + 2µ) ∂u ∂x ) (a1) + ∂ ∂y (µ ∂u ∂y ) (a2) + ∂ ∂z (µ ∂u ∂z ) (a3) + ∂ ∂x (λ ∂v ∂y ) (a4) + ∂ ∂y (µ ∂v ∂x ) (a5) + ∂ ∂x (λ ∂w ∂z ) (a6) + ∂ ∂z (µ ∂w ∂x ) (a7) = (2.20a) ∂ ∂y ((λ + 2µ) ∂v ∂y ) (b1) + ∂ ∂x (µ ∂v ∂x ) (b2) + ∂ ∂z (µ ∂v ∂z ) (b3) + ∂ ∂y (λ ∂u ∂x ) (b4) + ∂ ∂x (µ ∂u ∂y ) (b5) + ∂ ∂y (λ ∂w ∂z ) (b6) + ∂ ∂z (µ ∂w ∂y ) (b7) = (2.20b) ∂ ∂z ((λ + 2µ) ∂w ∂z ) (c1) + ∂ ∂x (µ ∂w ∂x ) (c2) + ∂ ∂y (µ ∂w ∂y ) (c3) + ∂ ∂z (λ ∂u ∂x ) (c4) + ∂ ∂x (µ ∂u ∂z ) (c5) + ∂ ∂z (λ ∂v ∂y ) (c6) + ∂ ∂y (µ ∂v ∂z ) (c7) = (2.20c) i -1 i i + 1 j -1 j j + 1 h x h x h y h y h x 2 h x 2 (i, j, k) k -0 k -1 k + 1 x y z Figure 2.4: Finite difference discretization for ∂ ∂x (Q ∂u ∂x ) i, j
,k (white boxes stand for the location of Q), black circles stand for the location of u associated and white circles stand for the location of left u).

i -

1 i i + 1 j -1 j j + 1 h x h x h y h y (i, j, k) k -0 k -1 k + 1 x y z Figure 2.5: Finite difference discretization for ∂ ∂y (Q ∂u ∂x ) i, j
,k (white boxes stand for the location of Q), black circles stand for the location of u associated and white circles stand for the location of left u).

There are two types of derivatives: one is ∂ ∂x (Q ∂u ∂x ) like a1, a2, a3, b1, b2, b3, c1, c2, c3 and another is ∂ ∂y (Q ∂u ∂x ) like a4, a5, a6, a7, b4, b5, b6, b7, c4, c5, c6, c7 in equation 2.20. There are several choices for the discritization process, like the Finite Element Method, the Finite Difference Method or the Finite Volume Method. Here we use the Finite Difference Method to discretize them (see Fig 2.4 and 2.5) [START_REF] Richard | Matrix iterative analysis[END_REF].

Assuming a 3D cubic domain Ω: (0, X) × (0, Y) × (0, Z) filled with identical but smaller cubic grids. The dimensions of each smaller cube are h x × h y × h z , where h x = X n x , h y = Y n y and h z = Z n z and n x , n y , n z are the number of grids in three directions respectively. Therefore, the relation between grid numbers and the geometry is:

(x i , y j , z k ) = (0 + i * h x , 0 + j * h y , 0 + k * h z ).
(2.21)

where i = 0, 1, 2, ..., nx, j = 0, 1, 2, ..., ny and k = 0, 1, 2, ..., nz. For the first kind of item ∂ ∂x (Q ∂u ∂x ), a short second order central discretization at point (i, j, k) is:

( ∂ ∂x (Q ∂u ∂x )) i, j,k = (Q ∂u ∂x ) i+1/2, j,k -(Q ∂u ∂x ) i-1/2, j,k h x + O(h 2 x ). (2.22)
with O(h 2 x ) representing an error from truncation of the Taylor expansions (about the truncation error, detail please see pages 19-21 in book [START_REF] Venner | Multi-level methods in lubrication[END_REF]). Furthermore, (Q ∂u ∂x ) i+1/2, j,k and (Q ∂u ∂x ) i-1/2, j,k are approximated again using the short second order central discretization:

(Q ∂u ∂x ) i+1/2, j,k = Q i+1/2, j,k (u i+1, j,k -u i, j,k ) h x + O(h 2 x ). (2.23) and (Q ∂u ∂x ) i-1/2, j,k = Q i-1/2, j,k (u i, j,k -u i-1, j,k ) h x + O(h 2 x ). (2.24)
where 2.22, 2.23 and 2.24 reduces to (also see Fig. 2.4):

Q i±1/2, j,k = (Q i, j,k + Q i±1, j,k )/2. Combining Equation
∂ ∂x (Q ∂u ∂x ) = Q i+1/2, j,k u i+1, j,k -(Q i+1/2, j,k + Q i-1/2, j,k )u i, j,k + Q i-1/2, j,k u i-1, j,k h 2 x + O(h 2 x ). (2.25)
For the second kind of item ∂ ∂y (Q ∂u ∂x ), a long second order central discretization is used:

( ∂ ∂y (Q ∂u ∂x )) i, j,k = (Q ∂u ∂x ) i, j+1,k -(Q ∂u ∂x ) i, j-1,k h y + O(h 2 y ). (2.26) in which (Q ∂u ∂x ) i, j+1,k and (Q ∂u ∂x ) i, j-1,k are discretized as: (Q ∂u ∂x ) i, j+1,k = Q i, j+1,k (u i+1, j+1,k -u i-1, j+1,k ) 2h x + O(h 2 x ). (2.27) and (Q ∂u ∂x ) i, j-1,k = Q i, j-1,k (u i+1, j-1,k -u i-1, j-1,k ) 2h x + O(h 2 x ).
(2.28)

A combination of Equation 2.26, 2.27 and 2.28 gives (also see Fig. 2.5):

∂ ∂y (Q ∂u ∂x ) = Q i, j+1,k (u i+1, j+1,k -u i-1, j+1,k ) -Q i, j-1,k (u i+1, j-1,k -u i-1, j-1,k ) 4h x h y +O[(h 2 x ), (h 2 y )]. (2.29)
The core idea of the finite difference discretization is to replace the derivatives like ∂ ∂x (Q ∂u ∂x ) and ∂ ∂y (Q ∂u ∂x ) with the discrete expressions omitting the truncation errors. Thus, a second order approximation for the above two derivatives is:

∂ ∂x (Q ∂u ∂x ) ≈ Q i+1/2, j,k u i+1, j,k -(Q i+1/2, j,k +Q i-1/2, j,k )u i, j,k +Q i-1/2, j,k u i-1, j,k h 2 x ∂ ∂y (Q ∂u ∂x ) ≈ Q i, j+1,k (u i+1, j+1,k -u i-1, j+1,k )-Q i, j-1,k (u i+1, j-1,k -u i-1, j-1,k ) 4h x h y        (2.30)
Such that, each derivative in Equation 2.20 can be replaced by the finite difference approximation.

MG methods

The above discretization results for each equation in Eqs 2.20 can be expressed in a simple matrix form, for example for the first equation:

L h u h = f h . (2.31)
where L h denotes a matrix on grid size h, u h denotes the unknowns and f h represents the right hand side term. If the size of L is small, then inversion is a direct choice: u h = L -1 f h . While for the large matrix L, which is the common case, the direct solver will be a disaster concerning computational time and required memory to store the matrix. An alternative is the iterative method. Iterative solvers like Jacobi type iteration or Gauss-Seidel type iteration usually are seen as "local" error smoothness, where "local" means the wavelength of the error is comparable to the local grid size, and they will become ineffective dealing with the relatively longer wavelength errors (an example is given in Section 2.3.2). Hence it appears a natural idea that the different wavelength errors should be treated separately (Fig. 2.6 shows a 1D example), which brings out the core of the MG methods.

The core of the MG methods is this: apply iterative solvers on grid Ω h , solve the local error comparable to the current grid size h, when convergence slows down, transfer error υ h to coarser grid Ω H , solve the approximated error equation L H υ H = f H , then when this problem is sufficiently solved, transfer this corrected error υ H back to the fine grid Ω h , and thus υ h is renewed and expected to have eliminated the long wavelength part.

Although the above description is straightforward, the real situation is more complicated and there are lots of issues that need to be specified. For example, the grids how to fit the real structure and how to choose the coarser grids, how to decide the operator L H on coarser grid, how to define the connection between two different levels of grid, how to select a proper relaxation scheme and etc. So in the following parts of this chapter, we will review some general treatments to specify the above mentioned details.

Grids

In continuum elastic mechanics, the computational domain is continuous and may have complicated boundaries. Usually, the first step in a numerical simulation is to map the target domain with suitable grids. If the boundaries are simple, then structure of grids are the natural and convenient choice (see Fig 2 .7). While for the domain with complicated boundaries, unstructure of grids are used (see a car example Fig 2 .3.1). In MG, unstructure of grids are difficult due to the construction of coarser grids. An alternative is to to use Algebraic MG (AMG) to generate coarser grids automatically [START_REF] Brandt | Algebraic multigrid theory: The symmetric case[END_REF][START_REF] Jw Ruge | Algebraic multigrid[END_REF].

Back to the objects studied in this thesis, since it is natural and reasonable to extract a cubic domain from a specimen, the structure of grids are the first choice. After generating the grids on the finest level Ω k , a sequence of coarser grids needs to be defined. Assuming the grid size on the finest level is

h k i (h 1 = h x , h 2 = h y , h 3 = h z ), on the next coarser level Ω k-1 it is h k-1
i , usually the ratio of mesh size h k i /h k-1 i is 0.5 since this ratio is the most efficient [START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF].

After defining a hierarchy of grids, we have decomposed the problem into different levels. The next step is to select a proper relaxation scheme on the highest level.

Smoothness solvers -iteration procedures

As mentioned above, the relaxation scheme on the finest grid is to smooth the local errors using iterative procedures like Jacobi type iteration or Gauss-Seidel type iteration. A typical iteration process is like this: give an initial approximation to every unknown, then using a prescribed order renew every point's approximation according to the residuals computed. The difference between Jacobi type relaxation and Gauss-Seidel type relaxation is the timing of the changes. In Jacobi type iteration, one obtains the new approximation simultaneously while in Gauss-Seidel type iteration, one obtains the new approximation by using the latest renewed values from surrounding points. Usually, Gauss-Seidel iteration has a better converged performance than Jacobi iteration and is selected as the smoothness solver in this thesis [START_REF] Trottenberg | Multigrid. Academic press[END_REF].

Adopting the simplified form in Equation 2.31, the Gauss-Seidel type relaxation scheme for the first equation in the Lamé equation 2.20 can be written as:

ūh i, j,k = ũh i, j,k + ωδ h i, j,k . (2.32) with δ h i, j,k = - f h i, j,k -L h ũh i, j,k 2 * (λ+2 * µ) h x 2 + 2 * µ h y 2 + 2 * µ h z 2 . (2.33) (a) (b) (c) (d)
Figure 2.9: Influence of lexicographic Gauss-Seidel relaxation on the error: a.error of initial guess (scale: 10 1 ), b.error after 10 relaxations on finest grid (scale: 10 0 ), c.error after the interpolation (scale: 10 -1 ), d.error after the 10 post relaxations on finest grid (scale: 10 -2 ).

and ω is the relaxation factor (0 ≤ ω ≤ 1), ūh i, j,k denotes the solution after relaxation at grid point (i, j, k) on level Ω h and ũh i, j,k denotes the solution before relaxation. The performance of the Gauss-Seidel relaxation depends on the relaxation order and the relaxation factor ω, which are discussed in detail in Section 2.5.

Here the relaxation order refers to the pointwise (lexicographic) or linewise relaxation order.

Transfer operators: restriction and interpolation

The connection between fine grid Ω h and coarse grid Ω H is through the transfer operators: the restriction operator I H h and the interpolation operator I h H . Once the convergence speed slows down on the fine grid, the approximated solution ũh and the residual r h (r h = f h -L h ũh ) are restricted through the operator I H h to the coarser grid Ω H . For example for the residual r, this process through a bi-linear restriction operator is:

r H i, j,k = (8 × r h 2 * i,2 * j,2 * k + 4 × (r h 2 * i-1,2 * j,2 * k + r h 2 * i+1,2 * j,2 * k + r h 2 * i,2 * j-1,2 * k + 1234r h 2 * i,2 * j+1,2 * k + r h 2 * i,2 * j,2 * k-1 + r h 2 * i,2 * j,2 * k+1 +)+ 2 × (r h 2 * i-1,2 * j-1,2 * k + r h 2 * i-1,2 * j+1,2 * k + r h 2 * i-1,2 * j,2 * k-1 + 1234r h 2 * i-1,2 * j,2 * k+1 + r h 2 * i+1,2 * j-1,2 * k + r h 2 * i+1,2 * j+1,2 * k + 1234r h 2 * i+1,2 * j,2 * k-1 + r h 2 * i+1,2 * j,2 * k+1 + r h 2 * i,2 * j-1,2 * k-1 + 1234r h 2 * i,2 * j-1,2 * k+1 + r h 2 * i,2 * j+1,2 * k-1 + r h 2 * i,2 * j+1,2 * k+1 )+ 1 × (r h 2 * i-1,2 * j-1,2 * k-1 + r h 2 * i-1,2 * j-1,2 * k+1 + 1234r h 2 * i-1,2 * j+1,2 * k-1 + r h 2 * i-1,2 * j+1,2 * k+1 + 1234r h 2 * i+1,2 * j-1,2 * k-1 + r h 2 * i+1,2 * j-1,2 * k+1 + 1234r h 2 * i+1,2 * j+1,2 * k-1 + r h 2 * i+1,2 * j+1,2 * k+1 ))/64.
(2.34)

One can find that restriction operator I H h combines 27 surrounding points for the central point (i, j, k) and is illustrated in Fig 2 .10.

After solving the approximated equation on a coarser grid Ω H , one expects that the long wavelength part of errors has been eliminated sufficiently and this corrected error should be transferred back to the approximated solution on fine grid Ω h to have a better approximation. This process is performed by the interpolation operator I h H and is the transpose of the restriction operator I H h :

I H h = ( h H ) d (I h H ) T . (2.35)
where d is the dimension of the problem and is 3 here, h H denotes the ratio of grid size between two adjacent levels.

Coarse grid operator

After defining the difference operator L h on the fine grid Ω h , the connection operator I H h and I h H , the issue remaining is to define a proper difference operator L H (also called coarse grid operator) on the coarser grid Ω H . A natural choice is to use the same difference operator as L h obtained on the fine grid Ω h . For some cases in this thesis, the direct difference operator exhibits good convergence performance while for some special cases, like the strongly discontinuous cases encountered later, it shows an unacceptable convergence speed or even diverges. Hence, for such cases, a more accurate coarse grid operator called Galerkin coarse grid operator should be used, which is:

L H = I H h L h I h H . (2.36)
This operator produces a 27-point stencil instead of the 7-point stencil in the direct difference operator. The disadvantages of the Galerkin coarse grid operator compared to the direct difference operator are the more complicated implementation and the heavier storage, which is discussed in detail in Section 2.7.3.

Chapter 2. Mechanics and numerical methods

Ω h : Ω H : pre-smoothing: L h ũh = f h residuals: r h = f h -L h ũh approximation corrected: ūh = ũh + υh post-smoothing: L h ūh = f h restriction: r H = I H h r h ûH = I H h ũh interpolation: υh = I h H υH smoothing: L H ũH = r H + L H ûH error corrected: υH = ũH -ûH Figure 2
.11: Two level MG scheme: a V cycle with the Full Approximation Scheme (FAS).

A typical two-grid cycle

By now, every part of a typical MG cycle has been illustrated. In this section, a two-grid cycle is demonstrated which combines all these parts (see Fig 2.11), that is:

(1) Pre-smoothing: on fine grid Ω h , relax v 1 times, obtain an approximation ũh to solution u h .

(2) Restricting: on fine grid Ω h , compute residuals r h = f h -L h ũh , transfer r h and ũh to coarse grid Ω H through the restriction operator I H h , obtain r H = I H h r h and ûH = I H h ũh on coarse grid Ω H .

(3) Smoothing: on coarse grid Ω H , relax v 0 times:

L H ũH = r H + L H ûH , obtain
a new approximation ũH .

(4) Interpolating: on coarse grid Ω H , compute the corrected error: υH = ũH -ûH , transfer υH to fine grid Ω h through the interpolation operator I h H : υh = I h H υH , obtain the corrected approximation: ūh = ũh + υh .

(5) Post-smoothing: on fine grid Ω h , relax v 2 times.

Pre-smoothing in

Step (1) is designed to eliminate the high frequency errors from the initial guess, which is illustrated in

Boundary conditions

In general, the boundary condition can be divided into 3 types: Dirichlet Boundary Condition -displacements specified on ∂Ω, Neumann Boundary Condition -stress specified on ∂Ω or the third Boundary Condition -displacements specified on a portion ∂Ω 1 and stress specified on a portion ∂Ω 2 . According to the problems encountered in the thesis, we outline 4 types of more specific boundary conditions and the related treatments in MG are illustrated. For convenience, all treatments are illustrated in 2D. However, the idea is as same as in 3D.

Homogeneous Strain Boundary Condition (HSBC)

In this case, the body is assumed to be subjected to a homogeneous strain field ε and produces the given strain on all boundaries. Then the displacements u on the boundaries can be specified through:

u i = εij x j , i, j = 1, 2, 3 (2.37)
where x j is the coordinate position. Hence HSBC [START_REF] Hashin | The elastic moduli of fiber-reinforced materials[END_REF] is a Dirichlet Boundary Condition. This boundary condition is often used to compute the elastic constants of a representative volume element (RVE) in composites, if the RVE is not periodic.

The treatments in MG are direct (see Fig 2.13): all the displacements can be specified on the boundaries (white circle) so that the relaxation scheme starts from the inner points (black circle) with the correct values on the boundaries, and unknowns on next coarse grid are still among the black ones.

Periodic Boundary Condition (PBC)

In composite materials, the same structure may be repeated thousands of times in different directions. It is a disaster for time and memory if simulating the entire structure. Therefore, it is natural to extract a RVE and thus a proper boundary condition is required. PBC [THKK00, XZE03] (see Fig 2 .15) is such a boundary condition that takes account of the periodic structure and only constrains the relations between boundaries according to a given strain field ε:

u i = εij x j + u * i , i, j = 1, 2, 3 (2.38)
where u * i is the unknown periodic displacement on the boundary ∂Ω and is dependent on the given global strain field εij and the structure. Note that the specific displacements u i can not be directly decided if only knowing εij without considering the structure.

The treatment for PBC in MG is a little complicated (see starts from the points on the boundary. The missing parts (dashed circles) of the difference operator for boundary points can be replaced by the corresponding parts inside the domain: note that the material properties (λ and µ) can be replaced directly by the corresponding ones while the unknowns (displacements) should consider the global strain ε before replacing (e.g. imaged point (-1,-1) connects to the inner point (7,7) and imaged point (9,7) connects to the inner point (1,7), and so does to the process of restricting, relaxing on coarse grid and interpolating,

Free Surface Boundary Condition (FSBC)

Under this condition, the body's boundary surfaces are supposed to be able to deform freely, which leads that the forces on these boundary surfaces are zero including the normal one and shear one. So it is a Neumann boundary condition. The displacements on the boundary surfaces are not known prior and solved through the stress equation 2.15 instead of the equilibrium equation 2.18:

σ i j = 0.
The treatment in MG for FSBC is a little different from the displacement boundary condition. First, the discretization form takes the second order upstream (backstream) approximation or long second order central approximation for different places, for example:

∂u ∂x =                -3u i +4u i+1 -u i+2 2h x i = 0 u i+1 -u i-1 2h x 0 < i < n x 3u i -4u i-1 +u i-2 2h x i = n x .
(2.39)

Then special attention should be taken when restricting or interpolating: for example in 2D, the bi-linear restriction operator I H h has:

for inner points (black circles in Fig 2 .16):

I H h =             1 2 1 2 4 2 1 2 1             (2.40)
for points on edges excluding vertices (e.g. red box in Fig 2 .16):

I H h =             2 2 4 4 2 2             (2.41)
for vertices (e.g. blue box in Fig 2 .16):

I H h =       4 4 4 4       (2.42)
At last, it is better to have several local relaxations for the neighboring points to the boundaries since the relaxation of inner points are sensitive to the accuracy of these boundary points.

Symmetric Boundary Condition (SBC)

For this condition, the body simulated is assumed to be one part of a symmetry one. On the symmetric plane, the displacements perpendicular to this plane are zero and the shear stresses on this plane is zero as well. Hence SBC involves both the displacement and stress boundary condition at the same time:

u i = 0, u i, j = 0.
(2.43)

The treatment in MG is the combination of HSBC and FSBC.

Local mode analysis

Local mode analysis in MG, which is called the Local Fourier Analysis (LFA), was first introduced by A.Brandt [START_REF] Brandt | Multi-level adaptive solutions to boundary-value problems[END_REF] in 1970s and later extended and refined [START_REF] Brandt | Multigrid solvers for non-elliptic and singular-perturbation steady-state problems[END_REF][START_REF] Brandt | Rigorous quantitative analysis of multigrid, i. constant coefficients two-level cycle with l_2-norm[END_REF]. The nature of LFA is that, under general assumptions, any general discrete operator, nonlinear, with nonconstant coefficients, can be linearized locally and can be replaced locally (by freezing the coefficients) by an operator with constant coefficients. Thus, general linear discrete operators with constant coefficients are considered in LFA. Formally, they are defined on an infinite grid [START_REF] Trottenberg | Multigrid. Academic press[END_REF].

In this section, we use local mode analysis to quantify the performance of MG methods in solving the Lamé system of equations 2.18. For convenience to demonstrate and analyze, we select the first equation in Equation 2.20 and discretize it using the Finite Difference Method like in Equations 2.30. The grids are assumed to be uniform rectangular which means h x = h y = h z , and material is homogeneous, where λ and µ are assigned to be constant. So the first equation in Eq 2.20 can be written as:

(2 + λ µ ) ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 + ∂ 2 u ∂z 2 = f. (2.44) with f = 0 - λ µ ∂ 2 v ∂x∂y - ∂ 2 v ∂x∂y - λ µ ∂ 2 w ∂x∂z - ∂ 2 w ∂x∂z . (2.45)
With Gauss-Seidel relaxation, we have:

ūh i, j,k = ũh i, j,k + ωδ h i, j,k . (2.46) with δ h i, j,k = - 1 8 + 2 λ µ [ f h i, j,k h 2 -(2 + λ µ )(ū h i-1, j,k -2ũ h i, j,k + ũh i+1, j,k )- (ū h i, j-1,k -2ũ h i, j,k + ũh i, j+1,k )- (ū h i, j,k-1 -2ũ h i, j,k + ũh i, j,k+1 )].
(2.47)

The error before the relaxation is defined by: υh i, j,k = u h i, j,kũh i, j,k .

(2.48) and the error after the relaxation is defined by: ῡh i, j,k = u h i, j,kūh i, j,k .

(2.49) Substituting Equation 2.48 and 2.49 in 2.46 and 2.47, we obtain:

ῡh i, j,k = u h i, j,k -(ũ h i, j,k + ωδ h i, j,k ) = υh i, j,k + ω 8 + 2 λ µ [(2 + λ µ )( ῡh i-1, j,k -2 υh i, j,k + υh i+1, j,k )+ 12345678910111213141( ῡh i, j-1,k -2 υh i, j,k + υh i, j+1,k )+ 12345678910111213141( ῡh i, j,k-1 -2 υh i, j,k + υh i, j,k+1 )] = (1 -ω) υh i, j,k + ω 8 + 2 λ µ [(2 + λ µ )( ῡh i-1, j,k + υh i+1, j,k )+ 12345678910111213141516171( ῡh i, j-1,k + υh i, j+1,k )+ 12345678910111213141516171( ῡh i, j,k-1 + υh i, j,k+1 )].
(2.50)

From Equation 2.50, one can find that the error after relaxation for a given point is affected by the error before relaxation from itself and the error before and after relaxation from its neighbors. Then for a given point i, j, k we expand Equation 2.50 in Fourier series:

ῡh i, j,k = 0<|θ m |≤π Ā(θ m )e ιθ m m m = 1, 2, 3 (2.51)
Similarly the error before relaxation can be expanded to:

υh i, j,k = 0<|θ m |≤π Ã(θ m )e ιθ m m m = 1, 2, 3 (2.52)
Substitute Equation 2.51 and 2.52 in 2.50 and eliminate the same terms, we get: (2.53)

Ā(θ) =(1 -ω) Ã(θ) + ω 8 + 2 λ µ [(2 + λ µ )( Ā(θ)e -ιθ 1 + Ã(θ)e ιθ
So the amplitude amplification factor of the relaxation for component θ is

ρ(θ) = Ā(θ) Ã(θ) = (1 -ω)(8 + 2 λ µ ) + ω[(2 + λ µ )e ιθ 1 + e ιθ 2 + e ιθ 3 ] 8 + 2 λ µ -ω[(2 + λ µ )e -ιθ 1 + e -ιθ 2 + e -ιθ 3 ]
.

(2.54)

If the Poisson ratio ν is a constant (0.3), then λ µ = 2ν 1-2ν is a constant ( λ µ = 1.5) as well. Here one assumes that the low frequency component (|θ i | ≤ π 2 ) can be eliminated sufficiently while the high frequency ones cannot be solved on the next coarser grid, which gives:

ρ = min(ρ(θ i )) π 2 ≤ θ i ≤ π. (2.55)
For pointwise relaxation, ρ = 0.704 when ω = 1. Therefore, pointwise relaxation is already a good choice. For example, one can obtain an elimination speed of 2.87 per V(2, 1) cycle. We observed this speed in our code in practice.

Local grid refinement

Local grid refinement (also called adaptive grid refinement) (see Fig 2.17) is an efficient tool in MG when a locally much finer grid is required. Compared to a model with a fine grid everywhere, the computing time and allocated memory are reduced dramatically in a locally refined model. A natural way to regard the local grid is to treat it as one part of the fine grid (see Fig 2 .18): we start from the hierarchy of global grids but on fine grid we focus only on certain domains, and relax only on these domains. Therefore, sub boundary conditions will be introduced (see Fig 2 .19) and can be seen as a "Dirichlet Boundary Condition" so that the relaxation step are carried out only on the sub domains excluding points on sub boudaries, so does one to compute residuals.

Strong discontinues -Jumping coefficients

It is natural for composite materials to have material properties that differ greatly in different parts, which leads to strongly discontinuous coefficients λ and µ. Such discontinuities can lead to slow convergence or even divergence when using standard MG methods. That is because, if λ or µ jump significantly, so does u i,i , while Figure 2.17: Local grid refinement.

Ω 1 Ω 2 Ω 2+1 Ω 2+2
Figure 2.18: Hierarchical grid. in standard MG, only smooth functions can be approximated properly on coarse grids, which means that the quantity transferred through the standard bi-linear restriction and interpolation operator (see Section 2.7.2) is not smooth on coarse grids. Therefore, it is the smoothed correction (λ + 2µ)u i,i or µu i,i that should be transferred instead of u i,i . Based on this, there are some special treatments for this kind of problem in MG methods.

Ω 1 Ω 2 Ω 2+1 Ω 2+2

Coefficients on coarse level

For convenient illustration, a uniform 3D cubic mesh is assumed, giving h = h x = h y = h z . According to the discretization forms in Eq. 2.30, coefficients Q i±1/2 locate between two mesh points (see the left figure in Fig. 2.20, and stand for the mesh points, stands for the coefficient location). In order to reflect the connections from fine level to coarser level as accurately as possible, these coefficients on the coarser level should be defined through the ones on fine level. Alcouffe et al. [START_REF] Re Alcouffe | The multigrid method for the diffusion equation with strongly discontinuous coefficients[END_REF] proposed several alternatives for the 2D case and based on that Hoekema et al. [START_REF] Hoekema | Multigrid solution of the potential field in modeling electrical nerve stimulation[END_REF] proposed a definition for the 3D case, which is adopted in this thesis. Fig. 2.20 shows an example to illustrate how the coefficients on the coarse level are constructed through the ones on its next finer level. The left figure in Fig. 2.20 shows mesh points on the finer level k, where are the points coinciding the ones on the coarse level k -1 and are the left mesh points on the finer level k. We select points A and B as the example to illustrate. Points A and B appear on the coarse level and the coefficient ( , see the right figure in Fig. 2.20) between A and B locates at the center. As shown in Fig. 2.20, the link between A and B is the boundary line between red rectangular loop (HGFBCDEA) and green rectangular loop (NMLBIJKA). Hence, on the coarse level is defined as the weighted average of , and on the finer level. If we use Q k-1 AB representing the coefficient between A and B on the coarse level k -1 and Q k OB representing the coefficient between O and B on the finer level k, then Q k-1 AB is defined as:

Level k Level k -1 O A B C D E F G H I J K L M N A B
Q k-1 AB = 1 2 δ(Q k AO + Q k OB ) + 1 8 [δ(Q k BC + Q k CD + Q k DE + Q k EA )+ δ(Q k BF + Q k FG + Q k GH + Q k HA )+ δ(Q k BI + Q k I J + Q k JK + Q k KA )+ δ(Q k BL + Q k LM + Q k MN + Q k NA )].
(2.56)

where

δ(Q k AO + Q k OB ) = 1 1 2 ( 1 Q k AO + 1 Q k OB ) .
(2.57) and

δ(Q k BC + Q k CD + Q k DE + Q k EA ) = 1 1 4 ( 1 Q k BC + 1 Q k CD + 1 Q k DE + 1 Q k EA ) .
(2.58)

Intergrid operators

As mentioned above, for strong discontinuous case, it is Q i±1/2 u i instead of u i that is smoothed and transferred from the fine level to the coarse level. We explain how to implement this process through the interpolation since restriction can be seen as the transpose of interpolation.

The interpolation process is divided into 4 steps [START_REF] Re Alcouffe | The multigrid method for the diffusion equation with strongly discontinuous coefficients[END_REF], as shown in Fig. 2.21.

(1) In the first step, the values of u on the coarse level k -1 points are directly given to the points on the next finer level k coinciding with the ones on the coarse level k -1, e.g (a) in Fig. 2.21:

u k O = u k-1 O .
(2.59)

(2) After the loop in Step 1, all the points on the fine level k that coincide with the ones on the coarser level k -1 have been given the value. The second step is to obtain the values in the points locating at the center of the link between two coarse level points through the bilinear interpolation, e.g (b) in Fig. 2.21:

u k O = Q k OA u k A + Q k OB u k B Q k OA + Q k OB . (2.60)
(3) After Step 2, the values in the points at the center of the "edge" are defined.

Then the values in the points at the center of a side of the "cube" are defined through the ones at the center of the "edge", e.g (c) in Fig. 2.21:

u k O = Q k OA u k A + Q k OB u k B + Q k OC u k C + Q k OD u k D Q k OA + Q k OB + Q k OC + Q k OD .
(2.61) (4) At last, in the cubic element, only the value at the center of the "cube" remains to be defined, which is obtained through its surrounding six points at the center of the surfaces of the "cube", e.g (d) in Fig. 2.21:

u k O = Q k OA u k A + Q k OB u k B + Q k OC u k C + Q k OD u k D + Q k OE u k E + Q k OF u k F Q k OA + Q k OB + Q k OC + Q k OD + Q k OE + Q k OF . (2.62)
It is obvious to see that if Q is constant, then the interpolation can be seen as the standard linear interpolation at Section 2.3.3.

The Galerkin coarse grid operator

For the strongly discontinuous case, as mentioned in Section 2.3.4, a more proper choice for the coarse grid operator is to use a Galerkin like operator (see Equation 2.36). The reason is direct: the problem with strong discontinuities on fine grid Ω h should be reflected well enough on the coarse grid Ω H , where the Galerkin operator keeps this property through a natural way.

Compared to the direct difference operator, the implementation of Galerkin operator is more complicated. For a given grid point (i, j, k) on coarse grid, the first step is to interpolate its value to its 27 surrounding points on fine grid through operator I h H . The next step is to compute the value of the points around these 27 points through L h I h H , which leads that for a 7-point fine grid operator , there are 73 points affected by these 27 points from interpolation through L h while for a 27-point one, there are 125 points in total affected. The last step is to restrict value on these affected points in second step to 27 points on coarse grid around the given point (i, j, k). In this context, L H can be unsymmetrical.

The disadvantages are: 1. time to generate the coarse grid operator. 2. memory to store the 27-point stencil. 3. In nonlinear problems, the coefficients of the coarse operator need to be recalculated after every linearization step.

Conclusions

In this chapter, we have reviewed some basic notations in both mechanics and MG methods. A MG algorithm was built for the Lamé equation. We discussed several kinds of boundary conditions suitable for real applications. With the local mode analysis, we are able to predict the optimal convergence factor for the Lamé equation. At last the special treatments in MG for strongly discontinues case are illustrated.

Chapter 3

Composite homogenization

Introduction

The global elastic behavior of composite materials can be predicted from the microstructure using homogenization methods. The purpose of this chapter is to validate the MultiGrid solutions both in micro-scale and macro-scale and present the application of MultiGrid methods in composite material homogenization. Firstly, in Section 3.2, combining all the techniques mentioned above, we deal with a simple but fundamental problem in composite mechanics: inhomogeneity within a halfspace, where the efficiency and accuracy of MG methods for the Lamé equation are examined. Secondly, several analytical homogenization methods available in the literature are reviewed in Section 3.3. Thirdly, a numerical homogenization model using MultiGrid is investigated in Section 3.4. The comparison between the analytical and numerical predictions are made for a simple structure case. At last, a more complex but common structure -woven material, is studied using MultiGrid homogenization.

Microscopic -Eshelby's solution

Equivalent inclusion method

In elastic mechanics, Eshelby's solution refers to the analytical solution of a problem involving elliptical inclusions or inhomogeneities within an infinite homogeneous body [Esh57, Esh59, Esh61], which is the basic model in composite structure. According to the definition of Eshelby, an inclusion describes a domain that its elastic properties are as same as the matrix but it is subjected to the so-called eigenstrains while an inhomogeneity refers to a domain that it has the different elastic properties from its surroundings. In this thesis, we consider the problem with inhomogeneities.

The local stress field inside the inhomogeneity can be obtained from Eshelby's "equivalent inclusion method", which is simulated by imaging a fictitious eigenstrain ε * in inhomogeneity. One considers a matrix Ω with stiffness tensor C 1 subject to a remote uniform strain ε ∞ and containing a elliptical inhomogeneity D with stiffness tensor C 2 . The stress field inside the inhomogeneity is defined as:

σ ∞ + σ = C 2 (ε ∞ + ε) = C 1 (ε ∞ + ε -ε * ) in D. (3.1)
where σ ∞ is the uniform stress obtained if there is no inhomogeneity inside the matrix, which obeys the Hooke's law: σ ∞ = C 1 ε ∞ . ε and σ are the additional disturbed strain and stress respectively due to the existence of the inhomogeneity. Eigenstrain ε * is related to the additional disturbed strain ε:

ε = Sε * (3.2)
S is the Eshelby tensor, which is determined by the shape of the inhomogeneity. The details of the determination of the Eshelby tensor components S i jkl can be found in [START_REF] Mura | Micromechanics of defects in solids[END_REF] . Thus, the analytic local stress field inside the inhomogeneity can be determined through equations 3.1 and 3.2.

Validation with MG solution

To validate the accuracy of MG methods in solving this fundamental problem, we assume an approximated model: a hard spherical inhomogeneity is embedded at the center of a relative larger matrix (see Fig. 3.1), where the ratio between the sphere radius r and the cube length L is r/L = 1/20, the ratio of Young's modulus between the inhomogeneity (E i ) and the matrix (E m ) is defined by E i /E m = 10 and the Poisson's ratio is set same ν = 0.3. A uniform traction strain ε ∞ zz = 0.01 is applied on the boundaries of the matrix. The local hierarchical grid refinement is adopted to obtain the high resolution: 409 points are mapped along a sphere diameter. Eshelby's solution indicates that the local stress field produced inside the inhomogeneity is uniform. A very good agreement is found in Fig 3 .2 (except for the adjacent points near the boundary between sphere and matrix, where these points are coupled with the points in the matrix), which gives a comparison between Eshelby's analytical solution and the MG solution. If we compare them more precisely e.g. at the center z = 0:

(σ ∞ + σ)/σ ∞ zz (σ ∞ xx + σ xx )/σ ∞ zz (σ ∞ yy + σ yy )/σ ∞ zz (σ ∞ zz + σ zz )/σ ∞ zz σ eshelby /σ ∞ zz
σ Eshelby xx = 0.582σ ∞ zz , σ MG xx = σ ∞ xx + σ xx = 0.584σ ∞ zz σ Eshelby yy = 0.582σ ∞ zz , σ MG yy = σ ∞ yy + σ yy = 0.584σ ∞ zz σ Eshelby zz = 1.662σ ∞ zz , σ MG zz = σ ∞ zz + σ zz = 1.665σ ∞ zz
where the maximum difference is 0.34%.

Analytical homogenization methods for composite materials

The microstructure determines the behavior of composite materials. However, it is unacceptable to take account of the full microstructure when analyzing a large scale object since the cost will be overwhelming both in time and memory. Therefore, a lot of effort has been payed by researchers to develop approximated models on a representative region of a composite material, which is usually called the Representative Volume Element (RVE). Researchers are sought to predict the overall behavior of the RVE, which usually refers to the homogenization process. If a RVE exists for a random composite material, Sab [START_REF] Sab | On the homogenization and the simulation of random materials[END_REF] has shown that the homogenized material properties of the RVE can represent the overall properties of the material.

The Eshelby's equivalent inclusion method (see Section 3.2) gives the analytical solution of the stress field in a heterogeneous material containing simple shape inhomogeneities (ellipsoid, cuboid, ...). Various rigorous methods for predicting the overall material properties have been developed according to the Eshelby's theory. These methods take account of an important parameter-the volume fraction, which denotes the ratio of volume between the inhomogeneity and the RVE and directly determines the effective material properties. These methods are reviewed in the following sections and will be compared with results obtained from the MG calculations.

Hashin-Shtrikman bounds (HS)

Hashin and Shtrikman [START_REF] Hashin | On some variational principles in anisotropic and nonhomogeneous elasticity[END_REF][START_REF] Zshtr Hashin | A variational approach to the theory of the elastic behaviour of polycrystals[END_REF][START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] proposed an approximate equation to estimate the upper and lower bounds for the effective stiffness tensor by assuming that the inhomogeneity is spherical and bonds perfectly to the matrix. The upper and lower bounds of the effective stiffness tensor are determined by:

C e f f (HS +) = C 2 + V 1 [(C 1 -C 2 ) -1 + V 2 S 2 C -1 2 ] -1 C e f f (HS -) = C 1 + V 2 [(C 2 -C 1 ) -1 + V 1 S 1 C -1 1 ] -1 . (3.3)
where V 1 is the volume fraction of matrix and V 2 is the volume fraction of inhomogeneity and V 1 + V 2 = 1, S 1 and S 2 are the Eshelby tensor of matrix and inhomogeneity respectively and C 1 and C 2 denote the stiffness tensor of the matrix and the inhomogeneity. This formula provides a good prediction if the ratio of elastic properties between matrix and inhomogeneity is not too large.

Mori-Tanaka method (MT)

The Mori-Tanaka method [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF] assumes that each inhomogeneity is embedded in an infinite matrix subject to a remote strain ε ∞ and that the strain produced inside the inhomogeneity ε (2) is connected to strain in the matrix ε (1) and the remote strain ε ∞ through the influence tensor T and the concentration tensor A respectively, which are defined as follows:

ε (2) = Tε (1) .
(3.4)

52

Chapter 3. Composite homogenization and

ε (2) = Aε ∞ . (3.5)
where T can be obtained through the Eshelby tensor S 2 of the inhomogeneity, stiffness tensor of the matrix C 1 and the inhomogeneity C 2 :

T = [I + S 2 C -1 1 (C 2 -C 1 )] -1 . (3.6)
I is the unit tensor and A can be written in terms of T, which is:

A = T[V 1 I + V 2 T] -1 . (3.7)
Then the effective stiffness tensor C e f f (MT ) is given by [START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF]:

C e f f (MT ) = C 1 + V 2 (C 2 -C 1 )A. (3.8)
The Mori-Tanaka method gives a good estimation if the inhomogeneity volume fraction V 2 is less than 50% [START_REF] Ponte | The effect of spatial distribution on the effective behavior of composite materials and cracked media[END_REF]. While for the higher volume situation (larger than 50%), the Mori-Tanaka method considers the higher volume part as the matrix material. In fact, for harder inhomogeneities, the results obtained from the Mori-Tanaka method correspond to the lower bound of the Hashin-Shtrikman method when V 2 < 50% and to the upper bound when V 2 > 50%. While for the soft case, the Mori-Tanaka's prediction corresponds to the Hashin-Shtrikman's upper or lower bounds if V 2 is smaller or larger than 50%.

Self-consistent method (SC)

The self-consistent method was originally developed for polycrystals [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF][START_REF] Budiansky | On the elastic moduli of some heterogeneous materials[END_REF][START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] by assuming that each inhomogeneity is embedded in an infinite matrix with an effective stiffness tensor C e f f . Therefore, the computation of C e f f is performed by an iterative procedure:

(1) Initialize the effective stiffness tensor: Cef f

(S C) = C 1 (2) Compute the influence tensor: T = [I + S 2 ( Cef f (S C) ) -1 (C 2 -Cef f (S C) )] -1 .
(3) Compute the concentration tensor: (4) Compute renewed effective stiffness tensor:

A = T[V 1 I + V 2 T] -1 .
C e f f (S C) = Cef f (S C) + V 2 (C 2 - Ce f f (S C) )A.
(5) Go back to step 2 using the renewed C e f f , redo step 3 and 4 until the difference between Cef f and C e f f converges below a set error.

The convergence speed is fast if the material properties of the inhomogeneity and the matrix do not differ too much.

Numerical homogenization using MG methods

Analytical homogenization methods mentioned above give a good estimation for the simple microstructure case. However, they are restricted to specific geometry, ratio of material properties and number of phases in the composite material. Traditional numerical simulation methods like the Finite Element method overcome these difficulties and are widely used to determine the effective coefficients of composite materials. In this thesis, MG methods are adopted since they can simulate the complicated structure using dense grids and maintain a lower computing time and memory allocation compared with other numerical methods.

A typical numerical homogenization procedure for computing effective elastic coefficients of a RVE is straightforward:

(1) Apply proper boundaries conditions like homogeneous strain boundary condition or periodic boundary condition on the boundary surfaces of the RVE. This step is usually divided into several parts to model the pure tensile or shear loading situation respectively. In 3D, there are 6 loading cases (3 tensile and 3 shear, see Fig 3 .3).

(2) According to the given boundary conditions, one obtains the microstructural stresses σ i j after numerical calculation and computes the global average stress σij .

(3) Compute strain energy based on the average or microstructural stress and strain obtained, which enables us to determine the relation of global stressstrain response due to the Hill-Mandel energy condition [Man66, Hil66, Hil67].

Coupling between microscopic and macroscopic

A numerical homogenization procedure links the stress-strain response in macroscopic scale and microscopic scale, where the core lies on the equivalence of strain energy between the actual microstructural case and the averaged homogeneous case, where the strain energy stored in the microstructure is:

U = 1 2 V σ i j ε i j dV. (3.9)
and in the homogeneous case:

Ū = 1 2 σij εij V = 1 2 C i jkl εkl εij V. (3.10) where εij = 1 V V ε i j dV σij = 1 V V σ i j dV.
(3.11)

If U = Ū, the effective elastic constants C i jkl can be determined through:

C i jkl = 2U εkl εij V . (3.12)
The equivalence between U and Ū has been examined by many researchers [START_REF] Sun | Prediction of composite properties from a representative volume element[END_REF]. A rapid review is given here: Equation 3.9 can be rewritten as:

U = 1 2 V σ i j (ε i j -εij + εij )dV = 1 2 V σ i j (ε i j -εij )dV + 1 2 εij V σ i j dV = 1 2 V σ i j (u i, j -ūi, j )dV + Ū = 1 2 V [σ i j (u i -ūi )] , j dV - 1 2 V σ i j, j (u i -ūi )dV + Ū.
(3.13) Using the fact that σ i j, j = 0 according to the equilibrium equation 2.9 and converting volume integral into a surface integral by using Gauss theorem, we can reformulate Equation 3.13 as:

U -Ū = 1 2 S σ i j (u i -ūi )n j dS . (3.14)
where S is the surface and n is the unit outward normal. For homogeneous boundary condition, we have: In this section, we use MG methods to model these two boundary conditions when solving the problem of a spherical inhomogeneity embedded in a cubic matrix (see Fig. 3.1). The ratio between the radius (r) of the sphere and the length (L) of the cube is 0.3, which gives the inhomogeneity a volume fraction V 2 = 11.31%. The sphere is 10 times harder than the matrix. A total 257 3 grid points are used on the finest level. Fig 3 .5 shows how the RVE deforms under two different types of boundary conditions. It can be seen that the boundary surface under the shear loading in PBC does not remain in one plane. In fact, HSBC for the shear loading case is a over-constrained boundary condition since it requires that the plane remains plane, which indicates that the inhomogeneity has the same deformation as the matrix. This leads to the upper bounds of stiffness tensor compared to ones under PBC, which has been indicated by many researchers [START_REF] Needleman | Comparison of crystal plasticity and isotropic hardening predictions for metal-matrix composites[END_REF][START_REF] Sun | Prediction of composite properties from a representative volume element[END_REF].

u i = ūi on S . ( 3
Tabl 3.1 and 3.2 give the results of the analytical prediction and the numerical calculation of the stiffness matrix respectively. Although HSBC leads to the upper bounds, it stays within range of HS bounds. The stiffness matrix obtained under either PBC or HSBC resembles the Self-Consistent scheme most, where the biggest difference is 2.45% and 1.66% respectively. Therefore, both PBC and HSBC can (C e f f i jkl ) MG(PBC) = i jkl 11 22 33 12 23 13 11 4.6276 1.8865 1.8865 0.0000 0.0000 0.0000 22 1.8865 4.6276 1.8865 0.0000 0.0000 0.0000 33 1.8865 1.8865 4.6276 0.0000 0.0000 0.0000 12 0.0000 0.0000 0.0000 1.3215 0.0000 0.0000 23 0.0000 0.0000 0.0000 0.0000 1.3215 0.0000 13 0.0000 0.0000 0.0000 0.0000 0.0000 1.3215

(C e f f i jkl ) MG(HS BC) =
i jkl 11 22 33 12 23 13 11 4.6480 1.8855 1.8855 0.0000 0.0000 0.0000 22 1.8855 4.6480 1.8855 0.0000 0.0000 0.0000 33 1.8855 1.8855 4.6480 0.0000 0.0000 0.0000 12 0.0000 0.0000 0.0000 1.3623 0.0000 0.0000 23 0.0000 0.0000 0.0000 0.0000 1.3623 0.0000 13 0.0000 0.0000 0.0000 0.0000 0.0000 1.3623 be chosen as boundary condition in the homogenization process for this particular case.

Mesh size

MG methods provide the possibility of using dense grids. In this section, the influence of the grid size on the accuracy is discussed. The same model as in Section 3.4.2 is used but under PBC only. A series of grid points from 5 3 to 513 3 is applied (see Fig 3 .6, grid lines are not drawn in 257 3 and 513 3 cases since they are too dense). 3.3, where a constant is found. It is also found that even for the coarsest grid (5 3 ), the MG results are still in the bounds of Hashin-Shtrikman. With the increase of grid density, C MG 1111 converges to the Self-Consistent result and C MG 1122 and C MG 1212 approximate the lower bound of the Hashin-Shtrikman or Mori-Tanaka solution and converge to the value that exceeds a little the Hashin-Shtrikman's lower bound. Another finding is that, on level 4, where grid points are 33 3 and there are about 20 points along a sphere diameter, the MG results has approximated the solution on level 7 with an difference smaller than 2%.

Volume fraction

In this section, the influence of the inhomogeneity volume fraction on the MG solution is studied. The model is the same as the one used in Section 3.4.2 but with the radius changing from 0 to 0.5L, which gives a volume fraction from 0% to 52.36%. PBC is applied only and 257 3 grid points are used in each case. liable as pointed out in Section 3.3.2. These analytical homogenization methods produce a similar prediction when the inhomogeneity volume fraction is small. For bigger volume fractions, the difference between the analytical predictions becomes pronounced.

Ratio

In this section, the influence of the ratio of Young's modulus between the matrix and the sphere is also investigated. The radius of the sphere is fixed to be 0.3L as in Section 3.4.3 and PBC is applied. The density of the grid points are chosen to be 129 3 according to the observation in Section 3.4.3. The soft and hard spherical inhomogeneities are both simulated. 3.5 Homogenization: woven structure

Model description

In this section, a more complex but common structure-a woven material is simulated (see Fig 3 .12). The dimension of the block selected is [24r, 24r, 6r], where r is the radius of the fiber. These fibers along two different directions twine each other and are distributed uniformly with the distance d = 3r. The block is full filled with the matrix and the fibers with a fiber volume fraction V f = 48.19%. The fibers are 10 times harder than the matrix. PBC is applied. Three normal loading (ε xx , ε yy , ε zz ) and two shear loading (ε xy , ε xz ) are simulated since the structure is only symmetric in xoz and yoz plan.

Numerical results

The symmetry of the structure makes it possible to simulated on the smaller scale. A subdomain outlined with red lines in Fig 3 .12 is selected as the comparison with the entire domain. The dimension of the partial domain is [6r, 6r, 6r]. The density of the grid points is 385 × 385 × 97 and 97 × 97 × 97 for the entire and partial structures respectively. The same PBC is applied on both structures. Table 3.4 shows the effective stiffness matrix for two cases under PBC. A good agreement is found where the maximum difference is 0.54%. Different stiffness ratio cases are simulated and Fig 3 .14 shows the changing of each component. It is clear to find that the increasing of the stiffness ratio does not produces the same increase amplitude to each component, where those components along the fiber directions (x or y) increase most.

(C e f f i jkl ) MG(entire) = i jkl 11 22 33 12 23 13 11 13.1963 4.1635 3.2331 0.0000 0.0000 0.0000 22 4.1635 13.1963 3.2331 0.0000 0.0000 0.0000 33 3.2331 3.2331 8.6466 0.0000 0.0000 0.0000 12 0.0000 0.0000 0.0000 3.6623 0.0000 0.0000 23 0.0000 0.0000 0.0000 0.0000 2.5577 0.0000 13 0.0000 0.0000 0.0000 0.0000 0.0000 2.5577

(C e f f i jkl ) MG(partial) =
i jkl 11 22 33 12 23 13 11 13.1978 4.1411 3.2283 0.0000 0.0000 0.0000 22 4.1411 13.1978 3.2283 0.0000 0.0000 0.0000 33 3.2283 3.2283 8.6685 0.0000 0.0000 0.0000 12 0.0000 0.0000 0.0000 3.6472 0.0000 0.0000 23 0.0000 0.0000 0.0000 0.0000 2.5594 0.0000 13 0.0000 0.0000 0.0000 0.0000 0.0000 2.5594 The cost of the MG calculation maintains small. For such a model with 7.2 million grid points (197 3 ), it is solved in 2 hours.

Conclusions

This chapter discuss the MG solutions both in micro-scale and macro-scale. At first in micro-scale, the MG solutions are compared with the classical Eshelby solution, where a high accuracy is found. Then in macro-scale part, the MG solver is developed for the application of homogenization process in composites. In this part, we first reviews the notation of the homogenization and the related analytical prediction methods in the literature. Then the numerical homogenization method using MG is presented. The influence of the boundary conditions, the mesh size, the volume fraction and the stiffness ratio on the MG solution is investigated. The comparison of the accuracy with the analytical homogenization methods are made. At last, the MG methods are used to solve the woven structure problems. The potential of using MG methods in homogenization process is revealed, which overcomes the limitation of the analytical methods while keeping numerical cost limited.

Chapter 4

Modeling composite lamination:

Free edge effects

Introduction

Laminated composite materials play an important role in modern industry since they can provide the required mechanical properties by modifying the fiber orientations and the stacking sequence. The pre-peg unidirectional layers consisting of high-modulus fibers and glued matrix are assembled into a layered structure. However, due to the mismatch of the elastic properties of the adjacent layers, a well known phenomenon called the free edge effect occurs at the inter-laminar interface in the vicinity of the free edges under uniaxial tension [START_REF] Pipes | Interlaminar stresses in composite laminates under uniform axial extension[END_REF][START_REF] Nj Pagano And R Byron | The influence of stacking sequence on laminate strength[END_REF]. Stress In this chapter, the MG method is applied to simulate this phenomenon. Firstly, the free edge effect is discussed in detail in Section 4.2 and a comparison is made between the Finite Element homogenized model and the MG micro-structural model. Then the real structure, obtained using X-ray tomography, is introduced into the model in Section 4.3. At last, the parameters that may influence this phenomenon are investigated in Section 4.4.

Free edge effects

Let us assume a laminated structure made of several layers under an uniaxial extensional strain ε zz applied at the ends, see Fig 4 .1. The layers are consid-ered to be bonded perfectly. In each layer, fibers are distributed uniformly with the same orientation angle θ with respect to the z direction. θ is chosen to be 0 • ≤ θ ≤ 90 • and can differ from different layers. The Classical laminated Plate Theory (CLPT) [START_REF] Robert | Mechanics of composite materials[END_REF][START_REF] Hong | Introduction to composite materials[END_REF], which is generally used in composite design, assumes a two dimensional plate state and only takes into account the intralaminar stresses, where in this context it refers to the two normal in-plane stresses σ zz and σ yy and the one shear in-plane stress τ yz . However, many failures happen at a loading level that is much lower than the one predicted by the CLPT. In fact, due to the discontinuities of the elastic properties in adjacent layers, the complex interlaminar stresses, which are the two shear components τ xy and τ xz and the normal one σ xx , emerge at the interface between two adjacent layers. Although these stresses (τ xy , τ xz and σ xx ) exist inside the layers as well, they are usually largest at the layer interfaces [START_REF] Robert | Mechanics of composite materials[END_REF].

The existence of the interlaminar stresses can be derived from the equilibrium conditions. For example, consider that the structure in Fig 4 .1 is sufficiently long in the longitudinal z direction and subjected to the uniform axial extension. According to the Saint-Venant principle, each stress component is independent of the z axis and solely depends on the x and y directions, this leads the third equation in equilibrium equation 2.18 to be:

τ xz,x + τ yz,y = 0. (4.1)
Since the free surface boundary condition is applied on the boundary plane (x, z), the shear stress τ yz decays rapidly near the free edge. This leads to the variation of τ yz in y direction and thus τ xz should also have a variation in x direction to satisfy the balance equation (4.1). This sharp variation leads to the singularity in the shear stress distribution, which was firstly observed by the numerical pioneering work of Pipes and Pagano [START_REF] Pipes | Interlaminar stresses in composite laminates under uniform axial extension[END_REF] in 1970. Such a stress concentration occurs at the vicinity of the free edges in the interface between the two different orientation layers, which dominates the laminated strength and life and is named the free edge effect by researchers. However, no exact analytical solution is known yet, due to the complexity and an alternative is the numerical simulation.

Homogenization model

The simulation of the laminated structure is usually carried out on two scales: microscopic and macroscopic. At the microscopic level, all the details including the fibers, matrix, voids... are modeled simultaneously, which, however, becomes impractical due to the complexity. An alternative is to consider each layer to be homogeneous anisotropic at the macroscopic level through the homogenization method mentioned in Chapter 3.

Each fiber layer with an orientation angle θ to the z axis can been regarded as the one transformed from a fiber layer where the fiber's principle direction coincides with the z axis. Thus, the stiffness matrix of the inclined fiber layer can be obtained from the one of the straight fiber layer through the stiffness transformation. The stress-strain relations for the straight fiber layer is an orthotropic one, which includes 9 constants:

                              σ xx σ yy σ zz τ xy τ yz τ xz                               =                               C 1111 C 1122 C 1133 0 0 0 C 1122 C 2222 C 2233 0 0 0 C 1133 C 2233 C 3333 0 0 0 0 0 0 C 1212 0 0 0 0 0 0 C 2323 0 0 0 0 0 0 C 1313
Then the stiffness matrix of the inclined fiber layer can be obtained using the coordinate transformation in the (y, z) plane and leads to:

                              σ xx σ yy σ zz τ xy τ yz τ xz                               =                               C1111 C1122 C1133 0 C1123 0 C1122 C2222 C2233 0 C2223 0 C1133 C2233 C3333 0 C3323 0 0 0 0 C1212 0 C1213 C1123 C2223 C3323 0 C2323 0 0 0 0 C1213 0 C1313 σ xx,x + τ xy,y = 0 τ xy,x + σ yy,y = 0 τ xz,x + τ yz,y = 0. (4.4)
Upon substitution of Equation (2.3), (4.3) and (4.4) and noting that the transverse displacements u and v are independent of z as well, the displacement-equilibrium equations for each layer are:

C1111 u ,xx + C1212 u ,yy + ( C1122 + C1212 )v ,xy + ( C1123 + C1213 )w ,xy = 0 C1212 v ,xx + C2222 u ,yy + C1213 v ,xx + C2223 v ,yy + ( C1122 + C1212 )u ,xy = 0 C1213 v ,xx + C2223 v ,yy + C1313 w ,xx + C2323 w ,yy + ( C1123 + C1213 )u ,xy = 0. (4.5)
Combined with the boundary conditions given, this set of equations can be solved by the MG method. This formulation was first derived by Pipe and Pagano [START_REF] Pipes | Interlaminar stresses in composite laminates under uniform axial extension[END_REF] and in their work the finite difference method with very coarse meshes is adopted due to the limitation of the computer capacity at that time. Wang and Crossman [START_REF] Wang | Some new results on edge effect in symmetric composite laminates[END_REF] developed a finite element method based on Pipe's formulation to solve this problem. With the development of computational ability, more accurate solutions have been obtained.

Comparison with the FE model

In this section, the laminated structure is considered both on microstructural and macrostructural scales by simulating a two-ply consisting of different fiber orientations with angles [-15 • / + 15 • ] to the z axis and the corresponding homogenized model. The microscopic model is solved by the MG method and the macroscopic one by the Finite Element method. Fig 4 .2 shows the microstructure where a staggered distribution of the fibers is assumed. The bulk dimensions compared to the fiber radius r are [±8r, ±8r, ±32r]. The fibers have a Young's modulus 10 times larger than the matrix denoted as E f /E m = 10. In each layer, the vertical distance between any two adjacent fibers is 2.2r. The two layers are assumed to be bonded perfectly.

To obtain the effective stiffness matrix of the selected structure, a homogeneous strain boundary condition (HSBC) is applied on one layer that has the same dimension but whose fiber directions coincides with the z axis. This homogenization process is conducted using the MG method. In this step, we can also select a repre- sentative volume element (RVE) to be analyzed, but a large scale dimension will be more convincing considering the latter real case. The results obtained in this step are transformed to conform with the selected structure using the stiffness transformation mentioned in Section 4.2.1. Then in the FE routine, we use the homogenized results under the following conditions: at the bottom surface (z = -32r): u = v = w = 0, at the top surface (z = 32r): u = v = 0, w = 0.01 × 64r, and on the other four surfaces the normal displacement is kept zero. These conditions are chosen because they approximate the actual free edge conditions. In the MG routine, we apply the same boundary conditions while using the ideal laminated structure (Fig 4 .2). Obviously, the displacement w distribution in FE results is smooth due to the homogenization while in the MG results the distribution varies according to the fiber location.

The variation of the displacement w near the interface induces a stress concentration in τ xz . From the above comparisons, we find that the global effects of both components match very well between FE and MG calculations. However, the MG simulation allows one to capture all the strain and stress heterogeneities generated by the micro-structure.

Numerical results

Since the random nature exists in the real structure case and its interface layer thickness is between 2r and 3r, the thickness of interface layer in two similar ideal cases with perfect cylindrical fibers and definite interface layer thickness are selected to be 2r and 3r respectively. Fig 4.9 shows the structure and displacement w distribution at the boundary surface: y = -32r for both the real structure case and the ideal structure case (interface layer thickness = 2r). The variation of displacement w along the x direction can also be observed in of the real structure case.

Fig 4 .10 shows the comparison of ∆w between the real case and two similar ideal cases, where ∆w is the average displacement w along the z direction minus its average along the x direction and describes the global displacement variation. These results are very similar.

Parameterized study

In this section, we investigate several parameters influencing the ∆w.

Orientation

The first parameter is the fiber orientation. We model three different cases: [-10 • /+ 10 • ], [-15 • / + 15 • ] and [-20 • / + 20 • ]. In all these cases, the thickness of the interface layer and fiber layer is kept at 2r and 6r respectively and the ratio of the Young's modulus between fiber and matrix is 10. Fig 4 .11 shows that the displacement variation ∆w increases with the angle. This confirms what is reported in the literature and indicates that a more severe mismatch of the properties between two layers yields a larger displacement variation. The mismatch depends on the angle and there is actually a maximum mismatch for an angle between 0 • and 45 • . The value of the angle for the maximum mismatch depends on the mechanical properties of the individual plies and should be obtained from the homogenized FE model by varying the orientation between the two plies. However, this investigation is beyond the scope of this thesis.

Interface layer thickness

The second parameter is the thickness of the interface layer. Four different thicknesses of the interface layer are simulated: 0r where two different fiber orientation layers touch each other at the interface, 1r, 2r and 5r. In all theses cases, the thickness of the fiber layer is kept at 6r, the fiber orientation is set at [-15 • / + 15 • ] and the ratio of the Young's modulus between fiber and matrix is 10. Fig 4 .12 shows that a thicker interface has a larger amplitude of the displacement variation but has a smaller slope. The amplitude of the displacement variation between the 2r case and the 5r case varies little while the slope varies much, which indicates that the 2r case is more likely to induce the delamination than the 5r case from this point of view.

Fiber layer thickness

The third parameter is the thickness of the fiber layer. Two different thicknesses of fiber layer are simulated: 5.5r and 8r. The thickness of the interface layer is kept to be 5r, the fiber orientation is set at [-15 • / + 15 • ] and the ratio of the Young's modulus between fiber and matrix is 10. Fig 4 .13 shows that the laminated structure has the same displacement variation for different thicknesses of the fiber layer.

Ratio

The fourth parameter is the ratio of the Young's modulus between fibers and matrix. In this comparison, the thickness of the interface layer and the fiber layer 14 shows that the tendency of the displacement variation is steeper and larger when increasing the ratio of the Young's modulus, which indicates that a stronger discontinuity between fiber and matrix is more likely to induce delamination.

E f /E m = 5 E f /E m = 8 E f /E m = 10

Conclusions

This chapter concentrates on the laminated structure and the related free edge effect. Both microscopic and macroscopic structures are investigated. A good agreement is found between the two scales. Furthermore, a real laminated structure is first simulated by the MG method. Several parameters that may influence the early material damage, including fiber orientation, interface layer thickness, fiber layer thickness as well as the different ratios of the Young's modulus between fibers and matrix, are investigated as well.

Introduction

In the previous chapters, the MG solver developed in this work has shown its efficiency in dealing with large scale problems. In this chapter, another challenging problem existing in large scale systems -the numerical simulation of crack propagation, is investigated. In the literature, different models have been proposed to deal with material degradation and failure:

In linear elastic fracture mechanics, the crack is driven by the concept of stress intensity factors that estimate how fast the stress increases in the vicinity of the crack tip. In this model, the displacement discontinuity across the crack surfaces and the stress singularity at the crack tip have to be captured by the numerical model. Finite element based techniques have been used either using remeshing or enrichment [MDB99, DMD + 00]. The main difficulty of this model is to handle the complex three-dimensional crack topology, branching, coalescence and nucleation. Furthermore, the evaluation of the stress intensity factors along a complex curved crack front in three-dimension remains a difficult question, as these quantities are intrinsically defined in two-dimension. The boundary element method has also been applied in this context because it reduces the complexity of the geometrical modeling by recasting the problem into a space of lower dimension.

Damage mechanics [START_REF] Lemaitre | Mechanics of solid materials[END_REF] is based on a constitutive modeling of progressive and diffusive material degradation. A loss of material stiffness is considered in order to obtain the softening response. The numerical version of this model suffers from a critical dependence on the finite element mesh. Damage tends to localize numerically into the smallest elements of the mesh. To circumvent this problem, and also to account for the effect of the material structure that prevents damage to localize in an infinitely thin zone, gradient models or non-local models have been proposed [START_REF] Rhj Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF][START_REF] Zdenek | Nonlocal integral formulations of plasticity and damage: survey of progress[END_REF][START_REF] Vincent Pensée | Micromechanical analysis of anisotropic damage in brittle materials[END_REF]. By introducing an internal length, these extensions of damage models lead to a saturation of the strai gradient which forces damage to diffuse over a zone with a thickness that equals the internal length.

Between these two models, one can mention the cohesive zone models [START_REF] Donald S Dugdale | Yielding of steel sheets containing slits[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF]. In this model, material degradation is assumed to occur within the interface. It has been motivated by the failure mechanisms of certain materials like polymers, in which the fracture process zone reduces to a prolongation of the actual crack over which non zero traction is applied on the crack faces. The numerical implementation of this model which is based on a non rigid behavior of the finite element edges [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Zhou | Dynamic crack propagation with cohesive elements: a methodology to address mesh dependency[END_REF], introduces a dependence to the mesh size and mesh orientation. In this version where all the element edges potentially break into complex crack morphology can be obtained. A mesh independent implementation of cohesive zone models has been established by introducing the traction separation law within the XFEM [START_REF] Jjc Remmers | A cohesive segments method for the simulation of crack growth[END_REF]. However, the difficulty of geometrical modeling of the crack remains in this case.

In recent years, a new approach called the phase field method has been developed [FM98, BFM00, BFM08, MWH10, MHW10, BHLV14], which is able to describe the arbitrary crack surface without the complicated mesh generation and can handle crack initiation. Therefore, it provides an opportunity to work with the MG method to obtain the same efficiency in crack simulation.

In this chapter, the MG solver is further developed, combined with a phase field method, to simulate crack propagation in brittle heterogeneous materials. At first, the phase field method is briefly reviewed in Section 5.2. Next, the phase field method used in this thesis is illustrated in detail. Then in Section 5.3, we discuss the explicit governing equations that work well in the framework of MG. An overall algorithm suitable for MG is built. After that, in order to improve the efficiency of the MG solver for solving the phase field equation, some modifications are made in Section 5.4. At last, we present three crack examples for simple structures in Section 5.5.

Crack simulation -a Phase field method

Modern failure mechanics is based on Griffith's monumental work in the 1920s [START_REF] Griffith | The phenomena of rupture and flow in solids. Philosophical transactions of the royal society of london[END_REF][START_REF] Griffith | The theory of rupture[END_REF], where the importance of potential energy in crack growth was recognized for the first time. Based on that work, a variety of computational methods have been developed to understand and predict the brittle crack behavior [START_REF] Lambert | Dynamic fracture mechanics[END_REF][START_REF] Ted | Fracture mechanics: fundamentals and applications[END_REF][START_REF] De | Fracture in quasi-brittle materials: a review of continuum damage-based approaches[END_REF][START_REF] Rabczuk | Computational methods for fracture in brittle and quasi-brittle solids: state-of-the-art review and future perspectives[END_REF]. However, these classical methods faces the difficulties mentioned above. The phase field method is developed in this context.

Phase field method review

The phase field method for brittle crack simulation can date back to the pioneering work of Francfort and Marigo [START_REF] Gilles | Revisiting brittle fracture as an energy minimization problem[END_REF], who proposed a variational formulation through energy minimization for brittle crack. Bourdin et al. [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] regularized the formulation by approximating the crack surface with a smooth function instead of a sharp crack surface topology. This regularized formulation extends the classical Griffith theory and provably converges (in the Γ-convergence sense [START_REF] Braides | Gamma-convergence for Beginners[END_REF]) to the sharp crack solution. What is more, the smooth description of crack topology does not require any prescription of crack geometry, thus allowing nucleation or branching.

Many crack models have been applied to the phase field method, including quasi-static and dynamic crack propagation [START_REF] Christopher J Larsen | Existence of solutions to a regularized model of dynamic fracture[END_REF][START_REF] Bourdin | A time-discrete model for dynamic fracture based on crack regularization[END_REF][START_REF] Schlüter | Phase field approximation of dynamic brittle fracture[END_REF]. The method has also been applied to some specific applications, i.e Verhoosel et al. [START_REF] Clemens | A phase-field model for cohesive fracture[END_REF] combined the cohesive fracture with the phase field method, Nguyen et al. [NYZ + 15, Ngu15] investigated the crack propagation in real strongly heterogeneous structure obtained using X-Ray tomography. growth, which is usually a complex task for the general case. In a regularized framework, the variable d is assumed to be continuous that it has an auxiliary state between fully broken and unbroken denoting d ∈ [0, 1] and a regularized parameter l denoting the crack actual width (see (5.2)

Regularized approximation of a diffusive crack topology

This approximate function is the solution of the following homogeneous differential equation:

d(x) -l 2 ∆d = 0 in Ω. (5.3)
where ∆ is the Laplacian operator. Note that d(0) = 1 and d(±∞) = 0 and Eq. 5.3 is the Eular equation associated with the variational principle, then:

d = Arg {inf lΓ l (d)} d ∈ {d|d(0) = 1, d(±∞) = 0} . (5.4)
And the crack surface functional Γ l (d) is:

Γ l (d) = 1 2l Ω d 2 + l 2 |∇d| 2 dΩ = Ω γ(d, ∇d)dΩ. (5.5)
γ(d, ∇d) is defined as the crack surface density function per unit volume and is rewritten as:

γ(d, ∇d) = 1 2l d 2 + l 2 |∇d| 2 .
(5.6)

If l → 0, the diffusive crack surface functional Γ l approaches the shape crack topology Γ. This derivation outlined above can be straightly extended to multidimensional cases [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF].

Griffith energy principle

Francfort and Marigo [START_REF] Gilles | Revisiting brittle fracture as an energy minimization problem[END_REF] proposed a variational approach to fracture mechanics introducing the following energy functional for the cracked solid:

E(u, Γ) = E u (u, Γ) + E s (Γ) = Ω\Γ W u (ε(u))dΩ + G c χ n-1 (Γ).
(5.7)

where E u (u, Γ) denotes the elastic energy stored in the cracked body and E s (Γ) represents the energy required to create the crack according to the Griffith criterion.

The term E u (u, Γ) is expressed as the integral of strain energy density W u (ε(u)) in the volume and the term E s (Γ) is written as the product of the fracture toughness G c (energy required to create a unit surface crack) and the Hausdorff surface measure χ n-1 . n is the dimension of the solid and let χ n-1 be the crack length (n = 2) or surface (n = 3). In the phase field framework, the term χ n-1 represents the crack surface functional Γ l (d). Substituting Eq. 5.5 in Eq. 5.7, we obtain:

E(u, Γ) = Ω W u (ε(u), d)dΩ + G c Ω γ(d)dΩ = Ω WdΩ.
(5.8)

where the free energy W = W u (ε(u), d) + G c γ(d).

A unilateral contact formulation

One of the problems using Eq. 5.8 in simulations is that the crack can be interpenetrated under compression mode, as observed in the numerical examples in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF]. The reason is that the different fracture behaviors like tension and compression are not treated differently in Eq. 5.8. To avoid such situations, some unilateral contact formulations have been proposed by researchers in literature [LRC09, AMM09, FRC10, MWH10, MHW10, AGDL15].

In this thesis, we chose the formulation proposed by Amor et al.

[AMM09] since it fulfills the requirements in this work. The main idea of this method is to decompose the trace of the strain tensor into a positive and negative part: trε = tr(ε) + + tr(ε) -.

(5.9)

where a ± = (a ± |a|)/2. The positive part tr(ε) + is believed to correspond to the material expansion in the region where the crack opens. The negative part tr(ε) -is expected to have no contributions to the crack opening and thus to prevent the crack interpenetration. Therefore, the strain energy density W u (ε(u, d)) is reformulated in order to distinguish the contributions due to compression, expansion and shear of the stain energy:

W u (ε(u), d) = 1 2 K 0 tr(ε) 2 -+ {g(d) + k} [ 1 2 K 0 tr(ε) 2 + + µε dev : ε dev ].
(5.10) K 0 is the bulk modulus: K 0 = λ + 2µ n and n = 2 for two dimensional case and n = 3 for three dimensional case. ε dev is the deviatoric part of the strain tensor: ε dev = ε -1 n tr(ε)I, where I denotes the n-dimensional identity tensor. The degradation function g(d) describes the degradation of the positive and shear part of the stored energy with evolving damage. g(d) is chosen to have a simple form as:

g(d) = (1 -d) 2 .
(5.11) which exhibits the following properties:

               g(0) = 1 g(1) = 0 g (1) = 0.
(5.12)

The first equation in Eq. 5.12 denotes the unbroken statue of the material and the second equation denotes the fully broken statue. The third equation guarantees that the elastic energy density function W u converges to a finite value when d approaches to 1. k is a small parameter which is chosen as small as possible to maintain well-posed for partly broken systems [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF].

Governing equations and the overall algorithm

Governing functions

The phase field d describes a continual material statue from fully broken to intact. The crack propagation process under the external loading is approached by the evolution of the phase field d. The whole process is divided into a number of small time steps with an increasing load. In each time step, the load is fixed and the deformation induced is regraded as quasi-static. Hence, the problem can be summarized to solve the two sub-equations: one is the stress equilibrium equation and the other is the gradient-type evolution equation for the crack phase field.

Stress equilibrium equations Observe Eq. 5.10 and rewrite it as:

W u (ε(u), d) = Ψ -(ε(u)) + {g(d) + k} Ψ + (ε(u)).
(5.13)

where

Ψ -(ε(u)) = 1 2 K 0 tr(ε) 2 -and Ψ + (ε(u)) = 1 2 K 0 tr(ε) 2 + + µε dev : ε dev .
Since it is a quasi-static process, we have σ = ∂W/∂ε = ∂W u /∂ε and substituted with Eq. 5.13:

σ(u, d) = ∂ ε Ψ -(ε(u)) + {g(d) + k} ∂ ε Ψ + (ε(u)).
(5.14)

To be more clear:

σ(u, d) =          {g(d) + k} [K 0 tr(ε)I + 2µε dev ] if tr(ε) ≥ 0 K 0 tr(ε)I + {g(d) + k} 2µε dev if tr(ε) < 0.
(5.15)

Then it yields the governing equations to determine the displacement field in the solid:

         ∇σ(u, d) = 0 in Ω σ(u, d) = f∂Ω or u = Ū∂Ω on ∂Ω. ( 5

Overall algorithm

In this work, the process of the crack propagation is considered to be static: the crack grows with the increase of the load and if the load stops increasing so does the crack. Hence, the load P is defined as a function of the time step t, that is:

P i = t i ∆P t i = 0, 1, 2, .... (5.24)
∆P is the unit load increment. In each time step, the algorithm is to solve the displacement field and phase field respectively.

(1) At time step t i-1 , u i-1 , d i-1 and H i-1 are known.

(2) At time step t i , calculate tr(ε), determine the stress forms according to Eq. 5.15, and compute the displacement field u i at each place x according to Eq. 5.16, solved with the MG solver.

(3) At time step t i , compute the positive energy part Ψ + at each place x according to Eq. 5.21, then renew energy history H i according to Eq. 5.22.

(4) At time step t i , compute the phase field d i at each place x according to Eq. 5.23, solved with the MG solver.

(5) At time step t i+1 , (•) i+1 → (•) i and come back to Step (2).

MultiGrid implementations for the phase field equation

In this section, some MG implementation details for the phase field governing equations 5.23 are discussed. The treatment of the modified displacement governing equations 5.16 is similar to the method outlined in Chapter 2 (see also [GRB + 16]).

In order to make explanations more clear, we illustrate these implementations in a one-dimensional context, but they can straightforwardly be extended to higher dimensions.

Difference operator and boundary condition

Consider a one-dimensional body occupying the domain Ω: 0 ≤ x ≤ 1. The first step in MG is to discretize the continual variable d(x) in Ω through a series of isolated point d j , where:

d j = d(x j ),
x j = jh, j = 0, 1, ..., n, and nh = 1.

(5.25) which denotes that j = 0 or = n on ∂Ω and d = 1, 2, ..., n -1 in Ω. Then using a second order accurate discretization, the first equation in Eq. 5.23 on interior point j ∈ [1, n -1] can be replaced by:

2(1 -d j )H j - G c l d j -l 2 d j-1 -2d j + d j+1 h 2 = 0, in Ω.
(5.26)

The boundary condition for d is controlled by the third equation in Eq. 5.23, which denotes: ∂d/∂x = 0, on ∂Ω.

(5.27)

There are several choices how to discretize this boundary condition. For example, for the discretization in point j = 0, one can use the second order downstream form (see p107 [START_REF] Venner | Multi-level methods in lubrication[END_REF]) with the points j = 0, j = 1 and j = 2. However, we have to point out that such kind of discretization will introduce local non-smooth errors at the neighboring points on the boundaries and one has to use some extra local relaxation cycles to smooth these errors. A better way to avoid this problem is to imagine a ghost point, e.g for the point j = 0, use the long second order central discretization to approximate the Eq. 5.27:

d 1 -d -1 2h = 0. (5.28) which implies d -1 = d 1 .
Then we can use the same discretization form for the boundary points as for the interior ones:

2(1 -d 0 )H 0 - G c l d 0 -l 2 d -1 -2d 0 + d 1 h 2 = 0 j = 0. (5.29)
where the item d -1 can be replaced by d 1 . This simple modification maintains the efficiency as in the Dirichelet boundary conditions.

Time efficiency improvement

A typical stress-strain response under tensile load in brittle or quasi-brittle materials before complete failure is shown in Hence, in order to reduce the computational time, different load increments should be chosen in these two regions. The larger load increment is applied in the elastic deformation process. In this work, we set a criterion for the variable d, e.g:

∆P =          ∆P 1 if d < 0.5 ∆P 2 if d ≥ 0.5.
(5.30)

where ∆P 1 = 50∆P 2 . This choice significantly reduces the computational time used for the elastic deformation process.

Another improvement is made by providing a better initial guess for the displacement u and phase field d in the next load increment. Assuming u i and d i are known for the current time step t i , then in the next time step t i+1 , we propose:

         u i+1 = u i + ∆P i+1 u i -u i-1 ∆P i d i+1 = d i + ∆P i+1 d i -d i-1 ∆P i .
(5.31)

This algorithm speeds up the convergence rate of the MG solver and hence reduces the computational cost.

Numerical examples

In this section, we present three numerical examples of three-dimensional crack propagations for simple structure cases. In all the following examples, since every geometry is considered to be repeated, we impose a periodic boundary condition (PBC) with a tensile load in z direction: εzz = 0.01 in each case. The other parameters are the same in all cases and outlined here: the Poisson ratio is chosen to be 0.3; the crack width l is chosen as 2 times of the finest mesh size.

Single inhomogeneity

In the first example, a cubic solid is embedded with a spherical inhomogeneity at the center is simulated (as same in Fig. 3.1). The ratio between the sphere radius r and the cube length L is r/L = 0.1. The Poisson's ratio is chosen to be 0.3 for both matrix and inhomogeneity. The mesh on the coarsest level is chosen to be 8 × 8 × 8 and 5 levels are used, thus leading to a mesh of 128 × 128 × 128 on the finest level and giving about 9.8 million DOFs in total. In this simulation, the rigid and soft inhomogeneity cases are modeled respectively.

Rigid case:

The ratio of Young's modulus between the inhomogeneity and the matrix is 10. Fig. 5.4 shows the evolution of phase field around the inhomogeneity, where the diffusive crack geometry can be found. In the static elasticity model as studied in the homogenization part of Chapter 3, stress concentrations appear at the two poles of the rigid sphere under tensile load in z direction. The results in Fig. 5.4 show that the value of the variable d converges to 1 starting from the places where stress concentrations take place. Due to the totally symmetric structure and load, the phase field calculated is symmetric as well. Fig 5.5 shows the three-dimensional crack evolution in the cube. It is clear to see that the crack initiates from the poles of the rigid sphere and propagates uniformly in a horizontal plane due to the symmetry.

Soft case: In this case, the inhomogeneity is soft and the ratio of Young's modulus between the inhomogeneity and the matrix is 0.1. For the soft case, the stress concentrations take place at the middle plane of z direction. As the results shown in Fig. 5.6, the crack initiates like a ring at the middle plane of the soft inhomogeneity. Then it propagates outward symmetrically.

Interaction between inhomogeneities

In this example, the multi-inhomogeneity case is simulated to investigate the interaction between close inclusions. The geometry and boundary conditions are demonstrated in Fig. 5.7. The dimension of the cube is [10r × 10r × 10r]. A mesh of 128 × 128 × 128 is used on the finest level as in the single inhomogeneity case. The spheres are equally distributed. In the rigid case, all three spheres are 10 times harder than the matrix and in the soft case all three spheres are 10 times softer than the matrix. According to the superposition principle, the stress concentration are enhanced: for the rigid case in the vicinity of the origin O while for the soft case it is the vicinity of the middle between O 2 and O 3 . The numerical results in Fig. 5.8 and 5.9 confirm this understanding very well.

Laminated composites

In this example, a laminated composite structure is investigated as shown in Fig. 5.10. This laminated structure consists of two fiber orientation layers. The fibers in different layers are perpendicular to each other. The dimension of the bulk is

Efficiency

The efficiency of the MG method in solving elliptical-type differential equations like Eq. 5.16 and Eq. 5.23 provides a fast alternative choice to simulate the crack propagation process. A summary of the computational cost for the three previous examples is listed in Table 5.1. A residual of 10 -7 is selected as the convergence criterion, which is much smaller than the discretization error in every example. The code is implemented using the C language and the calculations are performed on a workstation with a 200 GB, 2.2 GHz, using a single core.

Conclusions

This chapter demonstrates the brittle crack propagation simulation combining a phase field method with the MG method. The basic notations of the phase field method are illustrated. The phase field method used in this work is reviewed. An overall algorithm is built. Some modifications are made in the MG implementations to improve the efficiency of the solver. To validate the MG solutions, three simple structure cases are modeled. The results show the ability of MG method to solve the crack propagation with the phase field method.

Chapter 6

Conclusions and perspectives

Conclusions

The work in this thesis is motivated by the large scale computation requirements in composite material simulations. The MultiGrid method is thought to be the most efficient numerical method to solve large scale problems, which has been proven in computational fluid dynamics (CFD). However, the research of the MG method in solid mechanics is far behind. Hence, the main effort concerns in this work an efficient MG solver for simulating the composite structures and predicting the composite behavior. The main contributions in this work can be resumed as follows:

1. An efficient MG solver is built for computing the elasticity problem in solids. This work partly follows Dr. H. Boffy's PhD work: " Techniques multigrilles et raffinement pour un modèle 3D efficace de milieux hétérogènes sous sollicitations de contact" [START_REF] Boffy | Techniques multigrilles et raffinement pour un modèle 3D efficace de milieux hétérogènes sous sollicitations de contact[END_REF]. However, the code is developed independently and is suitably for simulating the composite structures.

2. The accuracy of the MG solutions is validated with analytical solutions in both microscopic and macroscopic aspects. At micro-scale, the MG solutions are compared with the classical Eshelby's solutions for the equivalent inclusion problem. At macro-scale, the MG solutions are compared with the ones obtained from the composite homogenization theories. The MG solutions are reliable in both scales. Several different types of boundary conditions are developed in the framework of MG. The MG method has shown the potential in dealing with the Chapter 6. Conclusions and perspectives composite homogenization process.

3. With adopting Alcouffe's technique [START_REF] Re Alcouffe | The multigrid method for the diffusion equation with strongly discontinuous coefficients[END_REF], which was proposed for the diffusion equation, the current MG solver is able to simulate large property ratio problems in solids. It extends the applications of the MG solver for simulating a variety of materials.

4. As one of the applications in this thesis, the MG solver is applied to investigate the free edge effect in laminate composites, which is usually investigated in a homogenized model due to the limitation of computational cost. In this problem, stress concentrations, induced by the mismatch of the material properties of the adjacent layers, appear on the boundary surfaces and thus result in the initial inter-laminar failures. A sharp displacement evolution in this region is interpreted as a signature of this stress concentration. With the MG solver, several parameters influencing this displacement "discontinuity" are investigated in a microscopic model instead of in a homogenized model. A real laminated structure obtained from using X-Ray tomography is analyzed with the MG method for the first time. This application proves the ability of MG method for solving large complicated problems.

5.

Then the MG solver is extended into the fracture mechanics field. The solver is used to simulate the brittle crack propagation combined with a phase field method. An overall algorithm is built and several modifications are made to ensure that the phase field method works well in the framework of MG. A robust combination of these models is built and is applied to simulate the process of crack propagation. Three examples are presented. This application shows the ability of the MG method for crack propagation simulation with the phase field method.

Perspectives

The work carried out in this thesis has shown the potential of the MG method for solving the large scale problems in solid mechanics. For future investigations, some suggestions can be summarized as follows:

1. Considering the advantages of the MG method in large calculations, it is expected that the MG model can simulate the crack in real measured structures us-ing tomography with phase field methods, e.g the real laminated structure used in Chapter 4. The main obstacle remaining is the sensitiveness of the MG solver to the boundary conditions.

2. Chapter 5 has shown that the phase field method can work well in the framework of MG. However, the phase field method implemented in Chapter 5 [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments[END_REF] cannot prevent material interpenetration when compression and shear load are involved. Hence. it is necessary to investigate other phase field methods, e.g [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF].

3. Parallel computing in MG is also an interesting direction since it can take advantage of multiple cores to significantly reduce the computing time. Avec le développement des techniques d'imagerie telles que la tomographie par rayons X au cours des dernières années, il est maintenant possible de prendre en compte la microstructure réelle dans les simulations des matériaux composites. Cependant, la complexité des composites tels que des fibres inclinées et brisées, les vides, exige un grand nombre des données à l'échelle microscopique pour décrire ces détails et amène ainsi des problèmes difficiles en termes de temps de calcul et de mémoire lors de l'utilisation de méthodes de simulation traditionnelles comme la méthode Eléments Finis. Ces problèmes deviennent encore plus sérieux dans la simulation de l'endommagement, comme la propagation des fissures.

  2.6 A 1D example of error decomposition: different wavelength components. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Structure of grids . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Unstructure of grids . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Influence of lexicographic Gauss-Seidel relaxation on the error:

  (a) suitable for PBC, (b) suitable for HSBC. 56 3.5 Deformation of the RVE under PBC (left) and HSBC (right): (1)(2) normal loading, (3)(4) shear loading. . . . . . . . . . . . . . . . . 58 3.6 The grids on different levels, from level 1 to level 8. . . . . . . . . 59 3.7 Comparison of elastic constants between the MG solution and the analytical predictions on different levels. . . . . . . . . . . . . . . 60 3.8 Comparison of C 1111 for different sphere sizes. . . . . . . . . . . . 61 3.9 Comparison of C 1122 for different sphere sizes. . . . . . . . . . . . 62 3.10 Comparison of C 1212 for different sphere sizes. . . . . . . . . . . . 62 3.11 Comparison of effective Young's modulus E * between the MG solution and the analytical predictions for different Young's modulus of the inhomogeneity E 2 . . . . . . . . . . . . . . . . . . . . . . . 63 3.12 Woven structure. . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.13 Deformation of the woven structure under PBC. . . . . . . . . . . 65 3.14 Comparison of the stiffness matrix for different stiffness ratio E 2 /E 1 . 66 4.1 laminated geometry [PP70]. . . . . . . . . . . . . . . . . . . . . 70 4.2 laminated structure composed of two layers with different fiber orientations (top layer +15 • , bottom layer -15 • ). . . . . . . . . . . . 74 4.3 Displacement w distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 4.4 Displacement w at the intersection between the two planes y = -8r, z = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Shear stress τ xz distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Shear stress τ xz comparison at the interface x = 0. . . . . . . . . . 4.7 Measured structure: yellow-fibers, light blue-matrix, dark bluevoids (top layer +15 • , bottom layer -15 • , E f = 10E m = 100E v ). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Structural mesh: yellow-fibers, light blue-matrix, dark blue-voids (top layer +15 • , bottom layer -15 • , E f = 10E m = 100E v ). . . . . 4.9 Structure and displacement w distribution at the boundary surface: y = -32r for the real case and the ideal case (interface layer thickness = 2r) respectively (from left to right: real structure, displacement w distribution for the real case, ideal structure, displacement w distribution for the ideal case). . . . . . . . . . . . . . . . . . . 4.10 Displacement ∆w comparison between the real case and the ideal case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Influence of fiber orientation on ∆w. . . . . . . . . . . . . . . . . 4.12 Influence of interface layer thickness on ∆w. . . . . . . . . . . . . 4.13 Influence of fiber layer thickness on ∆w. . . . . . . . . . . . . . . 4.14 Influence of the ratio of the Young's modulus on ∆w. . . . . . . . 5.1 Crack topology approximation in a 1D case: (a) Sharp crack; (b) Diffusive crack with the length scale l [MWH10]. . . . . . . . . . 5.2 Crack topology approximation in a 2D case: (a) Sharp crack; (b) Diffusive crack with a length scale l [MWH10]. . . . . . . . . . . 5.3 A typical displacement-force curve for a quasi-brittle material under tensile load [Ngu15]. . . . . . . . . . . . . . . . . . . . . . . 5.4 The evolution of the phase field around the rigid inhomogeneity d under a periodic tensile load εzz . . . . . . . . . . . . . . . . . . . 5.5 Crack propagation around a rigid inhomogeneity under periodic tensile load εzz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Crack propagation around a soft inhomogeneity under periodic tensile load εzz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Geometry and boundary conditions for the three inhomogeneities case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Crack propagation around three rigid inhomogeneities under periodic tensile load εzz . . . . . . . . . . . . . . . . . . . . . . . . . . 102 5.9 Crack propagation around three soft inhomogeneities under periodic tensile load εzz . . . . . . . . . . . . . . . . . . . . . . . . . . 103 5.10 Geometry of the laminated case. . . . . . . . . . . . . . . . . . . 104 5.11 Crack propagation at the interface between two fiber layers under periodic tensile load εzz . . . .

Figure 1 . 1 :

 11 Figure 1.1: Materials used in the Boeing 787 Dreamliner [DS15].

Figure 1 . 2 :

 12 Figure 1.2: Meshes generated for a plate with a circular inclusion: (a) finite element mesh, (b) finite difference mesh.

Figure 1 . 3 :

 13 Figure 1.3: Laminate composite.

  Figure 1.4: The hierarchy meshes generated in MultiGrid for a plate with a circular inclusion.

Figure 2 . 1 :

 21 Figure 2.1: Elastic body deforms under external forces.

Figure 2 . 2 :

 22 Figure 2.2: Two-dimensional geometric deformation of an infinitesimal material element.

Figure 2 . 3 :

 23 Figure 2.3: Cartesian Cauchy stress components in three dimensions.

Figure 2

 2 Figure 2.6: A 1D example of error decomposition: different wavelength components.

Figure 2 .

 2 Figure 2.7: Structure of grids

Figure 2 .

 2 Figure 2.8: Unstructure of grids

Fig 2 .

 2 Fig 2.9 a and b show errors before iteration and after 10 Gauss-Seidel iterations for the relaxation of the Lamé equation. It is clear to see that local strongly oscillation errors are smoothed efficiently by Gauss-Seidel relaxation while the longer wavelength oscillation errors remain.

Figure 2 .

 2 Figure 2.10: Restriction operator I H h and its weighting [Bof12]

  Fig 2.9 (a) and (b): short wavelength fluctuations are smoothed quickly and long wavelength ones remain. Through restricting related components to the coarse grid Ω H , relaxing approximated equation on coarse grid and interpolating corrected errors to the fine grid Ω h , the previous long wavelength fluctuations are removed in Step (3) (see Fig 2.9 (c)). However, it will bring back some short wavelength fluctuations from Step (3) (see Fig 2.9 (c)). Therefore, Step (5) is designed to eliminate these new local fluctuations and it becomes smooth again (see Fig 2.9 (d)).This scheme can be used recursively for any number of levels (see an example in Fig 2.12).

Figure 2 .

 2 Figure 2.12: Different types of MG relaxation cycle: a 4-level case.

Figure 2 .

 2 Figure 2.13: Homogeneous Strain Boundary Condition.

  Fig 2.15): first, we donnot have any explicit boundary conditions for points on the boundary (white circles) and what we have is the Lamé equation for all points, the relaxation scheme

Figure 2 Figure 2 .Figure 2 .

 222 Figure 2.14: Periodic Boundary Condition.

  1 )+ 123456789101112131415161( Ā(θ)e -ιθ 2 + Ã(θ)e ιθ 2 )+ 123456789101112131415161( Ā(θ)e -ιθ 3 + Ã(θ)e ιθ 3 )].

Figure 2 .

 2 Figure 2.19: Sub-boundary (blue line) in a local hierarchical grid.

Figure 2 .

 2 Figure 2.20: Coefficient (right figure, black box, at center between A and B) construction on the coarser level k -1 (blue circles: points on the fine level k; black circles: fine level k points coinciding with ones on the former coarser level k -1).

Figure 2 .

 2 Figure 2.21: Steps of the interpolation (black circles: fine level points, coinciding with ones on former coarser level; green circles: fine level points between two coarse level points, center of edges; red circles: fine level points at the center of surfaces; blues circles: fine level points at the center of cube).

Contents 3 . 1

 31 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2 Microscopic -Eshelby's solution . . . . . . . . . . . . . . . . . . 48 3.2.1 Equivalent inclusion method . . . . . . . . . . . . . . . . . 48 3.2.2 Validation with MG solution . . . . . . . . . . . . . . . . . 49 3.3 Analytical homogenization methods for composite materials . . 50 3.3.1 Hashin-Shtrikman bounds (HS) . . . . . . . . . . . . . . . 51 3.3.2 Mori-Tanaka method (MT) . . . . . . . . . . . . . . . . . . 51 3.3.3 Self-consistent method (SC) . . . . . . . . . . . . . . . . . 52 3.4 Numerical homogenization using MG methods . . . . . . . . . . 53 3.4.1 Coupling between microscopic and macroscopic . . . . . . 54 3.4.2 Boundary conditions: HSBC and PBC . . . . . . . . . . . . 55 3.4.3 Mesh size . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.4.4 Volume fraction . . . . . . . . . . . . . . . . . . . . . . . . 61 3.4.5 Ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.5 Homogenization: woven structure . . . . . . . . . . . . . . . . . 64 3.5.1 Model description . . . . . . . . . . . . . . . . . . . . . . 64 3.5.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . 65 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

D : C 2 Figure 3

 23 Figure 3.1: Cubic matrix containing a spherical inhomogeneity.

Figure 3 . 2 :

 32 Figure 3.2: Stress component comparison between the Eshelby solution and the MG solution.

Figure 3

 3 Figure 3.3: Boundary conditions applied on RVE.

2

 2 Figure 3.4: Different types of RVE: (a) suitable for PBC, (b) suitable for HSBC.

Figure 3 . 5 :

 35 Figure 3.5: Deformation of the RVE under PBC (left) and HSBC (right): (1)(2) normal loading, (3)(4) shear loading.

Fig 3. 7 Figure 3

 73 Fig 3.7 shows the comparison of the elastic constants on different levels with the above analytical predictions. The model studied is isotropic giving C 1111 = C 1122 +2×C 1212 in analytical predictions. However, the MG solution does not satisfy this relation perfectly since the discretize structure is not symmetric completely

Figure 3 . 8 :

 38 Figure 3.8: Comparison of C 1111 for different sphere sizes.

  Fig 3.9 and 3.10 show the comparison between the MG solution and the analytical prediction for different size of spheres. The MG solution coincides with the lower bound of the Hashin-Shtrikman or Mori-Tanaka prediction well until the inhomogeneity volume fraction exceeds 50%, where the Mori-Tanaka prediction becomes unre-

Figure 3 . 9 :

 39 Figure 3.9: Comparison of C 1122 for different sphere sizes.

Figure 3 .Figure 3

 33 Figure 3.10: Comparison of C 1212 for different sphere sizes.

Fig 3 .

 3 Fig 3.11 shows the comparison of effective Young's modulus between the MG solution and the analytical predictions for different cases. Red lines are on behalf

  Fig 3.13 shows the deformation of the block under

Figure 3 .

 3 Figure 3.13: Deformation of the woven structure under PBC.

Figure 4

 4 Figure 4.1: laminated geometry [PP70].

Figure 4 . 2 :

 42 Figure 4.2: laminated structure composed of two layers with different fiber orientations (top layer +15 • , bottom layer -15 • ).

Fig 4. 3

 3 Fig 4.3 shows the displacement w distribution on the boundary surfaces, where the left one shows the results of the FE calculation using a homogeneous structure and the right one shows the results of the MG calculation using the ideal laminated structure. A displacement variation along the x direction is found at the plane y = -8r, the plane visible on Fig 4.3, because of the mismatch induced by the different fiber orientations. If we select a line on this plane, such as z = 0 i.e. in the middle section of the sample, Fig 4.4 shows a good agreement of the local displacement w distribution along the x direction between FE and MG calculations.

Figure 4 . 3 :

 43 Figure 4.3: Displacement w distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers).

Figure 4 . 4 :

 44 Figure 4.4: Displacement w at the intersection between the two planes y = -8r, z = 0.

  Fig 4.5 shows the shear stress τ xz distribution on the boundary surfaces. The shear stress τ xz shows concentration where two fibers with different orientations are close. Fig 4.6 shows the comparison of the local shear stress τ xz at two different positions in the interface x = 0, which shows that τ xz increases sharply when approaching to the boundary. A good agreement between FE and MG is found for both positions while MG results shows the detailed local stress distribution.

Figure 4 . 5 :

 45 Figure 4.5: Shear stress τ xz distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers).

Figure 4

 4 Figure 4.6: Shear stress τ xz comparison at the interface x = 0.

Figure 4 . 7 :

 47 Figure 4.7: Measured structure: yellow-fibers, light blue-matrix, dark blue-voids (top layer +15 • , bottom layer -15 • , E f = 10E m = 100E v ).

Figure 4 . 8 :

 48 Figure 4.8: Structural mesh: yellow-fibers, light blue-matrix, dark blue-voids (top layer +15 • , bottom layer -15 • , E f = 10E m = 100E v ).

Figure 4 . 9 :

 49 Figure 4.9: Structure and displacement w distribution at the boundary surface: y = -32r for the real case and the ideal case (interface layer thickness = 2r) respectively (from left to right: real structure, displacement w distribution for the real case, ideal structure, displacement w distribution for the ideal case).

  Figure 4.10: Displacement ∆w comparison between the real case and the ideal case.

  • orientation = ±15 • orientation = ±20 •

Figure 4 .

 4 Figure 4.11: Influence of fiber orientation on ∆w.

  = 0r interface layer thickness = 1r interface layer thickness = 2r interface layer thickness = 5r

Figure 4 .

 4 Figure 4.12: Influence of interface layer thickness on ∆w.

Figure 4 .

 4 Figure 4.14: Influence of the ratio of the Young's modulus on ∆w.

Figure 5 . 1 :

 51 Figure 5.1: Crack topology approximation in a 1D case: (a) Sharp crack; (b) Diffusive crack with the length scale l [MWH10].

Figure 5 . 2 :

 52 Figure 5.2: Crack topology approximation in a 2D case: (a) Sharp crack; (b) Diffusive crack with a length scale l [MWH10].

  Fig 5.1(b) and Fig 5.2(b)). This auxiliary variable d represents the diffusive crack geometry and is defined as the crack phase field [MWH10]. Miehe et al. [MWH10] derived the regularized crack surface functional Γ l in the one dimension case for the first time and then generalized it to the multidimensional case. A brief review of this derivation is listed here. Consider a crack located at x = 0 and approximate it by an exponential function (Fig 5.2(b)): d(x) = e -|x|/l .

  .16) f∂Ω and Ū∂Ω are the surface forces and displacements on boundary surfaces respectively. Solving this modified static equilibrium equation is like what we do in the previous chapters.

Figure 5

 5 Figure 5.3: A typical displacement-force curve for a quasi-brittle material under tensile load [Ngu15].

  Fig 5.3. The material goes through two main processes: (a)-(b) linear elastic deformation; (b)-(d) nonlinear plastic deformation.

Figure 5 .

 5 Figure 5.4: The evolution of the phase field around the rigid inhomogeneity d under a periodic tensile load εzz .

Figure 5 . 5 :

 55 Figure 5.5: Crack propagation around a rigid inhomogeneity under periodic tensile load εzz .

Figure 5 . 6 :Figure 5 . 7 :

 5657 Figure 5.6: Crack propagation around a soft inhomogeneity under periodic tensile load εzz .

Figure 5 . 8 :

 58 Figure 5.8: Crack propagation around three rigid inhomogeneities under periodic tensile load εzz .

Figure 5 . 9 :

 59 Figure 5.9: Crack propagation around three soft inhomogeneities under periodic tensile load εzz .

Figure 5 .

 5 Figure 5.10: Geometry of the laminated case.

Fig 5 .

 5 Fig 5.11 shows that under tensile load in z direction, the crack starts at the intersections of the interface between two fiber layers. With the load increasing, the isolated cracks merge.

Figure 5 .

 5 Figure 5.11: Crack propagation at the interface between two fiber layers under periodic tensile load εzz .

Contents 6 . 1

 61 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.2 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

  FOLIO ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : GU DATE de SOUTENANCE : 29/09/2016 Prénoms : Hanfeng TITRE : Multigrid methods for 3D composite material simulation and crack propagation modelling based on a phase field method NATURE : Doctorat Numéro d'ordre : 2016LYSEI090 Ecole doctorale : Mécanique, Energétique, Génie Civil, Acoustique (MEGA) Spécialité : Génie Mécanique RESUME :

  1.1 Materials used in the Boeing 787 Dreamliner [DS15]. . . . . . . . 1.2 Meshes generated for a plate with a circular inclusion: (a) finite element mesh, (b) finite difference mesh. . . . . . . . . . . . . . . 1.3 Laminate composite. . . . . . . . . . . . . . . . . . . . . . . . . 1.4 The hierarchy meshes generated in MultiGrid for a plate with a circular inclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Elastic body deforms under external forces. . . . . . . . . . . . . 2.2 Two-dimensional geometric deformation of an infinitesimal material element. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Cartesian Cauchy stress components in three dimensions. . . . . . 2.4 Finite difference discretization for ∂ ∂x (Q ∂u ∂x

  . . . . . . . . . 2.11 Two level MG scheme: a V cycle with the Full Approximation Scheme (FAS). . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Different types of MG relaxation cycle: a 4-level case. . . . . . . 2.13 Homogeneous Strain Boundary Condition. . . . . . . . . . . . . . 2.14 Periodic Boundary Condition. . . . . . . . . . . . . . . . . . . . 2.15 Treatments for PBC in MG. . . . . . . . . . . . . . . . . . . . . . 2.16 Free Surface Boundary Condition. . . . . . . . . . . . . . . . . . 2.17 Local grid refinement. . . . . . . . . . . . . . . . . . . . . . . . . Hierarchical grid. . . . . . . . . . . . . . . . . . . . . . . . . . . 40 2.19 Sub-boundary (blue line) in a local hierarchical grid. . . . . . . . 41 2.20 Coefficient (right figure, black box, at center between A and B) construction on the coarser level k -1 (blue circles: points on the fine level k; black circles: fine level k points coinciding with ones on the former coarser level k -1). . . . . . . . . . . . . . . . . . 42 2.21 Steps of the interpolation (black circles: fine level points, coinciding with ones on former coarser level; green circles: fine level points between two coarse level points, center of edges; red circles: fine level points at the center of surfaces; blues circles: fine level points at the center of cube). . . . . . . . . . . . . . . . . . . . . 44 3.1 Cubic matrix containing a spherical inhomogeneity. . . . . . . . . 49 3.2 Stress component comparison between the Eshelby solution and the MG solution. . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3 Boundary conditions applied on RVE. . . . . . . . . . . . . . . . 53 3.4 Different types of RVE:

	x	List of Figures
	2.18

  . . . . . . . . . . . . . . . . . . . . 105 Background and motivation . . . . . . . . . . . . . . . . . . . .

	xvi		List of Tables
	l			crack width	Nomenclature Chapter 1
	L h , L H		matrix on fine and coarse levels	Introduction
	N	γ, τ		shear strain and stress number of degree of freedom
	ε, σ r h , r H		strain and stress residual on fine and coarse levels
	U	ε *		eigenstrain strain energy
	ε ∞ , σ ∞ u, v, w Contents displacements in x, y, z directions remote uniform strain applied and uniform stress induced
	C V 1 , V 2	1.1	stiffness tensor volume percent of matrix and inclusion
	C 1 , C 2 AMG		stiffness tensor of matrix and inclusion Algebraic MultiGrid
	∆w CLPT		global averaged displacement in z direction Classical laminated plate theory
	∆ DoF		Laplacian operator Degree of freedom
	δ i j FAS		Kronecker delta Full approximation scheme
	Γ FDM		crack surface topology Finite difference method
	λ, µ FE(M)		Lamé coefficients Finite element (method)
	S FSBC		Ehselby tensor Free surface boundary condition
	H HS		strain history function Hashin-Shtrikman bounds
	∇ HSBC		difference operator Homogeneous strain boundary condition
	d LFA		phase field Local Fourier analysis
	E 1 , E 2 MG		Young's modulus of matrix and inclusion MultiGrid
	E f , E m , E v Young's modulus of matrix, fibers and voids MT Mori-Tanaka method
	G c PBC		fracture toughness Periodic boundary condition
	h, H PDE		mesh size on fine and coarse levels Partial difference equation
	I h H RVE		interpolation operator Representative volume element
	I H h SBC		restriction operator Symmetric boundary condition
	SC		Self-consistent method
		K		bulk modulus

1 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.3 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

Table 3 .

 3 2: MG results of stiffness matrix under PBC and HSBC respectively.Figure 3.6: The grids on different levels, from level 1 to level 8.

	Level	1	2	3	4	5	6	7	8
	Points	5 3	9 3	17 3	33 3	65 3	129 3	257 3	513 3
	p	0.0406 0.0524 0.0768 0.0928 0.0942 0.0964 0.0981 0.0982
	Table 3.3: The difference between C MG 1111 and C MG 1122 + 2C MG 1212 on different levels.

Table 3 .

 3 4: Comparison of the stiffness matrix between the entire domain and the partial domain. Comparison of the stiffness matrix for different stiffness ratio E 2 /E 1 .

		5		C * 1111 /C matrix 1111 C * 3333 /C matrix 3333
				C * 1122 /C matrix 1122
	C * i jkl /C matrix i jkl	3 4		C * 1133 /C matrix 1133 C * 1212 /C matrix 1212 C * 1313 /C matrix 1313
		2
		1	1	4	8	12	16	20
					E 2 /E 1
	Figure 3.14:

Table 5 .

 5 1: Summary of computational cost for three examples.[9.6r × 9.6r × 4.8r]. The meshes on the coarsest and finest level are set to be 8 × 8 × 4 and 128 × 128 × 64 respectively, where 5 levels are used, thus resulting in about 4.9 million DOFs on all levels. The fibers are 10 times harder than the matrix.

	Examples	DOFs	Memory used	CPU time in each increment	Number of increments	Total time
	Single inhomogeneity	9.8 million	600 MB	9.8 min	58	9.5 h
	Three inhomogeneities	9.8 million	600 MB	10 min	65	10.8 h
	Laminated structure	4.9 million	310 MB	5.1 min	80	6.8 h

  Par conséquent, il est nécessaire d'étudier des méthodes numériques plus efficaces pour ce genre de problèmes à grande échelle. La méthode Multigrille (MG) est une méthode qui peut être efficace parce que son coût de calcul est proportionnel au nombre d'inconnues. Dans cette thèse, un solveur de MG efficace pour ces problèmes est développé. La méthode MG est appliquée pour résoudre le problème d'élasticité statique basé sur l'équation de Lamé et aussi le problème de la propagation de fissures basé sur une méthode de champ de phase. La précision des solutions MG est validée par une solution analytique classique d'Eshelby. Ensuite, le solveur MG est développé pour étudier le processus d'homogénéisation des composites et ses solutions sont comparées avec des solutions existantes de la littérature. Après cela, le programme de calcul MG est appliqué pour simuler l'effet de bord libre dans les matériaux composites stratifiés. Une structure stratifiée réelle donnée par tomographie X est d'abord simulé. Enfin, le solveur MG est encore développé, combinant une méthode de champ de phase, pour simuler la rupture quasi-fragile. La méthode MG présente l'efficacité à la fois en temps de calcul et en mémoire pour résoudre les problèmes cidessus.MOTS-CLÉS : composites, matériaux composites, Multigrille, stratifié, rupture quasi-fragile, méthode de champ de phase

	Laboratoire (s) de recherche : Laboratoire de Mécanique des Contacts et des Structures
	UMR CNRS 5259 INSA de Lyon
	18-20 rue des Sciences
	69621 Villeurbanne Cedex FRANCE
	Directeur de thèse: LUBRECHT Antonius. A
	Président de jury :
	Composition du jury : VENNER Cornelis. H	DUMONTET Helene
	DINI Daniele	CHATEAUMINOIS Antoine
	RETHORE Julien	SAINSOT Philippe
	LUBRECHT Antonius. A	BAIETTO Marie-Christine

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2016LYSEI090/these.pdf © [H. Gu], [2016], INSA Lyon, tous droits réservés

Phase field equation According to the second law of thermodynamics, the dissipation φ in the Clausius-Duhem inequality under the isothermal process is expressed as:

(5.17)

Written in detail:

(5.18)

Substitute the σ = ∂W/∂ε and obtain:

Consider the crack propagation process is irreversible so that ḋ ≥ 0 and use the fact: k ≈ 0, we have:

where

(5.21)

It is clear that the value of the phase field d is determined by the positive energy part Ψ + (ε(u)). Note the fact that the crack emerged already will not disappear even if Ψ + (ε(u)) becomes very small when unloading. Therefore, Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] proposed a strain history functional H to govern the maximum of elastic energy over history:

The governing equations of the phase field at time step t can be concluded as follows:

• n = 0 on ∂Ω.

(5.23)