Failure under loading is one of the most severe problems and can lead to unpredictable losses and casualties. There is hence a motivation to study failure in material and structural viii

General introduction

Context

Composite materials have been used with outstanding success over several decades thanks to their superior mechanical performance, high damage tolerance, excellent corrosion resistance, and stable physical properties in the current industrial domain. Based on recent data, the application of composites in Boeing airplanes grew from 1 percent in 1969 to 50 percent in 2010. However, despite using 2D laminates for a long time, their use in many structural applications has been limited by manufacturing problems and inferior mechanical properties. Consequently, over the past 30 years, considerable attention has been given to developing advanced polymer composites reinforced with 3D fiber architectures, called textile composites. In textile composites, on the one hand, the 3D mechanical properties are improved thanks to the complex interlaced structure. On the other hand, textile reinforcements can be produced in massive quantities at a reasonable cost by using modern, automated manufacturing techniques. Accordingly, they are attracting growing interest from both academia and the industry.

Motivations

At the design stage of textile composites, the mechanical behavior and the effective properties of the material are generally not known precisely. Experimental tests are usually expensive, and it is sometimes impossible to obtain all properties. Hence much effort goes into developing mathematical tools to get the performance of the composites. The analytical models can generally give quick and effective results, but these coarse approximations show a deficiency in studying complex behavior. Therefore, numerical approaches such as Finite Element Methods (FEM) are introduced to simulate realistic models and complex behaviors.

The conformal mesh, the most used technique in FEM, can discretize the actual shape of the model. However, in the case of complex structures, like textile composites, interpenetration zones and non-physical voids are present between yarns. Although algorithms exist to establish full conformal meshes for these structures, their computational cost is hefty. Consequently, the voxel-based mesh is an exciting alternative to the conformal mesh as it does not require an advanced algorithm and can be used to create the model directly from the image.

vii

CONTENTS

Based on the voxel type mesh, an alternative numerical solver called the Fast Fourier Transform (FFT) based method was initially proposed by [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF]. It can offer the results in good agreement with FEM voxel-based models. Compared with FEM, the FFT method does not need to assemble the global stiffness matrix, which facilitates its parallel implementation over high-performance computation systems, and complex mathematical operations can be easily transformed into simple multiplications in the Fourier space. These characteristics make the FFT possess much higher computational efficiency than conventional FEM for voxel-type meshes. Thus, developments have increased interest in this method for multi-scale material modeling, including this Ph.D. thesis.

In general, the principle of this thesis is to push forward the application of the FFT method in heterogeneous materials such as fibrous composites concerning elastic modeling and damage modeling. This thesis will focus on the microscopic level, where the material is composed of a matrix and fibers.

Objectives

Elastic modeling

Since the introduction of the FFT method, many improved schemes have been proposed and compared with FEM. However, to the best of our knowledge, the performance of the different FFT algorithms is often assessed quantitatively only concerning the macroscopic responses. Although in nonlinear analysis, like damage analysis, the local response is essential (interface response in particular), quantitative comparisons between FFT and FEM have rarely been reported. Therefore, at first, a quantitative comparison between different schemes of FFT and FEM is reported.

Secondly, despite all the advantages of the FFT method, it also presents several drawbacks. Numerical artifacts in the form of spurious oscillations are among the critical issues of this method for solving multiphase elastic problems such as numerical homogenization. For decades, several causes of oscillations have been proposed in the literature, but their impacts have rarely been reported. In this thesis, a deep analysis is accomplished by evaluating the impact of each cause quantitatively.

Thirdly, due to the use of a voxel-based discretization, in heterogeneous materials, this thesis shows that the irregular discretization of the interface is the dominant cause of oscillation. Special treatments such as the composite voxel method or the neighbor voxel average method are presented to deal with these oscillations. An improved composite voxel method relying on a signed distance function is proposed to alleviate the implementation difficulty of the composite voxel method for non-parametrized interface representations such as those obtained from images.

CONTENTS

engineering, including in the case of composite materials. Among damage modeling techniques, the phase-field gains increased attention because it allows dropping the requirement of re-meshing or front tracking methods to follow the evolution of discontinuities. Indeed, the evolution of the phase-field variable completely describes the crack and the damaged regions of the material by regularizing the sharp cracks as a diffusive crack band. Meanwhile, the displacement field and crack propagation are computed simultaneously by minimizing the total energy with a variational approach to brittle fracture.

In the phase field, the characteristic length (l c ) is one of the essential parameters widely discussed in the literature. Even if it is a numerical parameter to describe the degree of crack diffusion, many researchers regard it as a material parameter since it can hugely impact the material mechanical behavior. In the literature, the study of l c is principally restricted to homogeneous materials. In this thesis, a study of the impact of l c in the heterogeneous medium is conducted. It is shown that the choice of the characteristic length for heterogeneous materials affects not only the macroscopic mechanical behavior but also the local crack propagation pattern. Consequently, a length-insensitive phase-field model proposed by Wu in 2017 is studied and implemented into an FFT solver. This implementation includes a modification as compared to previous works to model heterogeneous materials correctly.

Outline

The structure of this thesis is as follows:

• Chapter 1: A general presentation of composite materials, analytical modeling techniques, numerical modeling techniques, and meshing techniques. The Fast Fourier Transform method is also introduced in this chapter, and the detailed mathematical description. The final part is the presentation of some drawbacks of FFT and existing improvements.

• Chapter 2: An introduction of different causes that lead to spurious oscillations and the improved schemes proposed in the literature specifically to attenuate those oscillations. A quantitative comparison between different FFT and FEM schemes and an in-depth evaluation of the impact of each cause is reported.

• Chapter 3: A neighbor voxels averaging technique and an improved composite voxel technique are presented, and the comparisons between each other and with the conventional composite voxel method in terms of oscillations reduction.

• Chapter 4: Different damage modeling techniques and various phase-field models are reviewed, along with their detailed mathematical descriptions.

• Chapter 5: The application of Miehe's phase-field model in the FFT solver and a deep analysis of the sensitivity to the choice of l c for fibrous materials on the micro-level are shown as well as some issues of the method.

• Chapter 6: The implementation of Wu's phase-field model into the FFT solver for heterogeneous materials is detailed and some illustrative numerical examples.

Chapter 1

Introduction of materials and modeling techniques

Chapter overview

In this chapter, first, a general introduction to composite materials and textile composites is presented in Section 1.2, including their advantages, drawbacks as well as their hierarchical structures to emphasize the importance of multi-scale modeling. Analytical modeling approaches are presented in Section 1.3. Numerical modeling techniques are shown in Section 1.4, as well as different meshing techniques such as the conformal mesh and the voxel-based mesh. The Fast Fourier Transform (FFT) method is presented in section 1.5 with its advantages, drawbacks, and recent improvements.

Material introduction 1.2.1 Composite materials

Using materials to make tools is a characteristic of human beings compared to other species. From wood, stone to modern materials, studying and improving materials are the main lines of human civilization evolution. However, commonly, each material has specific properties that make it interesting for certain applications only. For example, concrete is cheap and relatively light, but it can break apart easily under tension. By contrast, steel is strong but expensive and heavy. An idea is to pour the concrete around pre-stressed metal bars to form a cheap, relative light, and strong material. This kind of material is called composite material, in which two or more different constituents are combined [START_REF] Milton | The Theory of Composites[END_REF]. Composite materials have been widely used in human civilization for thousands of years, like the bow in ancient China, which is made of wood and leather in which the wood offers rigidity, and the leather offers toughness.

In modern industrial domains (aerospace, automotive, marine, civil infrastructure, chemical processing equipment, and sports), composite materials have been used with outstanding success over several decades due to their superior mechanical performance, high damage tolerance, excellent corrosion resistance, and stable physical properties [2][3][4]. Based on recent data, the application of composites in Boeing airplanes grew from 1 percent in 1969 to 50 percent in 2010 [5]. From the view of reinforcement type, composites can be separated as particle reinforced, fiber-reinforced, and structural composites. According to the fiber architecture, structural composites can be subdivided into laminated composites and textile composites in which the fibrous reinforcements are interlaced in multiple directions [3,4].

However, despite using 2D laminates for a long time, their uses have been limited by manufacturing problems and inferior mechanical properties. Due to the requirement of dedicated labor (manual lay-up of plies) and expensive facilities, the manufacturing of laminates can be costly. Besides, the wide application of prepreg and fabric plies adds to the difficulty of molding complex shapes. As a result, numbers of machined laminate parts need further processes to build up a complex part. These manufacturing issues limit the widespread use of laminates in the structures of airplanes [3].

To overcome the problems with the manufacturing and mechanical properties of laminates, over the past 30 years, considerable attention has been given to the development of advanced polymer composites reinforced with 3D fiber architectures, which are called textile CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES composites. They are composed of textile reinforcements combined with a binding matrix (usually polymeric), and the term "textile" is used to describe an interlaced structure consisting of yarns.

In textile composites, on the one hand, the 3D mechanical properties are improved thanks to the complex interlaced structures. On the other hand, textile reinforcements can be produced in massive quantities at a reasonable cost by using modern, automated manufacturing techniques. That makes textile composites at the center of the cost and performance spectra, offering significant opportunities for new applications of polymer composites. Thus, they are attracting growing interest from both the academic domain and industry [6].

Following the fabrication techniques, the textile composites can be classified into several families, such as woven fabrics (weaving process), braided fabrics (braiding process), noncrimp fabrics (stitching or knitting processes), and nonwovens [3, 7-9]. Among all the above techniques, weaving is most used thanks to its high production speed and great flexibility to produce a diverse range of 3D composites [7, 8].

Hierarchy of textile composite material

Fiber is the most basic material, which is usually defined as textile raw material. In the manufacturing process, fibers are assembled into yarns by twisting or non-twisting techniques and then fibrous plies and textile composite. Consequently, an important feature of textiles is their hierarchical nature, which can be distinguished into three levels and associated scales:

• Fibers at the microscopic scale.

• Yarns at the mesoscopic scale.

• Fabrics at the macroscopic scale.

Each scale is classified by a characteristic length which can be 0.01 mm for fiber diameter, 0.5-10 mm for yarn diameter and repeating unit cells, and 1-10 m and above for textile structures. Each level is also characterized by dimensionality, where fibers and yarns are mostly one-dimensional while fabrics are two-or three-dimensional. Comprehending this hierarchy is helpful to understand the structure and the mechanical behavior at each scale and hence eases the material design as well as its optimization. [7] At the design stage of textile composites, figuring out the mechanical properties is an essential step. Although WWFE can provide them authentically [10], it is often on a large scale (e.g., macroscopic scale). When meeting issues that need to be studied at a fine scale (e.g., local mechanical or damage fields), the experimental tests are usually expensive and sometimes even impossible. Hence the effort goes into the development of mathematical tools to predict the performance analytically or numerically.

CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

Analytical modeling techniques for composites 1.3.1 Overview

Many simple yet effective analytical methods have been proposed in the literature to compute the properties of heterogeneous materials. This is achieved using a multi-scale approach [START_REF] Abdulle | The heterogeneous multiscale method[END_REF], which, for textile composites, generally consists of two homogenization steps that are presented as follows:

• The first is at the microscale level to obtain the properties of yarns.

• The second is at the mesoscale level from yarns to the final composite properties.

They correspond to the hierarchy of textile composite introduced in Section 1.2.2.

Micromechanics

Homogenization at the microscale can be achieved by micromechanical models [START_REF] Aboudi | Micromechanics of composite materials[END_REF][START_REF] Böhm | A short introduction to continuum micromechanics[END_REF][START_REF] Nemat-Nasser | Micromechanics : overall properties of heterogeneous materials[END_REF][START_REF] Markov | Elementary Micromechanics of Heterogeneous Media[END_REF]. The most used one is the rule of mixtures, which incorporates the Voigt [16] and Reuss [17] models and considers that the materials are arranged in parallel along the longitudinal direction and series in the transversal direction. Next, the Chamis model [START_REF] Chamis | Mechanics of Composite Materials: Past, Present, and Future[END_REF][START_REF] Hopkins | Unique set of micromechanics equations for high-temperature metal matrix composites[END_REF] considers a square reinforcement instead of a cylindrical one hence improving the calculation of the transverse and shear modulus while the longitudinal properties stay the same as in the rule of mixtures. Finally, Eshelby derived analytical solutions for a set of problems involving ellipsoidal elastic inclusions in an infinite elastic matrix [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion, and related problems[END_REF][START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF].

However,it should be kept in mind that these solutions are often referred to as non-interacting approximations and are only valid for small inclusion volume fractions. Generally, they give dependable results only for an inclusion volume fraction V i f < 0.1. In this case, each inclusion can be assumed to be loaded by the same far-field stress σ 0 , so that the contributions of each inclusion into additional strain can be treated separately [START_REF] Kelly | An introduction to composite materials[END_REF]. Models for reinforcement volume fraction exceeding a few percent must account for the interactions between local constituents. This is the case of the Mori-Tanaka scheme [23], which, following the work of Eshelby, approximates the far-field stress acting on an inclusion by an appropriate matrix overall stress. Besides, the self-consistent model is also under the category of interacting models, where the stress acting on the inclusion stays the far-field stress, but the matrix is replaced by an effective medium.

Moreover, many other micromechanical models have been proposed during the past century, such as the Hashin and Rosen models [24], which have been improved in [START_REF] Loos | Chapter 2 -composites[END_REF] for transverse shear properties, the Halpin-Tsai model [26], and the bridging model [START_REF] Huang | Micromechanical prediction of ultimate strength of transversely isotropic fibrous composites[END_REF][START_REF] Huang | Simulation of the mechanical properties of fibrous composites by the bridging micromechanics model[END_REF]. A comparative study can be found in [START_REF] Younes | Comparative Review Study on Elastic Properties Modeling for Unidirectional Composite Materials[END_REF][START_REF] Abaimov | On the closed form expression of the Mori-Tanaka theory prediction for the engineering constants of a unidirectional fiber-reinforced ply[END_REF][START_REF] Wang | A review of analytical micromechanics models on composite elastoplastic behaviour[END_REF].

Mesoscopic scale homogenization techniques

At the mesoscopic scale, the analytical modeling research in predicting mechanical properties of textile composites is encouraged and supported by the need for accurate and less CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES complicated models that would consume less computation time as compared to numerical methods [START_REF] Hallal | Review and comparative study of analytical modeling for the elastic properties of textile composites[END_REF][START_REF] Crookston | A summary review of mechanical properties prediction methods for textile reinforced polymer composites[END_REF][START_REF] Ansar | Modeling strategies of 3d woven composites: A review[END_REF]. The availability of dependable and flexible 3D models is necessary for research and development in advanced industries. Available analytical models for the mesoscopic scale can be divided into the following categories: Classical Laminate Theory (CLT), iso-strain/iso-stress models, mixed iso-strain-stress models, inclusions methods, and the Method Of Cell (MOC).

For 2D woven fabrics, three basic models have been proposed successively based on the pioneering work of Ishikawa and Chou [START_REF] Ishikawa | Elastic Behavior of Woven Hybrid Composites[END_REF][START_REF] Ishikawa | Stiffness and strength behaviour of woven fabric composites[END_REF][37][38] by using the CLT: the mosaic model, the undulation model, and the bridging model. The idea of iso-strain/iso-stress models was initially proposed in [START_REF] Khokar | 3d-weaving: Theory and practice[END_REF][START_REF] Kregers | Determination of the deformability of threedimensionally reinforced composites by the stiffness averaging method[END_REF] to discretize a Representative Volume Element (RVE) into asymmetrical ideal cross-ply laminate blocks and then assemble them under parallel/series conditions. This idea works well to predict the properties of 2D woven composites [START_REF] Scida | A micromechanics model for 3D elasticity and failure of woven-fibre composite materials[END_REF]. To overcome some drawbacks of single iso-strain or iso-stress assemblage, researchers have proposed mixed iso-strain-stress models, such as the parallel-series (PS) scheme [START_REF] Naik | An Analytical Method for Thermoelastic Analysis of 3D Orthogonal Interlock Woven Composites[END_REF][START_REF] Naik | Analytical study of strength and failure behaviour of plain weave fabric composites made of twisted yarns[END_REF] and the 3SHM scheme [44,[START_REF] Hallal | Improved analytical model to predict the effective elastic properties of 2.5d interlock woven fabrics composite[END_REF]. The PS scheme proposes to first assemble the sub-volumes across the loading direction under iso-strain condition (parallel) and then along the loading condition under iso-stress condition (series). As for the 3SHM scheme, the technique is clearly explained in Fig. 1.2.

The most used inclusion model is based on the Mori-Tanaka (M-T) model, where the continuous yarns are divided into sub-volumes that are considered inclusions with circular or elliptical cross-sections. Then the sub-volumes are assembled based on their local orientations (these orientations can be approximated by statistics) [START_REF] Gommers | Modelling the elastic properties of knitted-fabric-reinforced composites[END_REF]47]. It is shown that the M-T model yields better results than the classical iso-strain model, especially with knitted fabrics [START_REF] Gommers | Modelling the elastic properties of knitted-fabric-reinforced composites[END_REF].

The last one is the MOC, of which the fundamental idea is the same as in the methods listed above. The RVE is decomposed into a fiber-matrix system, where the properties are evaluated by the Chamis model and then upscaled with a complementary energy minimization technique [START_REF] Prodromou | The method of cells and the mechanical properties of textile composites[END_REF]. Details are shown in Fig. 1.3. As a summary, a general comparison of mesoscopic scale analytical modeling approaches can be found in [START_REF] Hallal | Review and comparative study of analytical modeling for the elastic properties of textile composites[END_REF].

Conclusions

A graphical overview is shown in Fig. 1. 4, where models are classified based on their predictive accuracy and computational efficiency. All these analytical methods are easy to implement and have a low computational cost.

Nevertheless, most analytical models are still based on the laminate theory and local orientation average approaches. Usually, they have good predictive capabilities in the fiberdominated directions (e.g., longitudinal direction) but not in the matrix-dominated directions, which is problematic for transverse, shear, and out-of-plane properties. Moreover, the real geometry of the complex structure of composites with their manufacturing defects, like voids generation during a manufacturing process [START_REF] Campbell | Chapter 6 -curing: It's a matter of time (t), temperature (t) and pressure (p)[END_REF], is neglected. Additionally, material damage is difficult to analyze using an analytical model. Figure 1.3: Scheme of decomposition and re-composition (multi-scale method) of the MOC [START_REF] Prodromou | The method of cells and the mechanical properties of textile composites[END_REF]. Figure 1.4: The capabilities of the various mechanical modeling approaches for composites [START_REF] Prodromou | The method of cells and the mechanical properties of textile composites[END_REF].

CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

Numerical modeling techniques for composites 1.4.1 Overview

As stated above, analytical approaches are practical but use restrictive assumptions. Thus, numerical methods such as the Finite Element Method (FEM) become increasingly attractive to understand the behavior of heterogeneous materials like textile composites [START_REF] Engquist | Heterogeneous multiscale methods: A review[END_REF]. The great advantage is that numerical modeling can be used to analyze a structure close to reality. Applying contemporary image-based characterization techniques (e.g., micro-CT) and the growing computational capability enable researchers to analyze exact structures and diverse mechanical behaviors. Nonetheless, a complete representation of the microstructure of a composite part in a single model is far beyond the current computational capability, as such a model would have to use micrometer-sized elements. Thus, multi-scale numerical (micro and meso-level) modeling strategies (as shown in Fig. 1.5) have gained popularity to bridge the gap between different material scales [51].

These strategies require a detailed description of the structures, as presented in the following.

Microscale geometry modeling

The microscale geometry of textile composites corresponds to a distribution of fibers and matrix in a domain. As mentioned above, modeling the whole domain is far beyond the current computational capability, so that an acceptable window that can stand for the characteristic size of the microstructure should be chosen for study. This small representative microstructure is called Representative Volume Element (RVE).

Accordingly, the key part of microscale modeling is the fiber (and voids if the voids are considered [52]) arrangement. Naturally, the first type of arrangements are the structured distributions shown in Fig. 1.6. The simplest among these microstructures are the periodic CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES hexagonal (PH0) and periodic square (PS0) arrangements. Their use to evaluate the effective properties of composites can be traced back to the 1960s [53]. Models with hexagonal symmetry (CH3, RH2, CH1, and PH0) have been shown to predict transversely isotropic elastic overall behavior more accurately than models with tetragonal (CS8, CS7, PS0) or mono-clinic (MS5) overall symmetry. The same phenomenon was observed in [START_REF] Ernst | Multiscale progressive failure analysis of textile composites[END_REF], where the authors concluded that the stress state in the square unit-cell is too homogeneous, so the contributions of the fibers to the overall response are underestimated. In contrast, the test results from the WWFE show that the overall predictive capability of the square unit cell is better than the hexagonal arrangement [START_REF] Ernst | Multiscale progressive failure analysis of textile composites[END_REF]. Although structured fiber arrangements can give reasonable estimates of the properties under longitudinal direction loading, the overall behavior under transverse mechanical loading depends markedly on the arrangement of the phase [START_REF] Böhm | A short introduction to basic aspects of continuum micromechanics[END_REF]. Furthermore, in many cases, simple periodic micro-geometries are insufficient to provide an adequate representation of composites as fibers are randomly distributed in reality . Besides, the distributions of microstresses and micro-strains in the composite also strongly depend on the fiber arrangement.

Based on all the above reasons, stochastic fiber arrangements are often introduced into the composite microstructure generation. It was confirmed in [START_REF] Ernst | Multiscale progressive failure analysis of textile composites[END_REF] that the properties estimated are improved with these arrangements, especially the shear stiffness. Recently, a study of the influence of the RVE size on mechanical properties has shown that the influence of the random distribution on global properties is negligible when the RVE size/fiber diameter ratio L/d > 22 [START_REF] Stefanou | Stochastic finite element analysis of composite structures based on mesoscale random fields of material properties[END_REF][START_REF] Savvas | Determination of rve size for random composites with local volume fraction variation[END_REF]. The remaining question is how to introduce this randomness.

Nowadays, increased algorithms are developed to introduce a random packing arrangement. They can be divided into two classes: statical construction and dynamic types. For constructive methods, the assemblies of fibers are accomplished without dynamic motion simula-CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES tion or fiber diameter resizing, e.g. Random Sequential Adsorption (RSA) [START_REF] Torquato | Random sequential addition of hard spheres in high euclidean dimensions[END_REF], Triangulation based method [START_REF] Cui | Analysis of a triangulation based approach for specimen generation for discrete element simulations[END_REF], Stienen model [START_REF] Hernandez-Aquino | Stochastic geometric modeling and analysis of non-uniform two-tier networks: A stienen's model-based approach[END_REF], Sedimentation method [START_REF] Bagi | A quasi-static numerical model for micro-level analysis of granular assemblies[END_REF], Layer-wise method [START_REF] Zsaki | An efficient method for packing polygonal domains with disks for 2d discrete element simulation[END_REF], and Nearest Neighbor Algorithm (NNA) [START_REF] Vaughan | A combined experimental-numerical approach for generating statistically equivalent fibre distributions for high strength laminated composite materials[END_REF]. As for dynamic approaches, the most representative instance is the Molecular Dynamic (MD) [START_REF] Lubachevsky | Geometric properties of random disk packings[END_REF] model.

Each method has its advantages and drawbacks. For example, the RSA approach is the first suggestion because of its robustness. In contrast, with its maximum packing density around 65% [START_REF] Chen | The transverse permeability of disordered fiber arrays: a statistical correlation in terms of the mean nearest interfiber spacing[END_REF], it cannot be used for high fiber volume fraction composites, like textile composites, where the tow fiber volume fraction can reach 75% [START_REF] Gommer | Quantification of micro-scale variability in fibre bundles[END_REF]. Except for RSA, NNA, and dynamic approaches, an issue for all methods mentioned above is that they cannot avoid the fiber contact problem. Recently, Liu [START_REF] Liu | Multi-scale damage modelling of 3D textile reinforced composites including microstructural variability generation and meso-scale reconstruction[END_REF] and Parvathaneni [START_REF] Parvathaneni | Characterization and multiscale modeling of textile reinforced composite materials considering manufacturing defects[END_REF] have pushed and accomplished magnificent work in their thesis. Liu proposed a 2D random microstructure generator based on the MD approach to control the local packing fraction, while Parvathaneni modified the NNA algorithm to have better computational efficiency and generate models that are statistically equivalent to the real microstructure.

Mesoscale geometry modeling

Mesoscopic scale modeling is an important bridge that connects the microscopic and macroscopic levels. Lomov et al. [START_REF] Lomov | Meso-FE modelling of textile composites: Road map, data flow and algorithms[END_REF] introduced a complete roadmap for numerical modeling of textile composites at the mesoscale: from the yarn geometrical modeling and meshing techniques to the boundary conditions and damage modeling. Because of the structural complexity of yarns, mesoscale modeling can be classified into two categories: idealized model and realistic model.

Idealized model

The "idealized model" means a simplification from the real shape. This kind of idealization can be used to estimate macro-mechanical properties of 3D orthogonal woven composites [START_REF] Yu | Numerical analysis of macro-scale mechanical behaviors of 3d orthogonal woven composites using a voxel-based finite element model[END_REF] but also to study the damage growth in woven composites [START_REF] Doitrand | Mesoscale analysis of damage growth in woven composites[END_REF]. Based on state of the art [START_REF] Ansar | Modeling strategies of 3d woven composites: A review[END_REF][START_REF] Lomov | Meso-FE modelling of textile composites: Road map, data flow and algorithms[END_REF][START_REF] Isart | Internal geometric modelling of 3d woven composites: A comparison between different approaches[END_REF][START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF], these geometrical simplifications can be summarized as follows:

• Simplified yarn cross-sections, such as ellipsoidal, lenticular, racetrack, etc.

• Constant yarn cross-section.

• Functional binder yarn propagation path.

• Alignment of warp and weft yarns. This kind of simplified geometries of textile composites can be easily constructed with available software: TexGen [START_REF] Lin | Modelling and simulating textile structures using texgen[END_REF] and WiseTex [START_REF]WiseTex[END_REF].

Textile composites can also be modeled by the Digital Element Method (DEM), which has first been proposed in [START_REF] Wang | Digital-element simulation of textile processes[END_REF] and [START_REF] Said | Kinematic modelling of 3d woven fabric deformation for structural scale features[END_REF]. In the DEM, each yarn is modeled as a bundle of 1D element chains, and yarns are connected with the use of frictionless pins (as shown in Fig. 1.7). Thus, some complex behavior, such as yarns sliding, can be represented. Moreover, the yarn undulations and varying cross-sections can be well-reproduced [START_REF] Green | Mechanical modelling of 3d woven composites considering realistic unit cell geometry[END_REF]. Besides, an analytical approach, which assumes that the yarn geometry can be generated by interpolating several cross-sections along two orthonormal planes, has been proposed in [START_REF] Isart | Internal geometric modelling of 3d woven composites: A comparison between different approaches[END_REF].

Realistic model

A "realistic" model should reproduce all characteristics of the textile geometry as close as possible to reality. With the growth of image treatment techniques, these models are usually constructed from an SEM image. The first attempt has been presented in [START_REF] Faes | Accuracy of 2d fe models for prediction of crack initiation in nested textile composites with inhomogeneous intrayarn fiber volume fractions[END_REF]. The main part of this work has been the image processing on a 2D microscope image, as presented in Fig. 1.8. The yarns perpendicular to the image have been figured out by block detection, and the undulated yarns have been modeled by segments that were designed to fit the neighboring above yarns cross-section closely.

A 3D model reconstruction method has been proposed by Naouar et al. in [START_REF] Naouar | 3d composite reinforcement meso f.e. analyses based on x-ray computed tomography[END_REF]. The authors applied segmentation slice by slice to identify the matrix, warp, and weft yarns and then assign different properties to these identified constituents.

Meshing technique

During several decades of application of FEM, the conformal (body-fitted) mesh has been the most desirable approach because it can capture the real geometry of the microstructure. For modeling textile composites, due to their complex shape, tetrahedra are typically used [START_REF] Ullah | A unified framework for the multi-scale computational homogenisation of 3d-textile composites[END_REF][START_REF] Doitrand | On the influence of fabric layer shifts on the strain distributions in a multi-layer woven composite[END_REF][START_REF] Espadas-Escalante | Mesoscale analysis of the transverse cracking kinetics in woven composite laminates using a phase-field fracture theory[END_REF] despite their lower performance as compared to hexahedra [START_REF] Wu | Meso-level finite element modeling method for mechanical response of braided composite tube with gradient structure in axial direction[END_REF]. Although several examples use hexahedral or hexahedral dominant meshes [START_REF] Matveev | Efficient meshing technique for textile composites unit cells of arbitrary complexity[END_REF][START_REF] Jacques | Application of periodic boundary conditions on multiple part finite element meshes for the meso-scale homogenization of textile fabric composites[END_REF][START_REF] Young | THE HYBRID APPROACH OF A 3D TEXTILE COMPOSITE FINITE EL-EMENT MODELLING TECHNIQUE AT MESO-SCALE LEVEL[END_REF], the meshing procedure is complex. In the conformal mesh, nodes of matrix and yarn elements at the interface between the two different material phases are shared. However, automatic meshing tools for modeling realistic and complex 3D woven architectures often lead to non-physical interpenetrations and voids at the contact zones between yarns (shown in Fig. 1.9). Consequently, manual control algorithms for conformal mesh generation have been developed, such as in [START_REF] Grail | Consistent Finite Element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements[END_REF]. These algorithms work in three steps: Figure 1.9: Illustration of non-physical voids and interpenetrations at the contact zone between two yarns [START_REF] Grail | Consistent Finite Element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements[END_REF].

• Step 1 (Fig. 1.10(a)): Draw the real shape of yarns with some techniques, such as describing the form of the yarn by cross-section and path equations [START_REF] Lomov | Meso-FE modelling of textile composites: Road map, data flow and algorithms[END_REF][START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF]. Hence for each cross-section curve and path curve, there will be an intersection point.

• Step 2 (Fig. 1.10(b)): Calculate the distance between intersection points among contacting yarns. If the distance is less than a specific value, this intersection point is defined as a contacting point.

• Step 3 (Fig. 1.10(c)): Use these contacting points as nodes and generate the mesh.

The above instance [START_REF] Grail | Consistent Finite Element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements[END_REF] has been adopted by Doitrand et al. [START_REF] Doitrand | On the influence of fabric layer shifts on the strain distributions in a multi-layer woven composite[END_REF], who modeled a multi-layer woven composite and studied the influence of layer shifts. Other conformal mesh generation techniques also exist, such as in [START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF][START_REF] Wintiba | An automated procedure for the generation and conformal discretization of 3D woven composites RVEs[END_REF], which are not discussed here. Furthermore, another technique called the Domain Superposition Technique (DST) has been proposed in [START_REF] Jiang | Development of Domain Superposition Technique for the modelling of woven fabric composites[END_REF] and applied for predicting the stiffness matrix of a weave-reinforced composite [START_REF] Gao | A multiscale method for predicting abd stiffness matrix of single-ply weave-reinforced composite[END_REF]. In the DST method, the matrix domain is the global domain leaving the yarns on independent meshes. A coupling technique between the domains ensures the continuity of displacements while the properties of the yarns are corrected due to the superposition.

In summary, the manual mesh control algorithms are usually complex and computationally demanding. Meanwhile, the DST method is neither straightforward to apply nor cannot be representative of interactions between two materials.

Consequently, a voxel-based (pixel in 3D) mesh may be preferred to overcome the difficulties of conformal mesh generation. According to [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF], the generation of a voxel mesh for textile composites can be divided into three steps:

• Step 1: Create a big voxel mesh box filled with the matrix material.

• Step 2: Establish the real geometry of yarns.

• Step 3: Assign all voxels inside the yarns geometry as yarns voxels determined by calculating the position of the voxel center point.

Fig. 1.11 illustrates voxels assignment on varied materials, with the matrix in gray and the yarn in red. Voxel-based meshes have two obvious advantages. The first is that there is no yarn interpenetration and thus no need for manual control. The second one is the straightforward establishment of a realistic model from images obtained using micro-CT or SEM, as described in Section 1.4.3. Hence, the voxel-based mesh is attracting increased attention, especially for complex structure modeling, such as the study of elastoviscoplasticity in polycrystalline [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF], the study of realistic waviness in textile composites [START_REF] Green | Mechanical modelling of 3d woven composites considering realistic unit cell geometry[END_REF], progressive damage in triaxially braided composite [START_REF] Zhang | Progressive damage simulation of triaxially braided composite using a 3d meso-scale finite element model[END_REF], etc.

Using a voxel-based mesh does not come without drawbacks:

• Since the element boundaries do not follow the yarn surfaces, artificial contacts between yarns are introduced, as shown in Fig. 1.12(a).

• Isolated matrix elements: the center of these elements are located within the matrix, but most of their volume belongs to one or several yarns, as shown in Fig. 1.12(a).

• Strong numerical oscillations lie along the interface between yarns, which can be due to the non-smooth discretization of the interface, as shown in Fig. 1.12(b).

(a) (b)

Figure 1.12: The drawbacks of voxel-based mesh: (a) isolated matrix elements and artificial contacts (yarns are represented in red and blue while the matrix is in gray), (b) strong numerical oscillations in the results [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF].

The first issue cannot be solved but can be reduced by refining the mesh. The second problem can be solved by reassignment:

• Re-calculation: calculate the volume fraction of yarns and matrix in these elements.

• Re-assignment: if the matrix volume fraction is lower than the volume fraction of at least one of the yarns in these elements, it is then assigned to the yarn with the biggest volume fraction in the element [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF].

CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

The third drawback, which is the most severe, can be reduced by the average smoothing method proposed by Fang et al. [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF]. In this method, for each voxel along the interface, the stress state is averaged with neighbor voxels in one or several layers with a weight function.

The results in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF] have proven that the oscillations can be significantly reduced by using this technique, and the damage patterns were also better predicted.

Conclusions

In the early 20 t h century, due to the lack of modern computers, many simple but effective and robust analytical approaches have been proposed by researchers for both microscopic and mesoscopic scale modeling of composite materials.

With the development of computers, more attention is transferred to numerical methods. One of the most used methods in the industry and research is the FEM. Thanks to the contributions of researchers and engineers over decades, the accuracy and fidelity of FEM have been proven in many domains. Relying on modern computers and imaging techniques, increasingly complex structures of materials can be represented realistically. This gives some challenges to meshing techniques. The conformal mesh is usually the first choice because it can capture the real shape of the model. However, for modeling complex structures, like textile composites, some drawbacks are hard to overcome, such as the interpenetration problem between yarns. Although some techniques exist for establishing full-conformal meshes, it is often at the cost of tedious programming and significant computation time.

To overcome the issues of conformal meshing, the voxel-based mesh can be preferred. It gets rid of the interpenetration problem and makes it easier to create realistic models directly from images. Despite it presenting some drawbacks, the voxel-based mesh is an attractive method for modeling complex structures, like textile composites.

The following section will present an alternative to the FEM with higher computational efficiency for voxel-based meshes.

CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

1.5 Fast Fourier Transform method (FFT)

Overview

As presented in Section 1.4.4, the voxel-based mesh is a simple way of representing a geometry while overcoming some problems such as interpenetrations. In addition, it can be created directly from an image. If image's resolution is too fine, such as an image output from micro-CT that can easily have over millions of voxels, the computational cost can be too large for FEM. Although parallelization techniques have been proposed, such as domain decomposition [START_REF] Butrylo | A survey of parallel solvers for the finite element method in computational electromagnetics[END_REF] and multigrid methods [START_REF] Zhu | An efficient multigrid method for the simulation of high-resolution elastic solids[END_REF], their applications are not straightforward [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF]. Therefore, there is a demand for approaches that can be easier to parallelize. In 1994, the Fast Fourier Transform (FFT) method was initially proposed by Moulinec and Suquet [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] as a voxel-based methodology that does not need stiffness matrix assembling, unlike conventional FEM.

In the FFT method, we can compute the overall properties of composites, as well as the local distribution of stresses and strains. Its analysis is based on the Lippmann-Schwinger's equation, which is solved iteratively by employing the Green operator of a homogeneous reference medium. As mentioned above, FFT algorithms require data sampled in a grid of regular spacing, allowing the direct use of digital images of the microstructure.

The local strain tensor is calculated by a convolution product with a fourth-rank Continuous Green Operator (CGO) and a polarized term. The convolution product is transformed into a simple tensorial product in Fourier space. The calculation of the polarized term and the CGO multiplication are local operations that can be easily parallelized [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF]. This is less straightforward for the computation of the FFT itself, but packages can easily be found with different parallel implementations of this operation (such as FFTW [103]). Consequently, large-scale simulations based on full-resolution images can be performed with the FFT method.

During the last decades, FFT-based methods have been applied to investigate a wide range of physical phenomena in heterogeneous media, such as eigenstrains / thermal strains [START_REF] Liu | Analytical solution for elastic fields caused by eigenstrains in a half-space and numerical implementation based on fft[END_REF][START_REF] Wang | Analytical solutions for elastic fields caused by eigenstrains in two joined and perfectly bonded half-spaces and related problems[END_REF], crystal plasticity [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF][START_REF] Ma | FFT-based homogenization of hypoelastic plasticity at finite strains[END_REF][START_REF] Suquet | Multi-scale modeling of the mechanical behavior of polycrystalline ice under transient creep[END_REF][START_REF] Otsuka | FFT-based modelling of transformation plasticity in polycrystalline materials during diffusive phase transformation[END_REF][START_REF] Paramatmuni | A crystal plasticity FFT based study of deformation twinning, anisotropy and micromechanics in HCP materials: Application to AZ31 alloy[END_REF][START_REF] Shawish | Full-field polycrystal plasticity simulations of neutron-irradiated austenitic stainless steel: A comparison between FE and FFT-based approaches[END_REF] and damage [START_REF] Li | A damage model of mechanical behavior of porous materials: Application to sandstone[END_REF][START_REF] Cao | A novel FFT-based phase field model for damage and cracking behavior of heterogeneous materials[END_REF]. A number of physics in different types of composites have been investigated, such as the thermoelastic properties of Alumina/Al composites [START_REF] Leclerc | A numerical investigation of effective thermoelastic properties of interconnected alumina/Al composites using FFT and FE approaches[END_REF], the Kapitza interface resistance of composite conductors [START_REF] Monchiet | FFT based iterative schemes for composites conductors with nonoverlapping fibers and Kapitza interface resistance[END_REF], the damage analysis of SiC/SiC materials [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF], the rate-dependent behavior of resin-bonded nonwoven structures [START_REF] Staub | Fast FFT based solver for rate-dependent deformations of composites and nonwovens[END_REF], the effect of heterogeneous interphase on unidirectional fiber composites [START_REF] Wang | Effect of heterogeneous interphase on the mechanical properties of unidirectional fiber composites studied by FFT-based method[END_REF], the damage of laminates [START_REF] Li | A non-local fracture model for composite laminates and numerical simulations by using the FFT method[END_REF] and the damage analysis of textile composites [START_REF] Wang | Progressive damage analysis of 3D braided composites using FFT-based method[END_REF][START_REF] Fang | A coupled FE-FFT multiscale method for progressive damage analysis of 3D braided composite beam under bending load[END_REF].

The FFT method proposed by Moulinec and Suquet relies on a fixed-point scheme ("basic scheme"), which consists in iteratively updating the polarized stress and reapplying the convolution with the CGO through FFTs until convergence. In the following, the mathematical description of the FFT method proposed by Moulinec and Suquet is presented.

Mathematical description of basic scheme

As shown in Fig. 1.13, the mechanical behavior of a composite material composed of fibers and a matrix is governed by the individual behavior of each constituent and by its mi- crostructure. Its effective response to a prescribed path of macroscopic strains or stresses can be determined numerically via the solution of a local problem on an RVE. This RVE is assumed to be a periodic unit-cell V composed of divergent phases that are perfectly bonded, with continuous displacements across interfaces.

The displacement, strain, and Cauchy stress fields are denoted u, ε, and σ respectively. Note that the local strain ε (x) and the displacement u (x) fields can be split into the spatial average and fluctuation terms, which are expressed as follows:

ε (x) = ε * (x) + 〈ε〉 or u (x) = u * (x) + 〈ε〉.x, (1.1)
where ε * (x) and u * (x) are the local fluctuations of the strain and displacement fields respectively with u * (x) being periodic, and the operator 〈-〉 is the spatial average. The macroscopic strain is hence a uniform overall strain that is defined as equal to the spatial average of the strain field ε (x) over the domain. Defining the local fourth-order stiffness tensor as C (x) and introducing a reference material with fourth-order stiffness tensor C 0 , Hooke's law for linear elasticity is written as:

σ (x) = C 0 : ε (x) + τ (x) , ∀x ∈ V, (1.2)
where τ (x) = (C (x) -C 0 ) : ε (x) is called the polarized tensor. Usually, the optimum choice of reference material is

λ 0 = max (λ i ) + min (λ i ) 2 , µ 0 = max µ i + min µ i 2 , (1.3)
with λ i and µ i the Lamé coefficients of the constituents. It can be shown [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF] that the equilibrium equation div

(σ (x)) = 0, ∀x ∈ V (1.4)
combined with the constitutive equation in Eq. (1.2) is equivalent to the periodic Lippmann-Schwinger equation

ε (x) = -Γ 0 (x) * τ (x) + 〈ε〉, ∀x ∈ V, (1.5)
where * denotes a convolution product and Γ 0 is the fourth-order tensor field that is called Green operator. Meanwhile, the strain ε (x) should satisfy the compatibility equation.

ε (x) = 1 2 ∇u (x) + (∇u (x)) T , (1.6) CHAPTER 1.

INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

The ∇ means gradient operation, and (-) T is the matrix transpose operator. Solving Eq. (1.5) is not an easy task in real space because the numerical computation of a convolution product is complicated. At the same time, the term Γ 0 is easier to calculate in Fourier space. Thus, Eq. (1.5) can be easily written in Fourier space as follows:

ε (ξ) = -Γ 0 (ξ) : τ (ξ) , ∀ξ ̸ = 0, ε (0) = 〈ε〉, (1.7)
where the convolution is transformed into a double contraction. Here ε, Γ 0 and τ are the Fourier transforms of ε, Γ 0 and τ respectively. The vector of frequency, denoted by ξ, varies in Fourier space.

Eq. (1.2) is nonlinear because of the polarization tensor. It can be solved using an iterative fixed-point scheme:

Initialization:

(a 0 ) ε 0 (x) = 〈ε〉, ∀x ∈ V (b 0 ) σ 0 (x) = C (x) : ε 0 (x) Iteration (i + 1): ε i and σ i are known (a) τ i = σ i -C 0 : ε i (x) (b) τ i = F τ i (c) ε i +1 (ξ) = -Γ 0 (ξ) : τ i (ξ) , ∀ξ ̸ = 0, ε (0) = 〈ε〉 (d) ε i +1 = F -1 ε i +1 (e) σ i +1 (x) = C (x) : ε i +1 (x) (f) Convergence test: ∥〈ε i +1 -ε i 〉∥ ∥〈ε i +1 〉∥ ≤ e (1.8)
Here ∥-∥ is the L 2 norm, -i is the iteration count, and F() and F() -1 are the Fourier transform and its inverse. The iterative algorithm (1.8) is stopped when the convergence condition is satisfied (e = 1.10 -12 in our calculations).

In numerical calculations, the DFT is used. The unit cell is discretized into a regular grid of N 1 •N 2 pixels in 2D cases and N 1 •N 2 •N 3 voxels in 3D problems. For 3D cases, the coordinates of each voxel labeled by i 1 , i 2 , and i 3 are

x (i 1 , i 2 , i 3 ) = i 1 - 1 2 T 1 N 1 , i 2 - 1 2 T 2 N 2 , i 3 - 1 2 T 3 N 3 , with i 1 = 1...N 1 , i 2 = 1...N 2 , i 3 = 1...N 3 ,
and T i is the period of the unit cell in the i t h direction. The definition of the discrete frequencies ξ i depends on whether N i is even or odd. For the odd discretization types as those used throughout this thesis, the discrete frequencies are given by:

ξ i = - N i -1 2 1 T i , ..., - 1 T i , , 1 T i , ..., N i -1 2 1 T i .
The definition for even N i can be found in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

The Green operator of Moulinec-Suquet ( Γ M-S 0 )

Accompanied by the proposition of FFT methods, the expression of the Green operator was also given in [START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF] and detailed in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. Following Eqs. (1.1), (1.2), (1.4) and (1.6), the expression of Γ M-S 0 (ξ) (the notation of the Green operator of Moulinec-Suquet in Fourier space) depends on the choice of reference material. In most studies, even for anisotropic materials, an isotropic reference material is chosen. In this case, Γ M-S 0 (ξ) is

Γ M-S 0,i j kl (ξ) = δ ki ξ l ξ j + δ l i ξ k ξ j + δ k j ξ l ξ i + δ l j ξ k ξ i 4µ 0 ∥ξ∥ 2 - λ 0 + µ 0 µ 0 λ 0 + 2µ 0 ξ i ξ j ξ k ξ l ∥ξ∥ 4 , (1.9)
where µ 0 and λ 0 are the Lamé coefficients of the reference material. The Dirac delta function is denoted by δ. Eq. (1.9) is valid only if ξ ̸ = 0. Otherwise,

Γ M-S 0 (0) = 0.
(1.10) Because Eqs. (1.4) and (1.6) are continuous, Eq. (1.9) is also continuous. In numerical calculations, the domain is discretized by voxels. Thus, the DFT is applied and the local fields u, ε and σ are always calculated at the centroid of each voxel. In the DFT, ξ is discretized in limited series where the high frequencies are ignored. Therefore, the Green operator of Moulinec-Suquet can also be called the "truncated Green operator".

Advantages and drawbacks of FFT

The advantages of the FFT method can be seen in Section 1.5.2. One is that global stiffness assembly is avoided in contrast to the FEM, another is that most operations are local, hence helping the parallelization. The two global operations, namely the Fourier transform, and its inverse can be easily parallelized with an available library, like FFTW [103]. As compared to the DFT, the FFT relies on a divide-and-conquer algorithm to reduce the computational complexity from O N 2 to O N log N , where N is the size of the DFT matrix. With this simplification, on the one side, the complexity of calculation is reduced, and on the other side, the DFT calculations become easy to be parallelized because the small matrix operations can be calculated separately and then be assembled.

According to Fig. 1.14, FFT can give comparable results as FEM and significantly reduce the CPU time. Thus, it can be competitive, especially in simulating complex structures (e.g., textile composites).

Despite all the advantages above, the FFT method presents several drawbacks:

• Fig. 1.15(a): The convergence of the basic scheme sharply decreases when the mechanical contrast increases between constituents, and convergence cannot be achieved for infinite contrasts. Simulating the voids content in the material can be difficult.

• Fig. 1.15(b): The convergence also strongly depends on the choice of the reference material introduced in the method.

• The FFT method can only work with voxel-based mesh, so the mesh cannot be conformal nor be adapted or refined only in some regions. Furthermore, it is not possible to duplicate nodes or introduce zero-volume interface elements to introduce cracks or contact zones. In other words, for a model resolved in FFT solver, the voids and precracks should also be meshed and assigned with certain properties, although this may not be physical.

• Fig. 1.15(c): The voxel-based models show strong oscillations around the interface between phases.

Another characteristic of the FFT method is that the local fields are intrinsically periodic. This facilitates the modeling procedure for periodic microstructures. However, in the case of non-periodic microstructures, some specific modification should be introduced to cut the periodicity.

Improvements of FFT methods

Since 1994, many researchers have contributed to the development of the FFT method and have proposed solutions to overcome its drawbacks. In 1999, Eyre and Milton [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] intro-CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES duced an accelerated scheme which replaces the operation ε i +1 (x) = -Γ 0 (x) * τ (x) i + 〈ε〉 by .11) In contrast to the basic scheme, in this accelerated scheme, ε i +1 satisfies the compatibility condition only at convergence. The convergence rate of this accelerated method is proportional to the square root of the contrast between the phases, and the optimum choice of reference material can be expressed as:

ε i +1 = ε i + 2 (C + C 0 ) -1 : C 0 : -Γ 0 * τ i -ε i + 〈ε〉 . ( 1
λ 0 = -max (λ i ) • min (λ i ), µ 0 = -max µ i • min µ i .
(1.12)

Eyre and Milton improved the convergence rate at high contrast scenarios, nevertheless, it cannot deal with infinite contrast problems that often occur in the presence of voids or pre-cracks.

In 2001, J.C. Michel et al. [START_REF] Michel | A computational method based on augmented lagrangians and fast fourier transforms for composites with high contrast[END_REF] proposed an augmented Lagrangian scheme to overcome the difficulty at infinite contrast. In this scheme, the value of C 0 should be determined experimentally for every studied microstructure and contrast value, which may be a limitation.

In 2010, Zeman et al. [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients[END_REF] and Brisard [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] proposed simultaneously the use of Conjugate-Gradient (CG) based solvers that can improve both the sensitivity to the phase contrast and to the choice of reference material.

Zeman et al. [START_REF] Zeman | Accelerating a FFT-based solver for numerical homogenization of periodic media by conjugate gradients[END_REF] found that the relation between microscopic strain fields and imposed macroscopic loading strain with the FFT method can be written in a matrix form as:

(I + B) : ε = 〈ε〉.
(1.13)

The Eq. 1.13 can yield at the m t h iteration a strain field in the form

ε m = m i =0 (-B) i : 〈ε〉 (1.14)
which is suitable for a CG solver. Brisard and Dormieux [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF] have proposed to use the Hashin-Strikmann energy principle to solve the problem instead of using the basic scheme, but the principle stays the same.

In 2013, Gélébart and Mondon-Cancel [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF] proposed a combination of the Newton-Raphson algorithm and the CG solver to accelerate the convergence for nonlinear behavior. Based on their results, it can be concluded that this combination is more efficient than the basic scheme for strongly nonlinear materials.

Conclusions

In this section, the FFT method is briefly introduced as well as its mathematical description. As shown in Section 1.5.2, the global assembly of the stiffness matrix is not necessary, and the all operations can be easily parallelized thus making the FFT solvers much faster than conventional FEM solver for voxel-based meshes. Although it also presents several drawbacks, improvements have been progressively proposed since 1994. Section 1.5.4 only includes a few of these improved schemes and many more can be found in the literature [START_REF] Monchiet | Numerical homogenization of nonlinear composites with a polarization-based fft iterative scheme[END_REF][START_REF] Kochmann | A simple and flexible model order reduction method for fft-based homogenization problems using a sparse sampling technique[END_REF][START_REF] De Geus | Finite strain fft-based non-linear solvers made simple[END_REF][START_REF] Wicht | On quasi-newton methods in fast fourier transform-based micromechanics[END_REF][START_REF] Wicht | Anderson-accelerated polarization schemes for fast fourier transform-based computational homogenization[END_REF].

Thanks to the superiority in computational efficiency compared to the conventional FEM solvers, the FFT solvers are progressively applied in various domains where voxel-based meshes are preferred. Apart from the applications presented in Section 1.5.1 in the form of published papers, many works have been carried out and published in Ph.D. thesis reports.

Chen [START_REF] Chen | Damage mechanisms in SiC/SiC composite tubes : three-dimensional analysis coupling tomography imaging and numerical simulation[END_REF] has studied the damage mechanisms in SiC/SiC composites tubes to help the development of new-generation nuclear reactors; M.T. Nguyen [START_REF] Nguyen | Contribution aux méthodes de calcul des propriétés élastiques et de transport des milieux hétérogènes par la Transformée de Fourier[END_REF] has studied the conductivity and the elastic properties of composites and polycrystals, and proposed the application of the FFT method to study the Stokes flow problem in a porous microstructure; Chariere [START_REF] Chariere | Développement de nouveaux matériaux polymères composites allégés à base de micro-sphères creuses modifiées[END_REF] has studied the mechanical behavior of composites with hollow spherical inclusions to help to manufacture lighter composites; Gallican [START_REF] Gallican | Homogenization estimates for polymer-based viscoelastic composite materials[END_REF] has studied the viscoelastic behavior of polymer-based composites; Eloh [START_REF] Eloh | FFT-based modelling of X-Ray Diffraction peaks : application to dislocation loops[END_REF] has studied the X-ray diffraction peaks of single crystals under external loading and the presence of linear defects such as dislocation loops, and D.V. Nguyen [START_REF] Nguyen | Comportement mécanique à température ambiante d'un revêtement de chrome déposé sur un substrat en alliage de zirconium[END_REF] has studied the mechanical behavior of chromium coating on a zirconium alloy substrate under the room temperature environment to help the design of nuclear reactors.

Chapter conclusions

In this chapter, the general introduction of material and modeling techniques have been presented. Analytical modeling is the most straightforward approach to predicting material behavior but shows a deficiency in studying complex behaviors. Numerical approaches such as FEM are more advanced and can simulate realistic models and complex behaviors.

In the FEM, the conformal mesh is the most used technique, as it can represent the real shape of the model. However, in the case of complex structures, such as textile composites, an interpenetration problem between different yarns often occurs. Even if solutions have CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES been proposed, their implementation complexity and computational cost are prohibitive. Voxel-based mesh, where the regular cube mesh is used, is attracting more attention since it can overcome some drawbacks of conformal mesh and have the capability to create the model directly from the image.

Based on the voxel-based mesh, a much higher computationally efficient solver relying on the Fast Fourier Transform (FFT) method was proposed in 1994. Unlike the FEM, the assembly of the global stiffness matrix is avoided in the FFT method. Likewise, most operations are calculated locally. These two characteristics can ease the massive parallelization to take the advantage of modern computing systems. Thus, nowadays, the FFT method is widely applied in various domains, and attracts many researchers worldwide.

Despite all the advantages for the FFT method, the drawbacks cannot be avoided. In this chapter, we have discussed the first two drawbacks and their improved algorithms. For the third drawback, the numerical oscillations, which is one of the focuses of this Ph.D. thesis, details will be presented in the next chapter with a deep analysis on different causes of numerical oscillations to identify the most severe one and a quantitative comparison among improved techniques to find an optimal algorithm.

Chapter 2

Insightful and quantitative analysis of the causes of numerical artifacts and improved schemes for FFT CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

Chapter overview

As discussed in Section 1.5.3, the numerical oscillation (artifact) in the presence of materials discontinuities (heterogeneity) is one of the severe problems shown in the FFT solvers, which affects their precision and fidelity. With the oscillated results, it is hard to do the proper analysis and understand the mechanisms. Therefore, analyzing the causes of oscillations and proposing improved methods to reduce them are important research topics on applying of FFT methods in mechanics.

This chapter will principally follow our publication [START_REF] Ma | Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods[END_REF]. Reviews of possible oscillation causes and improved schemes are presented respectively in Sections 2.2 and 2.3. Section 2.4 introduces the mathematical description of improved techniques chosen for analysis. The quantitative comparison of different improved schemes is presented in Section 2.5. Finally, The deep analysis of each cause of oscillations is shown in Section 2.6.

Oscillation causes reported in the literature

The first cause is the Gibbs phenomenon or Gibbs effect, initially proposed by Henry Wilbraham in 1848 and rediscovered by J. Willard Gibbs in 1898 [START_REF] Gibbs | Fourier's series[END_REF] and well known in signal processing. It was reported in the mechanical field in [START_REF] Eloh | Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains[END_REF][START_REF] Ma | Fft-based solver for higher-order and multi-phase-field fracture models applied to strongly anisotropic brittle materials[END_REF][START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF][START_REF] Kaßbohm | Improved approximations of Fourier coefficients for computing periodic structures with arbitrary stiffness distribution[END_REF]. Fourier series approximate periodic functions by summing numbers of normal trigonometric functions with different frequencies and and up to a given truncation frequency. With a higher truncation, more frequencies can be included in the series, and the function's approximation can be closer and closer to the original function. Nevertheless, Gibbs found that the error between the function and its approximation decreases but does not reach zero when the function is discontinuous, even for a nearly infinite number of frequencies in the series. Instead, the error approaches a finite limit of around 9% of the discontinuity amplitude (as shown in Fig. 2.1). In heterogeneous materials like composites, as the constituents generally have different properties, local fields may be discontinuous, and this may lead to the Gibbs effect.

Figure 2.1: Illustration of Fourier series approximation using N frequencies for a square wave [START_REF] Ma | Fft-based solver for higher-order and multi-phase-field fracture models applied to strongly anisotropic brittle materials[END_REF].

CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

The second cause, proposed in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF][START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Eloh | Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF] is that the Discrete Fourier Transform (DFT) cannot satisfy Shannon's theorem for mechanical fields. Although the Fourier transform requires the frequencies (sampling frequencies) to cover the range from -∞ to +∞, the DFT needs a finite frequency interval from a numerical point of view.

Shannon's theorem, widely applied in signal processing, image processing, and wave recording, requires the sampling frequency to be at least twice the cut-off frequency (i.e., the frequency above which the Fourier transforms of the local fields vanish) to recover the full information. That is why, in the sound wave recording, a music CD is recorded with a sampling frequency of 44.1 kHz to match the human ear capability of which the cut-off frequency is around 20kHz.

It is not possible, however, to define a cut-off frequency for mechanical fields [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. That means there is no resolution of the model (generally, the sampling frequency is equivalent to the model resolution), even if it is very fine, that can fully recover the local fields information by using the FFT solver. Therefore, the Green operator proposed by Moulinec and Suquet is always "truncated", leading to the numerical artifacts.

The third cause is the hourglass effect, which is considered as one of the most important origins that can lead to numerical oscillations in FEM, and which has also been observed for the FFT method in [START_REF] Leuschner | Fourier-Accelerated Nodal Solvers (FANS) for homogenization problems[END_REF]. In Fig. 2.2, a 2D finite element with reduced integration is presented with the center lines as dash lines. The intersection of those lines is the Gauss point used for reduced integration. If this element is submitted to a shear strain with no volume change as shown in Fig. 2.2, the lengths of the dash lines do not change nor the angle between them, which means that all components of stress at the Gauss point are zero. This mode of deformation is thus a zero-energy mode because no strain energy is generated by this element distortion, leading to non-physical results. That is called the hourglass effect. Since the FFT method computes all fields at the center of the voxels, which can be considered as hexahedral finite elements with reduced integration formulation, it may be prone to the same hourglass effect as the FEM with reduced integration. Finally, the fourth cause is linked to the nature of the discretization type used in the FFT method, the voxel type mesh. Similar effects were observed in the FEM simulation of textile composites and discussed in [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF], where oscillations were formed mainly due to the zig-zag interfaces between phases (as shown in Fig. 2.3).

In general, possible origins reported in the literature that led to non-physical artifacts in FFT solvers can be summarized as follows: • Gibbs effect,

• Non-satisfaction of Shannon's theorem,

• Hourglass effect,

• Zig-zag interfaces between phases.

Improved techniques proposed in the literature

Over a decade, many researchers have worked to reduce these non-physical artifacts in FFT solvers. In the Fourier transform, the low frequency represents the global change, and the local fluctuations are always represented by high frequencies. Therefore, a simple approach, proposed in [START_REF] Shanthraj | Numerically robust spectral methods for crystal plasticity simulations of heterogeneous materials[END_REF], is to use low-pass filtering to dampen the results. However, the FFT already contains less information as compared to the continuous Fourier transform and lowpass filtering would only lead to more information loss.

The second strategy in [START_REF] Kaßbohm | Improved approximations of Fourier coefficients for computing periodic structures with arbitrary stiffness distribution[END_REF] is to use "attenuation factors" based on the DFT of linear spline functions, which means an interpolation operation is carried out to fill the missing information between the discrete values.

On the other hand, researchers found that the numerical oscillations in the FFT can be due to the creation of the "truncated" Green operator. In the FFT-based algorithm, the applicability of the DFT is based on the assumption that the function is periodic, which is not the case in a space-limited computational domain. Therefore, another idea is to create a quasi-continuous Green operator to match the perfect continuous Green operator as close as possible.

CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

Consistent periodized discrete Green operators, established based on this idea, were proposed in [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF][START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF] for FFT methods based on the Hashin-Shtrikman variational formulation and in [START_REF] Eloh | Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains[END_REF] for those based on the Lippmann-Schwinger equation. In the DFT, the sampling frequency ξ varies from -∞ to +∞, but this is not possible in practice. However, thanks to the periodicity, which is an intrinsic property of the Fourier transform, the frequency ξ can be replaced by ξ = mN + l , where l ∈ [0, N -1] and m varies from -∞ to +∞. Here N is the resolution of the model. For example, for a two-dimensional model, the consistent Green operator can be written as follows:

Γ cons 0 ξ i j = i j sinc πi N sinc π j N ∞ m,n=-∞ (-1) m+n (mN + i ) 2 nN + j 2 Γ 0 ξ (mN+i )(nN+ j ) , (2.1)
where m and n have the same meaning. Compared to the Green operator proposed by Moulinec and Suquet, Eq. (2.1) gets rid of the dependency on the discretization of the real model, which can be refined as needed by increasing m or n.

According to [START_REF] Eloh | Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains[END_REF], the consistent Green operator can reduce the oscillations, and good results can be obtained when m is superior to 81. However, there are two main drawbacks. The first one is that the peak stress computed using this operator is lower than that computed using the M-S (Moulinec and Suquet) operator. It is also lower compared to the analytical solution. Another issue is the computational cost, as the calculation of the consistent Green operator is much more computationally demanding, especially for 3D cases. Therefore, a filtering approach was finally used in [START_REF] Brisard | Combining Galerkin approximation techniques with the principle of Hashin and Shtrikman to derive a new FFT-based numerical method for the homogenization of composites[END_REF].

Borrowing ideas from FEM and finite differences, another strategy is to compute Discrete Green Operators (DGO) based on finite difference discretizations. This technique was firstly reported in [START_REF] Müller | Mathematical vs. Experimental Stress Analysis of Inhomogeneities in Solids[END_REF] and followed by Dreyer [START_REF] Dreyer | Approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials[END_REF] and Brown [START_REF] Brown | Discrete fourier transforms and their application to stress-strain problems in composite mechanics: A convergence study[END_REF] where a centered finite difference scheme was considered. Backward and forward finite difference schemes were published by Willot in [START_REF] Willot | Fast Fourier Transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media[END_REF]. Because of its simplicity and robustness, this research direction has been followed by many authors in the past few years and improved approaches have been proposed [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Willot | Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields[END_REF][START_REF] Schneider | Computational homogenization of elasticity on a staggered grid[END_REF][START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF][START_REF] Vidyasagar | Predicting the effective response of bulk polycrystalline ferroelectric ceramics via improved spectral phase field methods[END_REF] and their efficiency has been proven on diverse applications [START_REF] Berbenni | A numerical spectral approach for solving elasto-static field dislocation and g-disclination mechanics[END_REF][START_REF] Berbenni | A Fast Fourier Transform-Based Approach for Generalized Disclination Mechanics Within a Couple Stress Theory[END_REF][START_REF] Djaka | Field Dislocation Mechanics for heterogeneous elastic materials: A numerical spectral approach[END_REF][START_REF] Lebensohn | Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms[END_REF].

Among the different schemes, the discrete Green operator calculated based on the centered finite difference scheme on a rotated grid, called "rotated scheme" [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], is one of the most promising schemes. This scheme was implemented in the open-source software AMITEX-FFTP [158] developed by Lionel Gelebart et al., and it was successfully applied in [START_REF] Djaka | Field Dislocation Mechanics for heterogeneous elastic materials: A numerical spectral approach[END_REF] for solving field dislocation mechanics problems and in [START_REF] Yan | Meso-scale modelling of 3D woven composite T-joints with weave variations[END_REF] for analyzing the failure mechanisms of a SiC/SiC composite. Another very interesting discretization scheme based on linear hexahedral element was proposed by Schneider et al. [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] and adapted both to the basic scheme and conjugate gradient solvers. This scheme was presented in two versions, one with reduced integration and the other with full integration, similarly to FEM. It was analytically demonstrated that the scheme with reduced integration is equivalent to the "rotated scheme" proposed in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF].

In this section, four improved schemes are introduced. The first two, low-pass filtering and interpolation schemes, may not be sufficient to fill up the lost information. The consistent

CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

Green operator can reproduce the continuous operator but it requires significant computational resources, especially in the 3D cases. The idea of calculating a DGO based on finite difference schemes seems to get a balance between the accuracy and the computational complexity. Therefore, the improved techniques proposed by Willot [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], and Schneider [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] are chosen for deeper analysis as well as the original basic-scheme [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF]. The mathematical description of the Willot [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] and Schneider [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] schemes will be presented in the next section.

Finite-difference based improved techniques

The green operator of Moulinec and Suquet proposed in [START_REF] Moulinec | A fast numerical method for computing the linear and nonlinear mechanical properties of composites[END_REF][START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF][START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] has already been presented in Section 1.5.2. The methods described here are the rotated scheme proposed by Willot and Schneider's scheme based on the linear hexahedral elements.

Willot's rotated scheme

Two dimensions

Inspired by the finite difference technique, Willot modified the equilibrium equation and the strain expression, shown in Eqs. (1.4) and (1.6). They are given as follows:

k * i (ξ) σ i j (ξ) = 0, ε i j = 1 2 k i (ξ) u j (ξ) + k j (ξ) u i (ξ) , (2.2)
where k * and k represent the discrete divergence and gradient operators, respectively. A hat ( -) indicates a variable in Fourier space while no hat means the real space. In the traditionally centered scheme, k is usually taken as

k WI-C i (ξ) = Jsin (ξ i ),
which corresponds to the divergence of stress and the gradient of displacement:

∂ j σ i j (x) ≈ σ i j x + e j -σ i j x -e j 2 , ∂ j u i (x) ≈ u i x + e j -u i x -e j 2 ,
where e i and e j represent the unit directions that can be referred to as e 1 and e 2 in Fig. 2.4, and J = -1. In [START_REF] Willot | Fast Fourier Transform computations and build-up of plastic deformation in 2D, elastic-perfectly plastic, pixelwise disordered porous media[END_REF], where the forward-and-backward finite difference scheme was proposed, the k can be written as:

k WI-W i (ξ) = e Jξ i -1,
with the divergence of stress and the gradient of displacement:

∂ j σ i j (x) ≈ σ i j (x) -σ i j x -e j , ∂ j u i (x) ≈ u i x + e j -u i (x) .
Integrating Eq. (2.2), the new Green operator reads:

Γ WI 0,i j kl = k i (ξ) k m (ξ) C 0,m j kn k * n (ξ) -1 k * l (ξ) sym , (2.3)
where the symbol sym indicates the minor symmetrization with respect to the indices i , j and (k, l ). As stated previously, this definition is only valid for non-zero frequency vectors so that Eq. (1.10) can still be applied. Also note that the Moulinec-Suquet Green operator can be recovered if k i (ξ) is set as k i (ξ) = Jξ i .

CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

However, the centered finite difference scheme cannot define the Green operator Γ WI-C 0 at three frequencies ξ = (π, 0), (0, π) and (π, π). At the same time, the backward-and-forward scheme can break the symmetry in certain models. Therefore, the rotated scheme is proposed to overcome that. In this scheme, the displacement field and the divergence of the stress field are evaluated at the four corners of the pixels while the strain and stress fields are evaluated at the centers of the pixels, where the rotation angle is normally 45 • . This is shown in Fig. 2.4, where the original basis (e 1 , e 2 ) superimposed and the rotated basis ( f 1 , f 2 ) can be expressed as follows: With this new basis, the local problem equations can be written as:

f 1 = e 1 + e 2 2 , f 2 = e 2 -e 1 2 .
σ i j (x) = C i j kl (x) ε kl (x) , (2.4a
)

σ i 1 (x) -σ i 1 x -2 f 1 + σ i 2 x + f 2 -f 1 2 -σ i 2 x - f 2 + f 1 2 , ( 2.4b
)

ε kl (x) = 1 2 2 u k x + f 2 2 -u k x - f 2 2 + u l x + f 1 2 -u l x - f 1 2 , ( 2.4c)
where x is the pixel center and the x ± f i / 2 are the pixel corners. Re-expressing Eq. (2.4) in the original Cartesian basis and applying the Fourier transform, k can be expressed as:

k WI-R i (ξ) = J 4 tan ξ i 2 1 + e Jξ 1 1 + e Jξ 2 .
(2.5)

If N i is even, the frequency can be ξ = (π, π). Thus the k WI-R cannot be defined. An essential treatment should be carried out that k WI-R = 0 and Γ WI-R 0 = 0.

Three dimensions

Following the same logic, in 3D, k WI-R i can be extended as With the assumption of isotropic reference material, the full version of the Green operator of Willot's rotated scheme Γ WI 0 in the 3D case can be given as:

k WI-R i (ξ) = J 4 tan ξ i 2 1 + e Jξ 1
Γ WI 0,i j kl = λ 0 + 2µ 0 r i r * l δ j k sym + λ 0 r i r * l s j k sym -Re r i r * j Re r k r * l -µ 0 r i r j r * k r * l µ 0 2 λ 0 + µ 0 -λ 0 r 2 1 + r 2 2 + r 2 3 2 ,
where

r i = k i /∥ k∥, r * i = k * i /∥ k * ∥, ∥ -∥ is the L 2 norm
and s is the symmetric second-order tensor:

s j j = 4Im r i r * k 2 , s j k = -4Im r k r * j Im r k r * i , i ̸ = j ̸ = k.
Here Re (.) and Im (.) are the real and the imaginary parts of the complex values, respectively.

Schneider's hex scheme

The work of Willot [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] is interesting not only for the new DGO based on the rotated scheme but also for the comparison between various finite difference schemes and how they lead to various discrete Green operators. Based on this work, Schneider et al. [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] proposed to use the finite element method instead of the finite difference method to construct the DGO. They compared two types of finite elements, viz. the classical trilinear hexahedral element with full integration and the one with reduced integration. Following the notations of [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], it is reminded that the coordinates of the integration points for the reference fully integrated hexahedral element [0, 1] 3 are

z β = a β 1 , a β 2 , a β 3 , β ∈ {0, 1} 3 with a 0 = 3 -3 6 , a 1 = 3 + 3 6 , CHAPTER 2.

INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

and the finite element basis shape functions φ α are of the form: Similarly to Eq. ( 2.3) the DGO for hexahedral elements is

φ α = 3 i =1 h α i (z i ) , α ∈ {0, 1} 3 , z ∈ [0, 1] 3 with h 0 (a) = 1 -a, h 1 (a) = a, a ∈ [0, 1] .
Γ S 0 = B K B H K M 0 B K -1 B H K ,
where B K is the discrete symmetric gradient operator and -B H K is the discrete divergence operator at voxel K. Note that B H K is the Hermitian conjugate (the transpose of the complex conjugate) of B K . The operator M 0 is equal to C 0 when using reduced integration. Otherwise, M 0 ∈ R 8×3×3×8×3×3 is defined as C 0 applied independently at each integration point

M 0,αi j βkl = C 0,i j kl , α = β 0, α ̸ = β , α, β ∈ {0, 1} 3 , i , j , k, l = 1 . . . 3.
The formula for the discrete symmetric gradient operator B K is similar to that in Eq. (2.2), whereas this symmetric gradient is now computed at each integration point β of each voxel K and thus involves a gradient operator k K,β . There are always eight symmetric gradient operators BK,β ∈ C 3×3×3 to compute per voxel:

• in full integration, B K = BK,β β∈{0,1} 3 ∈ C 8×3×3×3 and B H K ∈ C 3×3×3×8 , • in reduced integration, B K = 1 8 β∈{0,1} 3 BK,β ∈ C 3×3×3 and B H K ∈ C 3×3×3 .
Consequently, we have

BK,βi j k = 1 2 k K,βi (ξ) δ j k + k K,β j (ξ) δ i k .

CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

The discrete gradient operator k K,β (ξ) ∈ C 3 is given by

k K,β,i (ξ) = α∈{0,1} 3        exp -2πJ 3 j =1 ξ j α j N j (-1) α i +1 k = 1 . . . 3 k ̸ = i h α k z β k       
, while satisfying Eq. (1.10).

Summary

In this section, two improved schemes, Willot's rotated scheme and Schneider's hex scheme, have been presented. They improve the procedure of Green operator calculation to reduce the spurious oscillations. It has been proven in [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] that the hex scheme of Schneider with reduced integration is equivalent to the rotated scheme of Willot. It is reminded that reduced integration is well known to produce numerical artifacts known as the hourglass effect in FEM computations. A similar difference between reduced and full integration when using FFT-based numerical methods has been reported in [START_REF] Leuschner | Fourier-Accelerated Nodal Solvers (FANS) for homogenization problems[END_REF].

Quantitative comparison between FFTs and FEMs

When dealing with complex analysis, like damage analysis, the local response is important (interface response in particular). As illustrated in Fig. 2.7, oscillations are also present even with the rotated scheme and FEM as well. Besides, the oscillations of FEM and FFT are so similar that they cannot be distinguished with only the figure.

Therefore, it is interesting to evaluate each cause of oscillations and analyze its impact, as proposed in this section through:

• Quantitative comparisons of the influence of different Green operators using FEM results as reference for different models proposed in [START_REF] Moulinec | A FFT-Based Numerical Method for Computing the Mechanical Properties of Composites from Images of their Microstructures[END_REF][START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF][START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF]) in order to conclude on an optimal choice.

• Identification and analysis of the causes of spurious oscillations in FFT and evaluation of their impacts quantitatively.

Introduction of test case models

This thesis focuses on numerical methods so that our conclusions should be as general and universal as possible. However, it is challenging to generate a universal model for complex structures like textile composites. To put in evidence all artifacts presented in section 2.2, three simple but representative models are proposed in Fig. 2.8, and we call them models A, B, and C. present in FFT with rotated scheme proposed in [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF] (left half part) and FEM (right half part) under tension in x-direction with the contrast between the fiber and matrix Young's moduli set to

E f /E m = 22.
Model A represents a case with two phases in parallel, thus without strain concentration points. In model B, a square inclusion is inserted into the model, with four strain concentration points. Note that the meshes of models A and B are intrinsically conformal if the fiber/inclusion volume fraction (V f ) is chosen appropriately. As for model C, it is the same as model B but with a circular inclusion. In model C, non-smooth zig-zag interfaces are present and the mesh is not conformal anymore if it is a voxel-based mesh.

Note that these are 3D models, even though only one voxel is used in the z-direction, and for model A, the width (size in y-direction) is also set as one voxel size. The fiber volume fraction is set to V A f = 0.6 for model A, V B f = 0.36 for model B and V C f = 0.55 for model C. For all three models, phase 1 (red part) is assigned as fiber, and phase 2 (blue part) is the matrix. The fibers are composed of E-Glass, while the matrix is Epoxy PMR-15 with properties as shown in Table 2.1.

Material type E(MPa) ν

E-Glass 72 000 0.22

Epoxy PMR-15 3750 0.375 Table 2.1: The elastic properties of fibers and matrix [START_REF] Abaimov | On the closed form expression of the Mori-Tanaka theory prediction for the engineering constants of a unidirectional fiber-reinforced ply[END_REF].

In the following, the FFT method using the CGO of Moulinec-Suquet is denoted as M-S, the method using Willot's rotated DGO (equivalent to Schneider's DGO with reduced integration) denotes W-S and the method using Schneider's fully integrated DGO is denoted as SF. Regarding the FEM, the voxel-based hexagonal FEM with reduced integration is called FEMR while FEMF refers to the fully integrated FEM with the same mesh type as FEMR. Furthermore, it should be pointed out that all FFT results shown below are computed using an inhouse python code while the FEM results are calculated using the Abaqus software.

Mesh convergence test

At first, a mesh convergence test is performed on the three models for both FFT and FEM. The mesh resolutions are denoted as N A , N B and N C for models A, B, and C, respectively. It is defined here as the number of voxels in the x and y directions. The resolutions used for model A are N A = 5, 15, 45, 135, and 405. The same resolutions are used for model B. The fiber volume fraction is exactly obtained for both models for all chosen resolutions. For model C, the target fiber volume fraction V C f = 0.55 can never be perfectly reached due to the intrinsic nature of voxel meshes, even with exceptionally fine resolutions. Therefore, the first step is to for model C verify the fiber volume fraction convergence by varying the mesh resolution. Then, for all models, the convergence in terms of macroscopic properties should be verified. Herein, FEMR is taken as the reference.

As shown in Figs. 2.9(a) and (b), the macroscopic properties for model A do not show any sensitivity to the model resolution. Different FFT methods as well as FEM show good consistency for all resolutions on model A with respect to E x and G x y , where E x means Young's modulus in x-direction and G x y stands for shear modulus in x y plane. As shown in Figs.

2.9(c) and (d), macroscopic properties for model B are converged when N B ≥ 135. The convergence for E x occurs with iteratively improving lower bounds for all methods except for SF which gives a higher bound. For G x y , all methods give higher bounds.

As aforementioned, because of the intrinsic characteristic of model C, the mesh convergence of fiber volume fraction is first verified, as presented in Fig. 2.10. The mesh convergence for fiber volume fraction begins at N C = 75 (which corresponds to the red circle in Fig. 2.10). Hence, the resolutions chosen for the tests of elastic properties convergence are N C =75, 225, 375, representing the small, medium, and large resolutions, respectively.

REP = E FFT -E FEMR E FEMR ,
where E stands for the elastic modulus (E x or G x y ).

As shown in Table 2.2, the properties of model C significantly depend on the fiber volume fraction (V C f ). From N C = 75 to 225, the V C f decreases from 0.5508 to 0.5502, and the E x value drops around 100 MPa, which cannot be neglected. From N C = 225 to 375, this change on the properties becomes more acceptable. Furthermore, the relative errors of all FFT methods are smaller than 1% for N C ≥ 225. Thus, for model C, the choice of resolutions in the following begins at N C = 225. Based on the tests above, we can see that FFT methods can accurately calculate the homogenized elastic properties. The error between FFT methods and FEM is typically less than 0.1%. Considering that the FFT methods have a significant advantage in computational efficiency, they can be a desirable alternative to FEM. This is only about the homogenized properties, however, and the local responses still need to be assessed. Hence, in the next section, the quantitative comparison of local responses will be presented to show the difference between FEM and FFT methods.

Quantitative comparison of local responses

In this part, all three models are subjected to a macroscopic strain along the x-direction (transverse direction) with a value 〈ε〉 xx = 0.001 while stress-free conditions are imposed in the other directions (〈σ〉 y y = 〈σ〉 zz = 〈σ〉 x y = 〈σ〉 xz = 〈σ〉 y z = 0). These mixed-type loading conditions are applied using the method presented in [START_REF] Kabel | Mixed boundary conditions for FFT-based homogenization at finite strains[END_REF]. It is reminded that periodic boundary conditions are intrinsic to FFT methods. Due to stress-free boundaries, Poisson's effect is active during loading.

Based on the mesh convergence test, the resolutions chosen for model A and model B are Here, the FEM results with voxel meshes are chosen as a reference to assess the performance of different FFT methods on local fields. Due to our remarks on the hourglass effect, we consider FEM with reduced integration (FEMR) and full integration (FEMF). The chosen error measure is the relative absolute error of the local strains: where FFTs denotes the result obtained using an FFT solver and FEMs is the result obtained from the FEM with reduced or full integration. ε xx is the strain component chosen for comparison for models A and B. Regarding model C, the interface radial strain (ε r r ) and the interface tangential strain (ε θθ ) are chosen, which can be expressed as:

N A/B =135,
RE = |ε FFTs -ε FEMs | |ε FEMs | , ( 2
ε r r = ε xx • cos (θ) 2 + ε y y • sin (θ) 2 + ε x y • sin (2 • θ) , (2.8a
)

ε θθ = ε y y -ε xx • sin (θ) • cos (θ) + ε x y • cos (θ) 2 -sin (θ) 2 , (2.8b)
where the components can be seen in Fig. 2.11. Note that all strain values along the interface are taken on the matrix side.

First, as shown in Fig. 2.12, the relative error between FFT and FEM results, both FEMR and FEMF, are null for model A. Indeed, FFT and FEM results are identical for this model, independently of Green operators or integration schemes. For models B and C, first, we consider both the maximum relative absolute error on local strain values along the interface and the relative absolute error on the maximum of local strains along the interface. The latter is relevant for investigating the influence of oscillations on damage initiation criteria. Meanwhile, the average relative error along the interface, in the fiber center as well as over the whole domain are also taken for comparison. All these relative errors are defined in Eq. (2.7). Clearly, W-S is the most accurate when compared with FEMR while SF is the most accurate when compared with FEMF. This confirms once again that the type of integration has an influence even for FFT methods. This is verified for both error measures. These two figures re-confirm that W-S is always more accurate when compared with FEMR, while SF is more accurate when compared with FEMF, except for the average relative errors at the fiber center.

We noticed that the W-S method produces more oscillations in the bulk (shown in Fig. 2.23), which however are less visible and can be neglected due to their sufficiently low values. Regarding the model C, the maximum and average relative errors are presented in Fig. 2.15 and 2.16 respectively. In general, the conclusions are the same as for model B. SF has the best results when compared with FEMF and it is W-S when compared with FEMR. Furthermore, Figs. 2.14 and 2.16 show that the average relative errors decrease when refining the mesh

Conclusions

Three FFT methods give identical results on model A. On models B and C, however, the SF is the most accurate when compared to full integration FEM and the W-S is more accurate when compared to the reduced integration formulation. Considering that SF has fewer oscillations (detailed in the next section) than others, in consequence, if a simulation needs high accuracy where a full integration is required in the FEM, the SF should be applied when the simulation is performed using an FFT solver.

However, applying SF requires eight times more memory than W-S and the computation time is also increased. In many cases, the full integration formulation is not needed, and the hourglass effect can be easily controlled by introducing "hourglass stiffness" as done in most FEM codes. Furthermore, the non-smooth interface has been proven to be the most critical cause of oscillations (detailed in the next section). Thus, a reduced integration formulation could be a better choice in many cases, where the W-S is, generally, more accurate than other FFT methods with an optimum computational speed.

In the next section, these three FFT methods as well as the FEM with both integration formulations are analyzed more thoroughly to investigate the causes of the oscillations.

Analysis of each cause of oscillations 2.6.1 Overview

The last part of our work in this chapter is an analysis of the oscillations to determine the impact of each cause quantitatively. In this section, the FFT methods and FEMs used in Section 2.5 are taken to make this evaluation. First of all, we need to establish a clear measure to quantify the oscillations.

Definition of the oscillations

To compare the influence of oscillations, the first step is to define the criterion to quantify them. Researchers often use the stress field to quantify the oscillations (e.g. [START_REF] Eloh | Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains[END_REF]). Even though the fiber and the matrix present the same amplitude of oscillations under some conditions, the oscillations in the fiber are often stronger than in the matrix because the fiber stiffness is normally higher than that of the matrix. Furthermore, when comparing different algorithms, researchers often put the color maps together and compare them visually. A better choice is to plot the curves of certain strain components of different methods along a direction. If there are some oscillations present in the curve, the curve will not be as smooth as it should be. This non-smoothness of the curve can be easily visualized and captured, and referred to as the oscillation degree of this curve. A detailed definition is presented in the following.

In Fig. 2.17, an illustrative reference curve without oscillations is shown as well as an oscillated curve. To quantify the irregularity of the oscillated curve (orange curve), we introduce a function D j standing for the differences in strain values along a curve. where ε j is the local strain at position j in the curve, and ε j +1 is the local strain value at the next position.

D j = ε j +1 -ε j , ( 2
Based on the Eq. ( 2.9), we get the D values of the oscillated curve and the reference. Generally, the more differences that two values have, the stronger are the oscillations, under the precondition that a reference curve is available. Therefore, we start by evaluating the smoothness of the FEMF method for model B. The values of D for half of the interface lines are shown in Figs. 2.18. We would expect the strain to increase increasingly intensely up to a maximum value at the fiber corner, and then decrease. Therefore, the values of D should be positive at the beginning and increase up to a maximum value, and then become negative at the fiber corner. This general tendency is confirmed in Fig. 2.18. The two large spikes are normal and correspond to the fiber corner where the maximum value is reached and D becomes negative.

However, we also see smaller spikes, which are more visible for the vertical interface line in Figs. 2.18(c,d). These can be considered as oscillations. Because their number stays nevertheless exceedingly small and they are localized around the fiber corners, it can be said that FEMF has almost no oscillation and can be considered as the reference.

Therefore, we can define the error measure as

MD j = |D FFTs/FEMR j -D FEMF j |, (2.10)
where |.| means the absolute value. Considering that the voxel position takes the same value independently of the method, the values of different methods can be directly compared.

If one method has fewer oscillations, its D j value should be closer to that of FEMF, which means that higher MD j values mean stronger oscillations, and vice-versa (as shown in Fig.

2.19

). Meanwhile, another function called MD is defined, which describes the average MD value along the oscillated curve. Lower MD values mean fewer oscillations in the strain variation curves.

Model and loading conditions

As introduced in Section 2.5.3, the mixed-type periodic boundary conditions are also applied in this section, with a macroscopic strain along the x-direction (transverse direction) of a value 〈ε〉 xx = 0.001 while stress-free conditions are imposed in the other directions (〈σ〉 y y = 〈σ〉 zz = 〈σ〉 x y = 〈σ〉 xz = 〈σ〉 y z = 0).

Based on the results of Section 2.5. As shown in Fig. 2.20, all methods show no oscillations for model A even for the original M-S algorithm, which confirms the results shown in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF]. As mentioned in the state of the art of this chapter, some researchers believe that one of the oscillations causes is the Gibbs phenomenon that describes the intrinsic defect of the spectral method when dealing with high contrast signals. This does not seem applicable to composites homogenization because there is no Gibbs phenomenon in model A.

Deep analysis of oscillations

To better visualize the oscillations in FFT and FEM, the local strain components are plotted for models B and C in Fig. 2.21, where the oscillations can clearly be seen.

For model B, Fig. Clearly, FFT methods based on the DGO reduce the oscillations significantly. Algorithm SF appears as the best choice in terms of oscillations reduction. W-S algorithm takes the second place and the M-S method is the worst. Furthermore, it should be paid attention that the FEMR gives remarkably comparable results as W-S while the oscillations are almost invisible on FEMF, where a full integration method is applied. This phenomenon confirms that the hourglass effect is also one of the causes of oscillation in FFT methods, as already reported in [START_REF] Leuschner | Fourier-Accelerated Nodal Solvers (FANS) for homogenization problems[END_REF]. This is because, like in the reduced integration FEM, the strains and the stresses calculated by FFT solvers are also computed at the centroid of the voxel, except for the SF method. This also explains why the latter shows fewer oscillations.

In general, the difference between M-S and W-S is the improvement of the finite-difference technique, which could be due to the non-satisfaction of Shannon's theorem, and the difference between W-S and FEMF is mainly due to the hourglass effect. However, as shown in Figs. 2.23(a,b) in black circles, the W-S method shows much more oscillations than others at the fiber center which could invite to a more advanced investigation. Fig. 2.22 also shows that the oscillations are reduced when refining the mesh, because refining the mesh increases the sampling frequency. It proves that the loss of high frequencies when discretizing the Fourier transform to a finite domain led to oscillations, which are reduced using a DGO.

Unlike models A and B, where both FFT methods and FEM rely on voxel meshes that are conformal, model C cannot rely on a mesh that is both voxel-based and conformal. Given that the non-smooth interface is one of the causes of oscillations, as shown in Fig. 2.24, all voxel meshes show strong oscillations. Thus, FEMF with a voxel-based mesh cannot be chosen as a reference anymore. Consequently, an FEM simulation with conformal mesh with full integration formulation is introduced to be used as a reference, so that Eq. (2.10) becomes

MD j = |D FFTs/FEMs j -D con f or mal j |.
Due to this choice of reference result, all calculations for model C do not share the same mesh, even for the same resolution. Therefore, a linear interpolation is applied on the conformal mesh to obtain the D j at the same positions as used in voxel meshes. In Fig. 2.25, the MD values for model C are presented. Herein, the radial strain component ε r r and the tangential strain component ε θθ are chosen, and D j is computed along the interface arc on the matrix side. This implies that the oscillations caused by non-smooth interfaces are much more critical than those originating from other causes such as the hourglass effect. Nevertheless, the oscillations present in model C do not seem to decrease with increasing resolutions (sampling frequencies). We did not investigate this further. Besides, it is not possible to explain why there are fewer oscillations with M-S in that regard. This is also left for future study.

Conclusions

The non-smooth interface is the most significant cause of oscillations that contributes more than five times than others. Note that non-smooth interfaces are more common in complex composite structures. The "truncated" Green operator could take the second place due to its higher impact than the hourglass effect but its influence on oscillations can be reduced significantly by using a discrete Green operator based on the finite-difference scheme. The hourglass effect takes the third place and could be overcome by using Schneider's discrete Green operator with full integration. However, Schneider's scheme means more computational resources occupation than Willot's scheme, and is therefore recommended only if a high precision is required. Finally, there is no evidence of the Gibbs phenomenon in our calculations, which proves the Gibbs phenomenon is not an issue at least for mechanical problems.

Chapter conclusions

In this chapter, we first have introduced a literature review of several causes that lead to oscillations in the FFT method. These causes can be divided into two categories. The first one regroups causes related to the intrinsic characteristics of the FFT, like non-satisfaction of Shannon's theorem. The second one regroups causes that can be also found for FEM with voxel-based meshes, like the hourglass effect or the non-smooth interface. Then, we have introduced some improved techniques proposed in the literature to reduce these oscillations. Some methods are easy to implement, such as low-pass filtering, but may lead to information loss. Meanwhile, there are also methods that work well, such as the consistent Green operator, but which are excessively demanding in terms of computational resources.

The idea of calculating a DGO based on finite-difference schemes seems to get a balance between the accuracy and the computational complexity.

Therefore, the improved techniques proposed by Willot [START_REF] Willot | Fourier-based schemes for computing the mechanical response of composites with accurate local fields[END_REF], and Schneider [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF] were chosen for deeper analysis as well as the original basic scheme. A quantitative comparison between FFT methods and FEMs was accomplished. Overall, the SF is the most accurate compared to the full integration FEM method and the W-S is more accurate compared to the reduced integration formulation. Consequently, if a simulation needs high accuracy where a full integration is required in the FEM solvers, the SF should be applied when the simulation is performed using an FFT solver.

However, applying SF means more computational time, which makes its application unessential in many cases. Besides, the hourglass effect can be easily controlled by introducing "hourglass stiffness" as done in most FEM codes. Thus, a reduced integration formulation can be a more optimum choice, where the W-S is, generally, more accurate than other FFT methods with an optimum computation speed.

Finally, we go back to the causes of the oscillations, where deeper analysis is performed to evaluate the impact of each of them. The non-smooth interface can be considered the most critical cause of oscillations. It is common in complex composite structures with voxel-based meshing. Therefore, an interface smoothing technique becomes an exigence to get cleaner local fields, as discussed in the next chapter.

Chapter overview

Among the causes of oscillations, the non-smooth is the most severe one. On the one side, it leads to local numerical oscillations at least five times stronger than others. On the other side, it is common in heterogeneous materials. Consequently, an interface smoothing technique is necessary to get more precise local fields.

As discussed in [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF], using Schneider's full integration scheme will cost around eight times the CPU time than Willot's scheme. Furthermore, because the hourglass effect has the weakest effect among all the causes, Willot's scheme is always chosen for this chapter's studies.

In this chapter, two interface smoothing techniques will be presented. The first one is the neighbour voxel average method (Section 3.2), which was initially proposed in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF] for a FEM solver and re-proposed for FFT solvers in this thesis [START_REF] Ma | Quantitative comparison between fast fourier transform and finite element method for micromechanical modeling of composite[END_REF]. Another one is the composite voxel method and, in particular, an enhanced composite voxel method proposed in this thesis (Section 3.3). A comparison between the different methods is performed in Section 3.4.

Neighbour voxels average method

Overview and description

The results of Doitrand et al. [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF] have shown strong oscillations in the FEM using a voxelbased mesh, and these oscillations seemed to be even stronger after mesh refinement, while there was no issue with a conformal mesh. That demonstrates that this issue is mainly due to the voxel-based mesh itself. The oscillations in FFT and FEM can also be called the "checkboard patterns" present in local strain and stress fields (as shown in Fig. 2.7). Clearly, the value in one voxel is larger than usual, while its neighbours have lower values. Based on the tests in Section 2.5.2, however, this does not affect the overall response. Therefore, to reduce this issue, a simple idea is to average the voxel values with surrounding voxels.

This idea of averaging was implemented in [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF] by calculating the stress at each integration point as an interpolation of the non-smoothed values at the surrounding integration points using a least-squares method. The degree of smoothing depends on the number of surrounding integration points that are taken into account. In [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF], all integration points that belong to elements that share a common node were taken into account. As an alternative, Fang [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF] proposed an approach by averaging the target voxel value with the neighbour voxels. In 2019, the present author proposed a similar method for FFT solvers in [START_REF] Ma | Quantitative comparison between fast fourier transform and finite element method for micromechanical modeling of composite[END_REF].

In general, the average equation in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF][START_REF] Ma | Quantitative comparison between fast fourier transform and finite element method for micromechanical modeling of composite[END_REF] can be presented as

σ = Ω σφ (r ) Ω φ (r ) , (3.1)
where σ is the Cauchy stress tensor, Ω is the averaging window, φ (r ) is the weight function, and r is the minimum distance from each surrounding voxel to the target voxel. As shown in Eq. (3.1), two important parameters should be defined: the size of the averaging window Ω, and the weight function φ.

As discussed in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF], for different voxel mesh densities and resolutions, it is difficult to determine a constant window size. In order to define this window size, we first define the notion of a layer. For a given voxel, layer L 0 is a set containing only this voxel. Then, layer L i +1 is the set of voxels whose boundaries intersect a voxel in layer L i . This intersection can be a face, an edge, or even a voxel corner. Note that all voxels chosen in layer L i are of the same phase as L 0 to ensure that the local fields may remain discontinuous across interfaces. The window size is defined as the layer count N L such that the neighbour voxels considered for the averaging are all voxels contained in the sets L i , i ≤ N L . A weight function can also be introduced to reduce the influence of voxels inside the window but far away from its center, the target voxel. In general, there are three types of weight functions (mentioned in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF]) that can be applied: constant, linear, and exponential functions, as φ

(L i ) = 1 (Constant), (3.2a) φ (L i ) = 1 - L i L max + 1 (Linear), (3.2b) φ (L i ) = exp L i L max (Exponential). (3.2c)
Therefore, the neighbour voxel average equation could be expressed generically as:

A orange = A orange φ (L 0 ) + VN L i j =1 A j φ (L i ) + ... + VN L max j =1 A j φ (L max ) 1 + VN L i * φ L i + ... + VN L max * φ L max , (3.3)
where A orange is the local value that needs to be averaged, and VN L i is the total voxel number of each layer taken into account. Note that a linear weight function is chosen in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF][START_REF] Ma | Quantitative comparison between fast fourier transform and finite element method for micromechanical modeling of composite[END_REF]. This is also the function that is chosen herein.

Test cases and conclusions

In this section, the neighbour voxel average with N L = 2 is applied on model C (presented in Fig. 2.8), where the fiber volume fraction is set as 55%, to test the effect of this approach. The materials properties are the same as in Table 2.1.

According to Fig. 2.25, the resolution N C = 255 is the one that shows the lowest oscillations. It is therefore chosen for the present study. The model is subjected to periodic boundary conditions with a macroscopic strain along the x-direction (transverse direction) of 〈ε〉 xx = 1.10 -3 and stress-free conditions in the other directions (〈σ〉 y y = 〈σ〉 zz = 〈σ〉 x y = 〈σ〉 xz = 〈σ〉 y z = 0).

As shown in Fig. 3.2 (b), the local stress fields show very strong oscillations without average, and the peak stress reaches over 30MPa while it is less than 25MPa in FEM results. Fig. 3.2(c) clearly shows the effect of averaging as the interface stress is smoother, and the peak stress is much closer to FEM results. In non-linear material modeling, the principal stresses and the von Mises stress often play a role in yield or fracture criteria. The results for these stresses are presented in Fig. 3.3. Once again, averaging reduces oscillations and compares well with FEM.

To determine the optimum layer count (N L ), a quantitative comparison is conducted. Our experiences show that the oscillations are reduced significantly when increasing the N L value at the beginning and become stable from N L ≥ 3. However, since the neighbour voxel average is only a post-processing operation, it is not coupled to the FFT solver, and there is no guarantee that the results will respect the equilibrium and compatibility equations. Therefore, it is preferred to keep N L as small as possible. For N L = 2, the abnormal stresses are sufficiently low. Thus, from our point of view, N L = 2 is the best choice for the neighbour voxel average method, which also confirms the conclusion in [START_REF] Fang | Smoothing artificial stress concentrations in voxel-based models of textile composites[END_REF] CHAPTER 3. INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTS

Composite voxel method

Overview

In a voxel mesh, the properties of each voxel are usually equal to the properties of the material, which is present at the center of the voxel [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF]. It is proper if the total volume of a voxel belongs to only one material. However, in the case of voxels shared by several phases and thus crossed by one or multiple interfaces, this kind of properties assignment is a rough approximation, and zig-zag interfaces are formed (shown in Fig. 3.4(a)). One of the reasons that limit the application of voxel-based meshes is that interfaces are not as smooth as in a conformal mesh that is closer to reality. More severely, as stated in the previous chapter, the zig-zag interface is one of the critical causes of oscillations near the interface both in the FFTs and the FEMs.

Apart from the neighbour voxel average method, which is a kind of post-processing, another idea is smoothing the interface by creating a new "phase" for the voxels crossed by the interface to model the interface. It should be noted that this composite voxel method is not related to the notion of interphase, which is not discussed nor modeled herein. The properties of these voxels are then homogenized from the local phases that share the voxels. This method is called the composite voxel method. As an illustration of this method, a 2D example is shown in Fig. 3.4(a), where two phases in blue and red are present, as well as a zig-zag interface. The orange line is the real interface. Meanwhile, in Fig. 3.4(b), apart from the blue and red voxels, it can be seen that all the voxels shared by two phases are shaded. They are called composite voxels, and their properties are equal to the homogenized properties of these two phases.

Conventional composite voxel method

The composite voxel approach was proposed by Kabel et al. in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF] and Gélébart in [START_REF] Gélébart | Filtering material properties to improve fft-based methods for numerical homogenization[END_REF].

The properties of a composite voxel are evaluated by a local homogenization based on the properties of phases that share the voxel and their volume fractions (shown in Fig. 3.4(b)).

CHAPTER 3. INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTS

Like the homogenization from microscopic to the macroscopic level, this local homogenization can also rely on several techniques: Voigt, Reuss estimates proposed in [START_REF] Gélébart | Filtering material properties to improve fft-based methods for numerical homogenization[END_REF] and Voigt-Reuss mean, or laminate theory as proposed in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF]. The first three can be expressed as:

C Voigt = N i =1 V ph i f C ph i (Voigt), (3.4a) C Reuss = N i =1 V ph i f C ph i -1 -1 (Reuss), (3.4b) C VR-mean = 1 2 C Voi g t + C Reuss (Voigt-Reuss mean), (3.4c)
where V

ph i f
and C ph i are the volume fraction and stiffness tensor of phase ph i inside the composite voxel, respectively. One drawback of the above mixing rules is that the orientation of the interface is neglected. In contrast, in the laminate mixing rule, the normal direction of the interface is taken into account. The formula of the laminate mixing rule is

P + K (C laminate -KId) -1 -1 = P + K (C -KId) -1 -1 , (3.5)
where C laminate is the homogenized elastic tensor of a composite voxel and K > 0 is a factor that should be chosen sufficiently large. The operator 〈-〉 is the volume average within the composite voxel. According to [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF], K should be larger than the largest eigenvalue of the stiffness matrices C (x) for all x in V. The fourth-order identity tensor is denoted as Id, and P is the fourth-order tensor that is expressed as follows:

P i j kl = 1 2 n i δ j k n l + n i δ j l n k + n j δ i k n l + n j δ i l n k -n i n j n l n m , (3.6)
where n is defined as the normal vector of the interface. Fig. 3.5 illustrates these three different mixing rules. In [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF], Hashin's model with inclusion (ph 1 ) and matrix (ph 2 ) was tested, and the results are shown in Fig. 3.6. When the inclusion is stiffer than the matrix, the Reuss rule gives the best results, while the Voigt one is the best in the opposite case. Meanwhile, the laminate rule gives a more balanced result. In [START_REF] Charière | Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites[END_REF], a combination of different mixing laws was studied, which, however, is not our purpose. Thus, the laminate rule is chosen herein. To calculate the volume fraction of each constituent in a composite voxel, a method called "sub-voxel approximation" has been proposed in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF]. As illustrated in Fig. 3 The last parameter that needs to be determined is the normal vector of a composite voxel (n). The method, connecting the center of mass of the dominant material with the center of the composite voxel (as shown in Fig. 3.8), was proposed in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF]. More precisely, assuming CHAPTER 3. INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTS W to be the composite voxel and S ⊂ W to be the dominant phase, the formula of the normal vector (n) is given as:

ñ := 1 |S| S xd x - 1 |W| W xd x, n = ñ ∥ ñ∥ , ( 3

Improved composite voxel method

One main drawback of the composite voxel method proposed in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF] is that a parametric representation of the interface and its geometry should be known in advance. When extracting a model directly from an image obtained from an SEM, where the real interfaces are not explicit nor parameterized, this approach can be difficult to apply. Besides, it is also difficult to describe the interface mathematically for some complex structures (e.g. textile composite and short fibers composite). Thus, for a composite voxel containing two or more phases, an approach mainly inspired from the level-set method is proposed hereafter to calculate the volume fraction of a phase and the normal vector.

In the conventional composite voxel method [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF], the identification of the composite voxel is determining whether a voxel contains more than one phase or not. In the proposed approach, instead of determining the number of phases in a voxel, a composite voxel zone is created by an artificial parameter, i.e. the composite voxel zone thickness l . All voxels with centers inside the zone are considered as composite voxels. Fig. 3.9 (red color voxels are the inclusion (ph 1 ), and blue ones are the matrix (ph 2 )) presents a simple example to describe this new approach for two phases. The steps are as follows:

• Identification (Fig. 3.9(a)): Identify the boundaries between each phase. Parametric representations of phases boundaries can be used if they are available. Otherwise,

CHAPTER 3. INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTS

which is more common, the zig-zag interface identification can easily be made with commercial software or algorithms.

• Calculation (Fig. 3.9(b)): For each phase ph i (ph 1 or ph 2 in the example), calculate the minimum distance d ph i from each voxel center to the phase boundary. Herein, a positive value is assigned for the voxels of the inclusion ph 1 and a negative value for the matrix ph 2 .

• Selection/assignment (Fig. 3.9(c)): Choose a composite voxel zone thickness value l . All voxels such that |d ph i | ≤ l 2 are considered as composite voxels (marked by a cross in the figure). of each phase ph i within the voxel:

Ṽph i f =              0, d ph i < -l 2 , 1 2   1 + 2d ph i l + sin 2d ph i π l π   , |d ph i | ≤ l 2 , 1, d ph i > l 2 .
If there are more than two phases, the real local volume fraction can be written as

V ph i f := Ṽph i f ph i Ṽph i f . (3.8)
Based on the properties of signed distance functions, the calculation of the normal vector (n) is given as:

ñ := ∇d P , n = ñ ∥ ñ∥ , with P = argmax ph i V ph i f .
(3.9)

In numerical calculations, the gradient (∇-) in Eq. (3.9) is calculated by classical centered finite differences. Finally, the homogenized properties of composite voxels are calculated using Eqs. (3.5) and (3.6). With this optimization, a parametric representation of the interface is not required, and there is no need for sub-voxels.

An added advantage of this enhancement is that it can easily be generalized to Schneider's full integration scheme. It is possible to compute the signed distance function at each integration point and then the local phase volume fractions and normal vectors. This can be interesting to smooth the transition over multiple integration points and reduce the oscillations. Due to the significant increase in computational cost when using eight integration points per voxel, we restricted our composite voxels analysis to Willot's scheme.

Comparisons of interface smoothing technique

In this section, three interface smoothing techniques, implemented in FFT solver with W-S scheme, will be compared in terms of the effect of oscillation reduction and the relative error to the results of FEM of reduced formulation with conformal mesh. These three methods are:

• Conventional composite voxel model proposed in [START_REF] Kabel | Use of composite voxels in FFT-based homogenization[END_REF] (TCV )

• Our optimized composite voxel model (OCV )

• Neighbour voxel average method (AVE)

For the improved method OCV , it is necessary to define the composite voxel thickness zone l . To be consistent with the conventional method TCV , we choose to set l as twice the voxel size, which ensures that there is always at least one layer of the composite voxel at the interface, and at most two layers (the one on the matrix side and the other on the fiber side). That ensures that all voxels at the interface are included in the composite voxel zone.

As for the AVE technique, the most suitable average layer count is N L = 2 as proven in Section 3.2.2. Thus, N L = 1 and 2 are both studied in the following.

In this part, the model C shown in Fig. 2.8(c) is chosen, and subjected to a macroscopic strain along the x-direction (transverse direction) with a value 〈ε〉 xx = 0.001 while stress-free conditions are imposed in the other directions (〈σ〉 y y = 〈σ〉 zz = 〈σ〉 x y = 〈σ〉 xz = 〈σ〉 y z = 0). These mixed-type loading conditions are applied using the method presented in [START_REF] Kabel | Mixed boundary conditions for FFT-based homogenization at finite strains[END_REF]. It is reminded that periodic boundary conditions are intrinsic to FFT methods. Due to stressfree boundaries, Poisson's effect is active during loading. Meanwhile, the resolutions chosen for model C are N C =225, 255, 305, 345, 375 and 455, and the quantitative definition of the oscillations is the same as in Section 2.6.2.

The first layer of voxels along the interface is occupied by the composite voxels, which correspond to neither matrix nor fiber material. Therefore, it is not possible to compare the values in this layer between the different techniques. That is why the first layer of voxels out of the composite voxels zone, denoted as interface 2nd to distinguish from the real interface, is chosen to make the comparison in Figs. 3.10, 3.11, 3.12, 3.13, 3.14, and 3.15. Since the composite voxels zone is very small compared to the model size, the layer of interface 2nd can also be considered as the interface response. As shown in Fig. 3.10, where the interface 2nd normal strain (ε r r ) and tangential strain (ε θθ ) are drawn, the composite voxel method can reduce significantly the oscillations in the presence of a non-smooth interface. Note that the detailed description of ε r r and ε θθ is introduced in Eq. (2.8) and in Fig. 2.11. For the TCV method, we tried using different numbers of sub-voxels. The response is converged at N sub = 1681. Since model C can be treated as a 2D model, that means a 41x41 division in sub-voxels in x x y direction. With our unoptimized Python implementations, the total computation time for the TCV method using N sub = 1681 is of the same order as CHAPTER 3. INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTS the one for the OCV method. More importantly, it is interesting that OCV reduces the oscillations more effectively than TCV for both strain components and all resolutions, and is also more general because it can easily deal with non-parametrized interfaces.

To compare OCV and AVE, the MD values are shown in Figs. 3.12. These two graphs show that OCV reduces the oscillations more effectively than the AVE if one layer is used whereas with two layers AVE becomes more effective. However, as illustrated in Fig. 3.13, the AVE cannot treat well the abnormal phenomenon on the normal interface 2nd strain at around 45 degrees, while the OCV does not have this kind of issue.

Besides, AVE is a kind of post-treatment, where the equilibrium equation and the compatibility condition are not verified. If a very high reduction of oscillations is required, a combination of these two methods can be applied. This possibility is not discussed herein, and the OCV seems to be a more optimal choice. In order to validate these interface smoothing techniques, we investigate the relative errors after applying these treatments. The reference is also the FEM of reduced formulation with conformal mesh.

The relative error for model C is first checked and shown in Fig. 3.14. For the normal strain error, apart from an abnormal result at N = 251 that has not been investigated yet, the different smoothing techniques effectively reduce the error and give comparable results. Regarding the tangential strain error, Fig. 3.14.(b) shows that the composite voxel method is still effective but not the neighbour voxels averaging. In addition, our proposed optimization of the composite voxel method does not deteriorate nor improve the accuracy significantly.

The second part is to verify the influence of the OCV technique on the two other models. For model A, the influence of OCV is null except for the voxels inside the composite voxel zone because of the homogenization in this zone. Regarding model B, the relative errors are shown in Fig. 3.15. The errors for M-S are also reported for comparison. We observe that AVE 2L increases the errors significantly making this only relevant for non-smooth interfaces and thus does not apply to general situations. The OCV technique also increases a little the relative error, which is quite normal because adding composite voxels artificially reduces the contrast between fiber and matrix. Nevertheless, the relative error of the OCV method is much lower than the AVE 2L method. Furthermore, it is very close to the M-S scheme.

In conclusion, the increase of the error for models A and B is acceptable for the improved composite voxel method that contributes to the universality of this approach in a general application.

Chapter conclusions

Figs. 3.2, 3.3, and 2.24 show that the zig-zag interfaces can lead to stronger non-physical oscillations and peak stresses than FEM conformal mesh. These strong and non-physical artifacts are problematic for damage and fracture modeling.

Therefore, this chapter has investigated two techniques to improve the zig-zag interfaces due to the voxel-based discretization. They are

• the neighbour voxels average method,

• the composite voxel method.

The first investigation is on the neighbour voxels averaging (AVE). Figs. 3.2 and 3.3 show that it reduces the oscillations significantly and leads to peak stresses much closer to FEM conformal mesh results, as well as smoother local fields. Regarding the layer count (N L ) used for averaging, our results show that the oscillation reduction becomes less effective when N L > 2 and almost non-effective when N L > 3. To minimize the computational cost, N L = 2 is the best choice for the neighbour voxels average method.

The second part is our study on the composite voxel method. An improved composite voxel method (OCV ) is proposed in this chapter. Compared with the conventional composite voxel method TCV , the OCV does not require parametric representations of interfaces and can hence be applied directly to digital images. Additionally, our results show that this simplified method is still at least as effective as the TCV in terms of oscillations reduction and accuracy.

When comparing the OCV with the neighbour voxels average method AVE, it is observed that OCV reduces the oscillations more effectively than AVE with one layer (N L = 1) but less effectively than AVE with two layers (N L = 2). In addition, we noticed an increase in the error when there are no zig-zag patterns for all these smoothing techniques. This increase is more significant when using AVE with two layers compared to OCV . Besides, Fig. 3. [START_REF] Böhm | A short introduction to continuum micromechanics[END_REF] shows that the AVE cannot treat well some abnormal phenomenons, while the OCV does not have this issue.

As a conclusion, both AVE and OCV have good oscillation reduction effect. On the one hand, since AVE is non-intrusive and easier to implement, even if it has some drawbacks on the accuracy, it can be used when only post-processing is possible. On the other hand, OCV shows less relative error than average methods, which makes it more suitable for complex analysis. It additionally ensures that the computed smooth fields satisfy equilibrium and compatibility equations.

Up to now, we have only discussed elastic behavior modeling using an FFT solver. In the next chapters, damage modeling in heterogeneous materials using the phase-field method and an FFT solver will be discussed.

Chapter overview

Structural failure under loading is one of the most severe problems as it can lead to unpredictable losses and casualties. Therefore, failure study is one of the most important topics nowadays in material and structural engineering. A general introduction of different failure and damage modeling techniques, as well as the phase-field method, is presented in Section 4.2 of this chapter. Then, a detailed mathematical description of the phase-field method is shown in Section 4.3.

Damage modeling of material 4.2.1 Phenomenological modeling

In the engineering domain, phenomenological failure criteria have been used for decades to estimate the feasibility of designs quickly. The maximum stress/strain criterion is undoubtedly one of the most straightforward criteria. It predicts failure when critical conditions in loading are larger than the strengths (tension, compression, shear) or maximum allowable strain (tension, compression, shear). However, real loading conditions are often very complex, especially for composite materials, and the maximum stress/strain criteria do not model mixed loading conditions. These can be modeled with the Tsai-Hill criterion.

Due to the quadratic form of the Tsai-Hill criterion, it can not distinguish between tension and compression. This is the reason why the Tsai-Wu criterion has been proposed in [START_REF] Tsai | A General Theory of Strength for Anisotropic Materials[END_REF].

For composite materials, a 3D failure criterion proposed by Hashin in [START_REF] Hashin | Failure Criteria for Unidirectional Fiber Composites[END_REF][START_REF] Hashin | Fatigue failure criteria for unidirectional fiber composites[END_REF] became more and more popular because it can predict the failure both for different constituents (matrix/fiber) and directions.

Fracture mechanics modeling

The failure criteria presented above can also be used to predict crack initiation, and they can be completed with a suitable crack propagation criterion. This subject has been extensively studied since the milestone work of Griffith [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] and Irwin [START_REF] Irwin | Analysis of Stresses and Strains Near the End of a Crack Traversing a Plate[END_REF], who established the theory of Linear Elastic Fracture Mechanics (LEFM). In the work of Griffith, a global energy approach is provided to find quantitative relations between the crack length, the resistance to the crack growth of the materials, and the criterion of crack propagation that regards fracture as a competition between the surface energy of propagation and the elastic energy stored in the bulk material. Griffith proposed the relation:

σ a = 2Eτ π , ( 4.1)
where a describes the crack length, and τ is the surface energy density. In general, the energy release rate G can be written as G = 2.τ for brittle materials. Under the critical case, where σ = σ c , we can get the critical energy release rate formulation:

G c = σ 2 c aπ E . ( 4.2)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Decades later, Irwin [START_REF] Irwin | Analysis of Stresses and Strains Near the End of a Crack Traversing a Plate[END_REF] introduced the Stress Intensity Factor (SIF), which is widely used in engineering and can be expressed as

K I = σ πa. (4.3)
In the case of mixed loading, the relation between G and K i is

G = K 2 I E ′ + K 2 II E ′ + K 2 III 2µ (4.4)
where E ′ = E for plane stress and E ′ = E 1-ν 2 for plane strain. When G > G c or K > K c , crack propagation is predicted. Based on the above equations, it can be seen that a pre-existing crack is needed, which means only crack propagation can be predicted with these theories.

With LEFM, it is possible to predict if the crack will propagate and also the crack propagation direction. The numerical implementation of this theory is quite complex as introducing and propagating a crack poses a great challenge in a mesh-based technique such as a standard finite element approach with a continuous displacement field. The displacement jump due to the crack can be captured by two methods in general.

The first one uses standard finite element polynomial approximation and relies on meshes that conform to discontinuities which requires duplicating nodes, and reconstructing or adapting the mesh [START_REF] Bouchard | Crack propagation modelling using an advanced remeshing technique[END_REF].

The second one is called enriched FEM and consists in enriching the polynomial approximation space so that the discontinuities can be modeled independently of the mesh, which is an attractive feature as no mesh modification is needed [START_REF] Fries | The extended/generalized finite element method: An overview of the method and its applications[END_REF]. This enrichment can be achieved, on the one hand, by adding special shape functions to the standard finite element (extrinsic enrichment); or on the other hand, by modifying some of the shape functions (intrinsic enrichment). The most popular of these enriched methods are the XFEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] and the GFEM [START_REF] Strouboulis | The design and analysis of the generalized finite element method[END_REF], where the main difference is that the XFEM only involves local enrichment while the GFEM enriches all the nodes in the discretization [START_REF] Shakoor | Computational methods for ductile fracture modeling at the microscale[END_REF]. One challenge for these enriched methods is that since the crack is not modeled by the mesh, it must be represented in a separate mesh or using implicit representations like signed distance functions. This can become a great challenge for complex crack propagation patterns, especially in 3D.

Cohesive Zone Model (CZM)

One of the drawbacks of LEFM described above is that a singularity is present at the crack tip, which leads to an infinite stress at this point. This is not physical and raises a problem for numerical modeling as any crack propagation criterion would be instantaneously satisfied whichever the loading with a sufficiently fine mesh at the crack tip.

Fortunately, under the framework of fracture mechanics, it is possible to introduce a Cohesive Zone Model (CZM) to remove this singularity. This approach was pioneered by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and Barenblatt [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF]. As shown in Fig. 4.1, the singularity at the crack tip is embedded in a process zone which corresponds to a prospective surface that can be separated under loading. In [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF], the forces that prevent the crack separation were related to atomic or molecular attractions. The creation of a new crack surface was interpreted as a breaking of atomic and molecular bonds in the material. Nowadays, the cohesive zone model is more and more used for simulating fiber/matrix debonding or yarns/matrix separation in composite materials [START_REF] Zhang | A novel interface constitutive model for prediction of stiffness and strength in 3d braided composites[END_REF][START_REF] Herráez | Transverse cracking of cross-ply laminates: A computational micromechanics perspective[END_REF].

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

A frictionless interfacial cohesive zone model is detailed in the following. Supposing an initial interface stiffness {k 1 , k 2 , k 3 }, the constitutive relationship at the interface can be described as:

   t 1 t 2 t 3    =   k 1 0 0 0 k 2 0 0 0 k 3      δ 1 δ 2 δ 3    . ( 4.5)
The softening initiation displacement, δ 0 1 (Traction), δ 0 2 (Shear), δ 0 3 (Shear), can be defined as:

δ 0 1 = N k 1 , δ 0 2 = S k 2 , δ 0 3 = T k 3 . ( 4
δ f 1 = 2G IC N , δ f 2 = 2G IIC S , δ f 3 = 2G IIIC T . (4.7)
With the help of δ 0 i and δ f i as well as the strengths, it is easy to establish the traction-separation law for each direction. Generally, to simplify the calculation, in the literature [START_REF] Zhang | A novel interface constitutive model for prediction of stiffness and strength in 3d braided composites[END_REF], it is usually assumed that k 1 = k 2 = k 3 = k, S = T, and G IIc = G IIIc . Under mixed loading, the effective relative displacement δ m can be easily expressed as the norm of the vector of displacement under different directions:

δ m = 〈δ 1 〉 2 + δ 2 2 + δ 2 3 = 〈δ 1 〉 2 + δ 2 shear , (4.8)
with 〈x〉 expressed as 〈x〉 = 0, x ≤ 0

x, x > 0 (4.9)

The damage evolution function of interface d is then given by

d = δ f m δ max m -δ 0 m δ max m δ f m -δ 0 m , d = [0, 1] (4.10)
where δ max m = max δ max m , δ m . The δ 0 m and δ f m denote the initial failure and final failure displacement of the surface, and can be obtained by adopting a quadratic stress criterion and the power interaction law of the energy, as detailed in [START_REF] Metoui | Separated représentations for th multiscale simulation of the mechanical behavior and damages of composite materials[END_REF][START_REF] Zhang | A novel interface constitutive model for prediction of stiffness and strength in 3d braided composites[END_REF][START_REF] Camanho | Mixed-mode decohesion finite elements for the simulation of delamination in composite materials[END_REF]. In Eq. (4.10), d = 0 means no interfacial damage while the interface will be totally damaged if d = 1. This damage parameter only describes the failure state of the interface, which is independent of the bulk material properties. However, the cohesive zone model only deals with the separation of the crack surface hence it should incorporate extra approaches, such as LEFM, to determine the propagation of the crack. Considering the exigence of a pre-defined crack for LEFM, from the fundamental point of view, LEFM and CZM require criteria to determine when/where a crack will initiate, grow, and in which direction.

The application of cohesive zone models in standard FEM usually involves inserting cohesive elements into the mesh. However, as most traction-separation laws have an initial part where the force increases monotonically with the displacement until a maximum value,

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

these insertions can modify the global stiffness of the material. This issue can be overcome by only allowing element opening once the critical criterion is achieved or inserting cohesive elements dynamically into the mesh during the simulation when the criterion is satisfied. This second approach requires a non-negligible modification of FE codes [START_REF] Shakoor | Computational methods for ductile fracture modeling at the microscale[END_REF]. This approach, however, has the advantage of solving another issue of cohesive elements, which is mesh dependence. Indeed, if cohesive elements are pre-inserted into the mesh, the crack can only propagate along element boundaries, which may lead to a non-physical zig-zag pattern.

Finally, implementing cohesive zone models for studying interface debonding involves fewer difficulties, as cohesive elements can be inserted only along the interface.

Continuum Damage Mechanics (CDM)

Another discipline called Continuum Damage Mechanics has been developed. If we look at the microscopic scale, right before the moment the material attains its strength, the macrocrack does not propagate directly. There is first an accumulation and evolution of local micro-defects, such as micro-cracks and micro-cavities. Because of these micro-defects, the surface for bearing the load reduces, and the material can not bear the load as before even if the effective macroscopic load increases on the material. From a global point of view, this phenomenon can be regarded as a degradation of local material properties. The degree of material degradation can be described by a continuous parameter d , which is called the damage variable. This variable d varies from 0 to 1, where 0 means intact state and 1 means fully degraded state. In [START_REF] Kachanov | Time of the rupture process under creep conditions[END_REF], the relation between stress (σ) and effective stress (σ eff ) can be written as:

σ eff = σ 1 -d . ( 4

.11)

In fact, the formula (1d ) can be generalized as a degradation function g (d ) that varies from 1 to 0 when d increases from 0 to 1. This degradation function or evolution law can be empirical or phenomenological [START_REF] Shakoor | Computational methods for ductile fracture modeling at the microscale[END_REF]. Under the assumptions of CDM, the material is still continuous under the degradation process, and the local mechanical fields are also continuous, which eases the implementation into standard finite elements in contrast to LEFM. In fact, under this continuous approach, damaged regions can naturally grow, branch, coalesce without any numerical difficulty. In general, numerically, the evolution and propagation of damage can be divided into four steps:

• Elastic part: The model is under mechanical loading and the damage initiation criterion is not satisfied, although local fields concentration point(s) may appear due to the micro-defects or the geometry.

• Damage initiation: When one or several points reach the damage initiation criterion, the damage variable becomes non-zero (d > 0) and can be calculated from a certain formulation.

• Damage evolution: With the loading increasing continuously, the damage at these points increases, and the stress fields decrease with the application of the degradation function until (d = 1), where the local stress fields vanish. An important property is that d cannot decrease.

• Damage propagation: When the local stresses fields at these points vanish, these points can not bear the load, which is transferred to other points. Consequently, the damage initiation criterion may then be satisfied at these other points. Progressively, the damage field will continue updating until the final failure of the model.

In [START_REF] Wang | Progressive damage analysis of 3D braided composites using FFT-based method[END_REF][START_REF] Matveev | Modelling of textile composites with fibre strength variability[END_REF], a simple damage model was applied for textile composites. The damage model utilized in these two papers is summarized hereafter. First, damage initiation criteria for yarns and matrix are established, where yarns utilize the maximum principal stress criterion while a modified von Mises yield criterion is applied for the matrix. Then, a value I can be calculated from these formulations to describe the difference between the actual load and damage criterion, if I < 1, it means the actual load is below the strength while I = 1 means damage begins, and I > 1 means damage evolving. With the values of I, the damage variables (d i ) of yarns and matrix in different directions can be calculated by

d i = 1 exp (-c 1 I i ) + c 2 , g (d i ) = 1 -d i , (4.12)
where parameters c 1 and c 2 are empirical constants. Generally, c 2 /c 1 = 1.62 can give close agreement with experimental tests. In [START_REF] Wang | Progressive damage analysis of 3D braided composites using FFT-based method[END_REF], constants c 1 = 8 and c 2 = 13 have been used. Because the d i in Eq. (4.12) can exceed 1, the degradation of elastic properties (E) can be expressed as: E = max (0.001, min (d i )) E 0 , (

where elastic properties with a superscript 0 correspond to the undamaged properties of constituents.

One of the drawbacks of this approach is that the numerical solution largely depends on the mesh size because the energy dissipated decreases and the local fields are increased with mesh refinement. Consequently, it is common to introduce a characteristic internal length l , which is a measure for the size of the localized area. The softening response after damage initiation is then characterized by a stress-displacement response rather than a stress-strain response. The crack band model falls under this category. It was first proposed in [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF], then followed by Lapczyk et al. [START_REF] Lapczyk | Progressive damage modeling in fiber-reinforced materials[END_REF], and Fang et al. [START_REF] Guo-Dong | Progressive damage and nonlinear analysis of 3d four-directional braided composites under unidirectional tension[END_REF]. In this method, the dissipated energy within an element of the mesh is set equal to the material fracture energy when the element has failed, which can be written as follows:

1 2 ε f eq σ f eq l 3 = G c l 2 (4.14)
where l is the characteristic length of the element, defined as the cube root of the element's volume in the case of cubic voxel elements, and ε f eq and σ f eq are the failure equivalent strain and equivalent peak stress. The equivalent displacement will vary with element size, and is defined as δ eq = ε eq l . (4.15)

Combining with Hashin's 3D damage initiation criterion (detailed in [START_REF] Zhang | A novel interface constitutive model for prediction of stiffness and strength in 3d braided composites[END_REF]), the evolution of the damage variable for each direction can be written as: In Eq. (4.16), δ i ,eq is the current equivalent displacement defined as in Eq. (4.15) and detailed in [START_REF] Guo-Dong | Progressive damage and nonlinear analysis of 3d four-directional braided composites under unidirectional tension[END_REF] for different directions. The displacement δ 0 i ,eq is the initial equivalent displacement at which the failure criterion is satisfied, and δ f i ,eq is the full equivalent displacement at which the material is completely failed. They are defined as:

d I = δ f i ,eq δ i ,eq -δ 0 i ,eq δ i ,eq δ f i ,eq -δ 0 i ,
δ 0 i ,eq = δ i ,eq / I i , (4.17)
δ f i ,eq = 2G I / σ i ,eq / I i . (4.18)
Here, I i is the value of the damage initiation criterion, G I is the critical energy release rate under mode I, and σ i ,eq is the current equivalent stress for a failure direction which is detailed in Table 1 of [START_REF] Guo-Dong | Progressive damage and nonlinear analysis of 3d four-directional braided composites under unidirectional tension[END_REF].

To assemble different damage variables into a matrix that describes a 3D damage field of the model, the Murakami damage model [START_REF] Murakami | Mechanical modeling of material damage[END_REF] can be adopted.

Non-local integration damage model

Apart from the crack band theory, another approach to solve the mesh dependency is the non-local damage model. Unlike the local model, the local fields at a point depend not only on the parameters at that point but also on the local fields in a finite neighborhood of that point. The solution may then depend on the size of the "neighborhood", but no longer on the mesh size.

The first idea is obviously taking an average over this finite neighborhood. This kind of work was pioneered by Baz ǎnt and Jir ásek [START_REF] Bazant | Nonlocal Continuum Damage, Localization Instability and Convergence[END_REF]. The non-local variable r is defined within a FE domain Ω by: r (x) = B(x) r y w xy dy

B(x) w x -y dy = 1 B (x) = {y ∈ Ω, |x -y| < l c }        , ∀x ∈ Ω (4.19)
where r can be damage variable or other local fields, and the size of the "neighborhood" is defined by a characteristic length l c . The weight function w xy decreases monotonically with increasing distance |x -y| to reduce the influence of points far away from x.

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Gradient-Enhanced Damage (GED) model

In order to reduce the computational cost due to the neighborhood search and non-local integration formulation, a so-called Gradient-Enhanced Damage model (GED) was proposed in [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]. If the weighting function w is chosen as a Gauss function, and r is chosen as local equivalent strain ε eq , the non-local integration formulation Eq. ( 4. [START_REF] Hopkins | Unique set of micromechanics equations for high-temperature metal matrix composites[END_REF]) can be easily rewritten in the form of a Taylor expansion, and the non-local equivalent strain ε eq can be approximated by solving a partial differential equation instead of computing an average. By omitting the higher order terms, an explicit form of non-local gradient formulation can be obtained:

r = r + l 2 c ∆r . (4.20)
First, high order gradient terms have been neglected in the explicit model. Second, the explicit formulation is said to be weakly non-local because the non-local strain in a point depends only on the local strain and its gradients at the same point. Third, the explicit gradient formulation imposes stronger continuity requirements on the displacements because of the gradient term [START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF]. To avoid these drawbacks, the implicit non-local formulation

r -l 2 c ∆r = r, ( 4.21)
is often preferred.

The gradient-based model has been the object of many propositions and novelties, such as in [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF][START_REF] Lorentz | A nonlocal damage model for plain concrete consistent with cohesive fracture[END_REF], which will not be discussed here. Alternative non-local damage models such as the thick level-set method [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] will not be presented here either.

The phase-field model, a modeling technique between fracture and damage mechanics, is attracting increased attention in the research domain. This is also one of the centers of this Ph.D. thesis. On the one hand, the phase-field can be under the category of gradient-based model. On the other side, it is closer to fracture mechanics. Therefore, it will be introduced in the following as an independent section.

Phase-field model

Similar to the CDMs, the phase-field models also utilize a damage variable d , which varies in the range [0, 1] to describe the failure state. The powerful idea is that this damage variable is also used to describe the cracks. Sharp cracks are regularized as diffusive crack bands by a function of d [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF]. On the one hand, the evolution of the phase-field variable itself completely describes the crack and the damaged regions of the material. On the other hand, phase-field models are closely related to the variational approach to brittle fracture [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF]. This approach enables to simultaneously get the displacement field and the cracks by minimizing the total potential energy of the cracking solid.

Table 4.1 shows the mathematical equations of the phase-field [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and the GED [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] models. Since they have similar mathematical formulations, it is natural to compare them, like the work in [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF][START_REF] Mandal | Phase field and gradient enhanced damage models for quasi-brittle failure: A numerical comparative study[END_REF]. The results of [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF][START_REF] Mandal | Phase field and gradient enhanced damage models for quasi-brittle failure: A numerical comparative study[END_REF] have shown that these two models are similar but also different. The similarities and differences can be summarized as two parts: the mathematics and the fundamental idea. 4.1: The PDE equation of gradient-enhanced damage model [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] and phase-field model.

Regarding the mathematics, based on Table 4.1, several conclusions can be raised:

• Left part of the equation: On the one side, in the phase-field model, the operator is applied to the damage variable itself. On the other side, the non-local variable r can also be the non-local damage variable in the GED [START_REF] De Borst | Gradient damage vs phase-field approaches for fracture: Similarities and differences[END_REF]. In a sense, they can be similar.

• Right part of the equation: In the phase-field model, when d = 1, its right part totally vanishes, which ensures a constant bandwidth of the crack. However, it is not the case in the GED models which results in damage widening.

• Considering the classical phase-field [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and GED [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] models, their damage profiles are different. As shown in Fig. 4.4, the damage diffusion in phase-field has an exponential form while it has a Gaussian form in the GED models.

.4: Typical damage profile of (a) phase-field model in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] and (b) GED model in [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF].

Recently, gradient damage models have been proposed that tend to be closer to the phasefield model, such as [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF]. They compare well in terms of equations and solutions. That is why the researchers consider the phase-field model as a gradient-type non-local damage model. However, their fundamental ideas are different. The gradient-type non-local damage model was proposed to average the local damage field to solve the mesh dependency, while the phase-field departs from the discontinuous description of the crack, where a distribution function is introduced to smear this discontinuity. For example, the regularization from the sharp crack to the diffusive crack is not introduced in [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF] even if the final equations are almost the same as the phase-field. More details can be found in the review paper [START_REF] Wu | Phase-field modeling of fracture[END_REF].

Despite the fundamental idea, generally, the phase-field method can be described as a kind of non-local gradient type model that uses a damage field d to regularize a sharp crack by a diffusive crack band, and this damage variable can be solved by minimizing the Griffith energy principle. The major characteristics of phase-field can be summarized as follows:

• The phase-field model is purely based on energy minimization and there is no need for a pre-defined crack, hence, the crack initiation, growth, and coalescence can be automatically determined.

• Multi-cracks merging and branching problems can be solved without additional effort.

• The model can easily be applied in 3D cases, and numerical implementation is straightforward.

• The model is non-local and intrinsically avoids mesh dependency.

In contrast, computational cost is one of the most severe drawbacks of phase-field, which can be dealt with using parallelization. The FFT method is easy to parallelize, hence, combining the phase-field with the FFT solver seems very interesting. This is the research track pursued in this thesis.

Many applications of the phase-field model can be found in the literature: [START_REF] Nguyen | Contributions to multiscale modelling of quasi-brittle damage in heterogeneous materials[END_REF] studies the damage and fracture behavior of quasi-brittle heterogeneous materials; T.H.N. Nguyen Hun [START_REF] Hun | Fracture modeling in clay materials under hydric shrinkageModélisation de fissure dans les matériaux argileux sous retrait hydrique : numerical models, comparisons with experiments and stochastic aspects[END_REF] studied crack propagation in clay materials under humidity to improve the understanding of degradation mechanisms in civil engineering structures; Rabette [START_REF] Rabette | Prédiction de la microfissuration par champ de phase et méthode FFT pour les matériaux énergétiques comprimés[END_REF] studied the micro-cracks inside compressed energetic materials (e.g., explosive powders), and B.T.

Vu [START_REF] Vu | Modeling and simulation of damage in anisotropic materials by the phase-field method[END_REF] studied the damage behavior of anisotropic composite materials, like, graphiteepoxy composites.

In the next section, a detailed mathematical description of different phase field methods is presented to explore the theoretical aspects of phase-field.

Phase-field theoretical aspects

Griffith theory

The phase-field model is based on Griffith's theory [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF]. In this energetic approach, the damage evolution and crack propagation are a consequence of the competition between the bulk energy stored in the body and the energy dissipation from the opening of the crack surface. From this point of view, in quasi-static loading, the total energy Π can be expressed as:

Π := Φ s + Φ d -P, (4.22)
where Φ s is the strain energy stored in the cracked body, Φ d is the energy dissipated for opening the crack surface, and P is the external loading. We consider a domain

Ω ⊂ R n (n = 1, 2, 3) that contains a crack set Γ with Γ ⊂ R n-1 ,
Φ s = Ω ϕ (ε (u) , Γ) dV, (4.23a)
Φ d = Γ G c dS, (4.23b)
P = Ω f * • udV + ∂Ω t t * • udS + ∂Ω u (σ • n) • u * dS, (4.23c)
where ϕ is the elastic strain energy density stored in the cracked body, G c is the critical energy release rate, σ and n are the stress tensor and the normal vector outside the boundary, and the dV and dS correspond to the volume and surface integration respectively.

The displacement field and the crack set can be solved by minimizing the total energy with the variational approach [START_REF] Bourdin | The variational approach to fracture[END_REF]:

(u (t ) , Γ(t )) = Arg {min [Π (u, Γ)]} (4.24)
under the constraint of irreversibility, which is stated as Γ ≥ 0, and under the boundary conditions: u(x) = u * at ∂Ω u , and σ • n = t * at ∂Ω t with n the outward unit normal vector to the external boundary ∂Ω, (∂Ω = ∂Ω u ∂Ω t ).

Regularization of sharp crack

The solid Ω with crack set Γ is presented in Fig. The influence of l c on the damage band diffusion is presented in Fig. 4.6. The smaller l c is, the narrower the crack band is. If l c = 0, a sharp crack is fully recovered. In most works, the function used to describe the diffusive crack state is one of: Miehe [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF], Pham [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF], and Wu [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF]. Actually, Pham's function was initially proposed for a gradient damage model. As explained in Section 4.2.5, his formulation is not much different from the phasefield method. Hence, Pham's method is included in this review but it is not discussed in details. The three diffusive crack functions are expressed as subject to the Dirichlet-type boundary conditions: d (x = 0) = 1 and d (x = ±∞) = 0. In consistence with the variational approach to fracture, the crack surface functional can be obtained from the minimization of the crack surface that yields the ODE shown in Eq. (4.26):

d M (x) = exp - |x| l c , Miehe, (4.25a
d (x) = arg { inf [Γ l ]} ,
with Γ l expressed as

Γ l (d ) = +∞ -∞ 1 2 1 l c d 2 + l c d ′ 2 dx = +∞ -∞ γdx,
where γ is the crack surface density. When the value of l c approaches zero, the regularized crack surface Γ l converges to a sharp crack. Meanwhile, the extension to higher dimension is straightforward:

γ (d , ∇d ) = 1 2 1 l c d 2 + l c (∇d ) 2 , ( 4.27)
with d ′ → ∇d . The crack surface density function in Eq. (4.27) depends on the diffusive crack function d , but this definition is not unique. In 2017, Wu [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF] proposed a generic form of the crack surface density function:

γ (d , ∇d ) = 1 c 0 1 l c α (d ) + l c (∇d ) 2 , (4.28)
where α (d ) is called geometric crack function and satisfies α (d = 0) = 0 and α (d = 1) = 1. Factor c 0 is the scaling factor which ensures the regularized crack at the fully damaged state represents a sharp crack [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF]. The expressions of α (d ) and c 0 are

α (d ) = ζd + (1 -ζ) d 2 , ( 4.29a
)

c 0 = 4 1 0 α (d )dd , (4.29b)
where ζ ∈ [0, 2] to guarantee α (d ) ∈ [0, 1]. In short, the α and c 0 values of the approaches presented in Eq. (4.23) are summarized in Table 4.2.

ζ α (d ) c 0 d (x)
Miehe [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF] There also exist other functions such as those proposed in [START_REF] Alessi | Gradient damage models coupled with plasticity: Variational formulation and main properties[END_REF] and [START_REF] Karma | Phase-field model of mode iii dynamic fracture[END_REF], which are not presented here.

Combination between regularized crack and Griffith theory

With the regularized crack field, sharp cracks are approximated by a diffusive band as following:

Γ ≈ Γ l (d ) = B γ (d , ∇d ) dV,
where B ⊂ Ω denotes the localization band over which the crack is smeared. In practice, this band is usually sufficiently small compared to the considered domain Ω. Therefore, the crack surface energy can be approximated as:

Φ d = Γ G c dΓ ≈ B G c γ (d , ∇d ) dV. (4.30)
The strain energy can be expressed as:

Φ s = Ω ϕ (ε (u) , d ) dV = Ω g (d ) ϕ 0 (ε (u)) dV, (4.31)
where g (d ) and ϕ 0 are the energetic degradation function and initial strain energy respectively. In general, the expression of the total energy (Eq. (4.22)) is where σ = ∂ϕ ∂ε denotes the Cauchy stress, and δu = 0 for ∀x ∈ ∂Ω u . In general, Eq. (4.33) can be divided into two parts δΠ (u, d ) P1 and δΠ (u, d ) P2 . The first one is the mechanical part without damage that is written as follows:

Π (u, d ) = Ω ϕ (ε (u) , d ) dV + B G c γ (d , ∇d ) dV - ∂Ω t t * • udS - ∂Ω u (σ • n) • u * dS, ( 4
δΠ (u, d ) P1 = Ω σδεdV- ∂Ω t t * • δudS.

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Applying the divergence theorem, δΠ (u, d ) P1 can be written as:

δΠ (u, d ) P1 = - Ω ∇ • σδudV + ∂Ω t σ • n -t * • δudS,
which is equivalent to the strong form

∇ • σ = 0 in Ω, σ • n = t * on ∂Ω t .
Regarding δΠ (u, d ) P2 , it can be written as:

δΠ (u, d ) P2 = B g ′ (d ) ϕ 0 δd dV + B G c c 0 l c α ′ (d ) δd dV + B 2G c l c c 0 ∇d ∇δd dV. (4.35)
Application the theorem of divergence yields

δΠ (u, d ) P2 = B g ′ (d ) ϕ 0 δd dV + B G c c 0 l c α ′ (d ) δd dV - B ∇ • 2G c l c c 0 ∇d δd dV+ ∂B 2G c l c c 0 ∇d δd • ndS. (4.36)
The surface integral vanishes by considering on the crack boundary the condition: ∇d (x) • n = 0. Thus, we have:

δΠ (u, d ) P2 = B g ′ (d ) ϕ 0 δd dV + G c B 1 c 0 l c α ′ (d ) δd dV - B 2l c c 0 • ∆d δd dV ,
in the case of a homogeneous material, for which G c and l c are constant over the domain Ω.

Due to the irreversibility of Γ, the rate Γl is elaborated as:

Γl = B γdV = B ḋ δ d γdV ≥ 0, (4.37)
by re-writing the term 1

c 0 l c α ′ (d ) -2l c c 0 • ∆d = δ d γ .
The irreversibility in Eq. (4.37) can be satisfied if the following conditions are met:

δ d γ ≥ 0, ḋ ≥ 0. (4.38)
Owing to the unilateral stationary condition of the total energy functional (δΠ (u, d ) P2 ≥ 0), the governing equations of the damage field are

-g ′ (d ) ϕ 0 -G c δ d γ = 0, ḋ > 0, -g ′ (d ) ϕ 0 -G c δ d γ < 0, ḋ = 0. (4.39)
Eqs. 4.38 and 4.39 require a derivable and decreasing degradation function g (d ). Re-writing f ϕ 0 , d = -g ′ (d ) ϕ 0 -G c δ d γ , the phase-field evolution equation can be summarized by the following Karush-Kuhn-Tucker conditions:

ḋ ≥ 0, f ϕ 0 , d ≤ 0, ḋ f ϕ 0 , d ≡ 0,
where f ϕ 0 , d ≤ 0 can be treated as the damage criterion. When ḋ ≥ 0, it follows from

f ϕ 0 , d = 0 that: G c δ d γ = -g ′ (d ) ϕ 0 .
The left part is called "crack resistance" and describes how hard it is to open the crack while the right one is called "crack driving force". The evolution of ϕ 0 can automatically determine the damage variable d , with no need to define pre-cracks.

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

To summarize, in the homogeneous material case, the local damage evolution equation of the phase-field is written as:

g ′ (d ) ϕ 0 + G c c 0 l c α ′ (d ) - 2G c l c c 0 • ∆d = 0, (4.40)
by enforcing ḋ ≥ 0. However, in the heterogeneous material case, G c (x) and l c (x) are functions of x for x ∈ Ω, so the term 2G c l c c 0 ∆d should be corrected as:

2G c l c c 0 ∆d = ∇ • 2G c l c c 0 ∇d , (4.41)
and the damage evolution function should be written as:

g ′ (d ) ϕ 0 + G c c 0 l c α ′ (d ) -∇ • 2G c l c c 0 ∇d = 0. (4.42)
Eqs. (4.40) and (4.42) are the keys for solving the phase-field evolution equation.

Energetic degradation function

Since α (d ), l c , c 0 have been determined in previous sections, the next variable that needs to be determined is g (d ), the energetic degradation function.

In general, this function has to satisfy the following conditions:

• g (0) = 1 (intact state) and g (1) = 0 represents fully broken state,

• g ′ (d ) = dg dd < 0, d ∈ [0, 1), • g ′ (1) = 0.
Based on the papers [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Lorentz | Convergence of a gradient damage model toward a cohesive zone model[END_REF], a general degradation function is introduced, which can be expressed as follows:

g (d ) := (1 -d ) p (1 -d ) p + Q (d ) = 1 1 + ω (d ) , ω (d ) = Q (d ) (1 -d ) p , ( 4.43)
with the exponent p > 0 and the continuous function Q (d ) > 0. Because Q (d ) is strictly positive, the following polynomial is considered:

Q (d ) = a 1 d + a 1 a 2 d 2 + a 1 a 2 a 3 d 3 + ... = a 1 d • P (d ), (4.44a) P (d ) = 1 + a 2 d + a 2 a 3 d 2 + .... (4.44b)
In this type of degradation function, the parameters a i and p need to be determined. In the following, a 1D bar problem is taken as example to illustrate the determination of these parameters.

The bar x ∈ [-L, L] is assumed sufficiently long so that the crack state is not affected by boundary effects. The bar is loaded at both ends by an increasing displacement u * along the x direction and the crack initiates at the symmetric point x = 0. Meanwhile, the crack band is limited in the interval [-D; D], with D < L.

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Because of the constitutive equation, the local strain field can be expressed as:

ε (d ) = σ E 0 g (d ) -1 = σ E 0 [ω (d ) + 1] , (4.45)
for ω (d ) = g (d )

-1 -1 defined in Eq. (4.43). Therefore, the displacement u * imposed at the extremity is given by

u * = σ E 0 L 0 g (d ) -1 dx = σ E 0 L + D 0 ω (d ) dx = σ E 0 L + 1 2 w (σ) ,
where w (σ) := 2σ

σ 2 ω ′ (d ) -A 0 α ′ (d ) -2l 2 c ∆d = 0, A 0 = 2E 0 G c c 0 l c . (4.46)
In the elastic stage, the crack fields are null along the bar, so ∆d = 0. Thus, Eq. ( 4.46) can be written as

σ 2 ω ′ (d ) -A 0 α ′ (d ) = 0, (4.47)
which it is only validated at the point the stress reaches is strength (σ = σ c ) according to Eq. (4.39). Thus, the strength (maximum stress) can be calculated as: . Thus, we have two cases to discuss:

σ c = A 0 α ′ (0) ω ′ (0) . ( 4
• ζ = 0, α ′ (0) = 0, we can get σ c = 0. The damage variable d > 0 since the beginning of the loading and increases gradually with loading increment.

• ζ > 0, α ′ (d ) = ζ > 0, we obtain σ c > 0.
The response of the phase-field has an initial elastic stage in which the material remains undamaged. In this case, the strength σ c can be expressed as: Denoting d * = d (x = 0) when the maximum damage state is reached, and assuming ∇d = 0 at this moment, the stress σ at this point is:

σ c = 2E 0 G c c 0 l c • ζ a 1 . ( 4
σ d * = A 0 α (d * ) ω (d * ) = σ c ω ′ (0) α ′ (0) α (d * ) ω (d * ) . ( 4
σ d * = σ c [ζ + (1 -ζ) d * ] (1 -d * ) p ζP (d * ) . (4.52) Since ∇ • σ = dσ dx = 0, σ (x) = σ (d *
). Inserting Eq. ( 4.51) into Eq. (4.50), we obtain:

d (d ) dx = - 1 l c H d , d * , with H d , d * = α (d ) - α (d * ) ω (d * ) • ω (d ), (4.53) if x is in the interval x ∈ [0, D].
If we inverse and integrate Eq. ( 4.53), we can obtain:

x d , d * = l c d * d H -1 d , d * d d , D d * = l c d * 0 H -1 d , d * d d , (4.54)
where d also describes the damage variable but distinguishes with the value d in the subscript of integration operator. Changing the variable x -→ d (x) and taking Eq. ( 4.54) into account, the apparent jump displacement yields:

w d * = 2σ E 0 D 0 ω (d ) dx = 2 2G c l c c 0 E 0 d * 0 α (d * ) ω (d * ) α d -α (d * ) ω d ω d d d . (4.55)
In the case that 0 < ζ ≤ 2, and by inserting Eqs. (4.29), (4.43), (4.51), Eq. ( 4.55) can be rewritten as follows:

w d * = 4G c ζ c 0 σ c d * 0 P (d * ) (1 -d * ) p • ζ + (1 -ζ) d ζ + (1 -ζ) d * - P d 1 -d p -1/2 d • P d 1 -d p d d . ( 4
w c = 2πG c c 0 σ c ζP (1) lim d * →1 1 -d * 1-(p/2) . (4.57)
According to Eq. (4.49), a 1 will vary and can be expressed as

a 1 = 2E 0 G c σ 2 c • ζ c 0 l c . ( 4

.58)

In the following, we will discuss the coefficients a i , and the degradation function g (d ) for the three abovementioned models.

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

The energetic degradation function of Miehe's model

As presented in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF], the energetic degradation function g (d ) is given as:

g (d ) = (1 -d ) 2 , (4.59)
where a i and p can be obtained from:

a 1 = 2, a 2 = -0.5, a i = 0, i ≥ 3, p = 2. (4.60)
In most implementations, a small positive numerical parameter k is introduced in the degradation function, e.g. g (d ) = (1d ) 2 + k in order to prevent the difficulty of solving the equilibrium equations with zero stiffness at points where d = 1. In general, the value of k should be chosen large enough so that the equilibrium equations are still well-posed but small enough so that the artificial stiffness that is added in totally damaged regions is negligible.

Because ζ = 0 for this method, there is no elastic stage, which means the damage becomes strictly positive right from the beginning of the loading. That means the peak stress can only be reached at the critical damage value (d c ).

∂σ ∂d | d =d c = 0, (4.61)
which can be considered as the strength. Based on [START_REF] Wu | Phase-field modeling of fracture[END_REF], for Miehe's model, d c is given as

d c = 1 4 which conducts to a strength σ c of σ c = 3 16 3E 0 G c l c . (4.62)

The energetic degradation function of Pham's model

The g (d ) of Pham's model as well as the a i values are the same as in Miehe's model. In this method, there is an elastic stage when ζ = 1. Integrating ζ = 1 into Eq. (4.49), the strength in Pham's method can be expressed as:

σ c = 3E 0 G c 8l c . (4.63)

The energetic degradation function of Wu's model

Wu's model incorporates the cohesive zone model, the stress softening law depends on the crack opening displacement. When the crack initiates, the initial slope k 0 is expressed as:

k 0 = lim d * →0 ∂σ ∂w = - c 0 4π • σ 2 c G c • ζ a 2 + p + 1 -1 3/2 ζ 2 (4.64)
with a 2 defined as

a 2 = 1 ζ - 4πζ 2 c 0 • G c σ 2 c • k 0 2/3 + 1 -p + 1 . (4.65)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Accordingly, for a given traction-separation law σ (w) with w c already known, the parameter a 3 in the quadratic polynomial P (d ) can be calibrated as:

a 3 =    0, p > 2, 1 a 2 1 ζ c 0 w c f t 2πG c 2 -(1 + a 2 ) , p = 2. (4.66)
Generally, the values of (a i ) i ≥3 can be omitted.

Linear softening case:

In this thesis, a linear softening law is considered, where the relation between σ and the crack opening displacement w is shown in Fig. 4.9(a). In this case, the initial slope is

k 0 = - σ 2 c
2G c , and if ζ = 2 and p = 2, it is easy to calculate:

a 1 = 4 π • E 0 G c σ 2 c l c , a 2 = - 1 2 , a 3 = 0 (4.67)
where the energetic degradation function is:

g (d ) = (1 -d ) 2 (1 -d ) 2 + a 1 d 1 -1 2 d (4.68)
which is shown in Fig. 4.9(b). Other softening laws can be found in [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF]. From a numerical point of view, the function derivative g ′ (d ) should be monotonically increasing, which equivalently means that the function g (d ) itself should be convex [START_REF] Lorentz | A nonlocal damage model for plain concrete consistent with cohesive fracture[END_REF]. This condition is equivalent to:

a 1 ≥ 3 2 . ( 4

.69)

In practical applications, l c is chosen as small as possible.

Combination between different diffusive crack and degradation functions

Generic degradation functions have been expressed in Eq. (4.43), where a 1 determines the independence of l c , and a 2 depends on the type of softening law.

In the case ζ ̸ = 0, where there is an elastic stage, it may seem a good idea to combine the different types of diffusive crack profile with l c independent degradation function to have different types of l c independent phase-field models. However, the work in [START_REF] Wu | Phase-field modeling of fracture[END_REF] shows that the feasibility of this idea is restricted due to the irreversibility of the damage profile ḋ > 0.

To prove this, we go back to the example of the bar studied in the previous section, where the crack band is limited in the interval [-D, D], with D the semi-crack bandwidth. Based on [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF], the initial semi-crack bandwidth (σ = σ c ) can be expressed as:

D 0 = l c 3 - c 0 π 2 σ 2 c 4k 0 ζ 2 G c (4.70)
where D 0 is a function of ζ. Meanwhile, the ultimate semi-crack bandwidth (d * = 1) can be derived with the help of Eq. (4.54): As for Wu's diffusive crack, the four standard softening laws can all be used. For Miehe's diffusive crack, establishing a l c independent phase-field would require more work.

D u = l c 1 0 H -1 d , 1 d d =      l c 1-ζ ln 2 1-ζ+2-ζ ζ , ζ ∈ (0, 1] l c 1-ζ π 2 -arcsin 2-ζ ζ , ζ ∈ [1, 2]

Choice of stored energy functional

The stored energy functional ϕ describes the energy stored in the body from the intact state to the fully cracked state, which are characterized by the initial free energy density function ϕ 0 and an energetic degradation function g (d ), where ϕ = g (d ) ϕ 0 . For an isotropic elastic body, the initial free energy density can be expressed as:

ϕ 0 (ε) = 1 2 ε : E 0 : ε = 1 2 λ 0 t r 2 (ε) + µ 0 ε : ε = 1 2 σ : S 0 : σ = ϕ 0 σ , (4.73)
where λ 0 and µ 0 represent the lamé coefficients; S 0 := E -1 0 denotes the fourth-order compliance tensor, and σ = E 0 : ε is the effective stress tensor. Besides, under complex loading conditions, it might be relevant to split the free energy density into separate tension and compression part to model the effect of damage only on the tension part. METHOD

Model without split

In the phase-field proposed in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], no split is considered and the free energy density is simply expressed as

ϕ + (d , ε) = ϕ (d , ε) , ϕ -(d , ε) = 0 =⇒ ϕ (d , ε) = g (d ) ϕ 0 (ε) . (4.74)
In this case, the stress field is:

σ = ∂ϕ ∂ε = g (d ) σ, σ = ∂ϕ 0 ∂ε = E 0 : ε, (4.75)
with the evolution of phase-field given by

G c δ d γ = -g ′ (d ) H, H = ∂ϕ ∂g = ϕ 0 (ε) . (4.76)
This type of method ignores the difference between the tensile and compressive behavior and predicts the same damage profile in regions under compression as well as under tension. That is not realistic for brittle or quasi-brittle fracture. Furthermore, it does model the unilateral effect under cyclic loading [START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF][START_REF] Reinhardt | Post-peak cyclic behaviour of concrete in uniaxial tensile and alternating tensile and compressive loading[END_REF], and the crack surfaces may be penetrated by existing cracks under compressive loading.

Miehe's model with split

Miehe [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF] applied a spectral decomposition to the strain tensor in order to fully distinguish the tension and compression strain:

ε = 3 n=1 ε n p n ⊗ p n = ε + + ε - (4.77)
where ε n and p n denote the eigenvalues and eigenvectors of the strain tensor respectively. Note that ε + and ε -represent respectively the positive and negative parts of the strain tensor:

ε + = 3 n=1 〈ε n 〉 + p n ⊗ p n , ε -= 3 n=1 〈ε n 〉 -p n ⊗ p n (4.78)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

with the bracket operators 〈x〉 + = x+|x| 2 , and 〈x〉 -= x-|x| 2 . The operator ⊗ denotes the outer product of two vectors. It then follows that:

ϕ ± 0 = 1 2 λ 0 (〈t r (ε n )〉 ± ) 2 + µ 0 ε ± : ε ± . (4.79)
As damage only affects the tension part, the stress tensor can be expressed as:

σ = g (d ) λ 0 〈t r (ε n )〉 + I + 2µ 0 ε + + λ 0 〈t r (ε n )〉 -I + 2µ 0 ε - (4.80)
with I the identity matrix. The decomposition can be performed by projection tensors P - and P + which are given as follows:

P + := ∂ [ε + (ε)] ∂ε
, and P -= I -P + (4.81)

with I the fourth-order identity tensor. The decomposition of strain can be given as:

ε + = P + : ε, and ε -= P -: ε. (4.82)
The detailed expression of the projection tensor can be found in [START_REF] Miehe | Comparison of two algorithms for the computation of fourth-order isotropic tensor functions[END_REF][START_REF] Miehe | Algorithms for computation of stresses and elasticity moduli in terms of seth-hill's family of generalized strain tensors[END_REF], and Appendix A of [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]. Besides, a shifted split algorithm has been proposed in [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure[END_REF] where the split procedure is performed in the previous loading increment to avoid costly linearization procedures when solving the displacement problem.

Wu's model with split

Similarly to the strain tensor decomposition, the effective stress σ can also be decomposed, as proposed in [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Wu | A novel positive/negative projection in energy norm for the damage modeling of quasi-brittle solids[END_REF]:

σ = 3 n=1 σ n p n ⊗ p n = σ + + σ -. (4.83)
The notations σ n , p n , and σ ± have the same meaning as in the strain decomposition. In general, this approach is simpler than strain decomposition and gives

σ + 1 = 〈σ 1 〉 + , σ + 2 = max σ 2 , ν0 σ 1 + , σ + 3 = max max σ 3 , ν0 σ 1 + σ 2 , ν0 σ 1 + . (4.84)
The orthogonal condition σ + : C 0 : σ -= 0 holds. Parameter setup ν0 = ν 0 1-ν 0 can be used in 3D and plain strain conditions, while ν0 = ν 0 should be used in plane stress condition. Alternative decompositions can be found in [START_REF] Wu | Reconsideration on the elastic damage/degradation theory for the modeling of microcrack closure-reopening (mcr) effects[END_REF]220] With the above split, the positive and negative parts of the energy are

ϕ ± 0 = 1 2 σ ± : S 0 : σ ± . (4.85)
As in the strain decomposition, the stress field is expressed as:

σ = g (d ) σ + + σ -, ∂ϕ ± 0 ∂ε = σ ± . (4.86)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

In the case of isotropic material, to simplify the calculation, a so-called equivalent effective stress σ eq + is defined in [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF] to represent σ + :

σ eq + = 1 1 + β c β c 〈σ 1 〉 + 3J 2 (4.87)
with

β c := σ compression c
σ tension c -1; σ 1 denoting the largest eigenvalue, and J 2 the invariant of the deviatoric part of the effective stress tensor. It then follows that

ϕ + 0 = 1 2 σ eq + 2 E 0 (4.88)
if β c → ∞, σ eq + = 〈σ 1 〉 + , which leads to the Rankine energy:

ϕ + 0 = 1 2 〈σ 1 〉 + 2 E 0 . (4.89)

Hybrid formulation

Different damage behaviors in tension and compression cannot be distinguished in the formulation of no split phase-field. Therefore, the split scheme is often preferred, although it leads to a nonlinear equilibrium equation which requires more computational effort for its solution.

To overcome the above issue, a hybrid formulation is applied in [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. The no split energy and split energy are applied in the stress field and the crack phase evolution, respectively and separately.

Although this hybrid formulation does not violate the second law of thermodynamics, as shown in [START_REF] Wu | Robust numerical implementation of non-standard phase-field damage models for failure in solids[END_REF], it breaks the variational consistency because the stress field and the crack phase-field use different energy functionals. It is still attractive from a computational point of view as the stress is simply given by σ = g (d ) E 0 : ε, which is a linear relation. In addition, Ambati et al. [START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF] confirmed that the crack patterns obtained with this hybrid formulation almost coincide with Miehe's model with single split, with slight differences in the load-displacement curve.

For synthesis, the energy ϕ 0 applied in the mechanical part can be expressed as:

ϕ 0 = 1 2 ε : E 0 : ε (4.90)
while different types of energy ϕ + 0 that can be used for the evolution of the phase-field are given by

ϕ + 0 =        ϕ + 0 (ε) , 1 2
σ eq + 2 E 0 , ϕ + 0 σ . (4.91)
In the unique energy strategy, whether there is a split or none, the displacement and crack fields can be obtained by solving:

(u, d ) = Arg min u,d Π (u, d ) (4.92)
As for the hybrid formulation, the minimization problem in Eq. (4.92) does not hold anymore for both fields (u, d ). A separated minimization function is considered instead:

   u = Arg min u Π (u, d ) d = Arg min d Π + (u, d ) (4.93)
where the energy is not split in the displacement minimization problem but only in the crack problem.

Irreversibility of crack phase-field

Damage models should reflect the physics, where a fully damaged material can not be damaged further, and the damage state cannot be recovered. Consequently, the boundedness d ∈ [0, 1] and irreversibility condition ḋ ≥ 0 should be met and dealt with carefully.

For the quadratic crack function α (d ) = d 2 the boundeness d ∈ [0, 1] is intrinsically satisfied. Therefore, in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], authors only enforced the irreversibility condition when the crack phasefield is close to one:

d (x, t > t 0 ) = 1, if d (x, t 0 ) ≈ 1. (4.94)
In order to prevent cracks from healing when the energy ϕ + 0 decreases, the energy entered at each load increment can be replaced by a local history field, as proposed in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF]:

H (x, t n ) := max 0≤t ≤t n ϕ + 0 (x, t ) . ( 4

.95)

This can make sure the energy entered in each increment is always the largest one among the 0 ≤ t ≤ t n . However, in other cases (ζ ̸ = 0), the bound d ∈ [0, 1] is not always satisfied. In fact, for the phase-field models where the elastic part exists, the damage can be calculated as negative when σ < σ c . This is the reason why extra constraints need to be added.

Eq. (4.39) can be rewritten as follows:

R = -g ′ (d ) ϕ 0 -G c δ d γ ≤ 0 (4.96)
where G c δ d γ denotes the crack resistance and -g ′ (d ) ϕ 0 is the driving force. If R = 0, it means the crack resistance is equal to the driving force where d > 0 and a correct value is calculated. If R < 0, it means the added energy is not sufficient to increase the damage state so d should be stable with respect to the previous increment. Finally, if R > 0, the added energy is larger than the crack resistance so the material cannot sustain the load anymore, therefore d should be set to 1 [START_REF] Amor | Approche variationnelle des lois de griffith et de paris via des modeles nonlocaux d'endommagement: Etude theorique et mise en oeuvre numérique[END_REF][START_REF] Francisco | Finite-Dimensional Variational Inequalities and Complementarity Problems[END_REF][START_REF] Farrell | Linear and nonlinear solvers for variational phase-field models of brittle fracture[END_REF], i.e.:

       0 ≤ d n < d n+1 < 1, R = 0, d n = d n+1 , R < 0, d n+1 = 1, R > 0.
(4.97)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

In [START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF][START_REF] Zhang | 3d micromechanical progressive failure simulation for fiber-reinforced composites[END_REF], the condition (4.97) was replaced by:

H (x, t n ) := max 0≤t ≤t n ϕ + 0 (x, t ) , 1 2
σ 2 c E 0 (4.98)
with the use of the Rankine formulation under the assumption of existence of the elastic stage. This ensures that d ≥ 0 as well as the irreversibility. However, if d > 1, we need force that d = 1.

Synthesis of different phase-fields

It should be remarked that all the descriptions above are based on the assumption that the model is loading-rate independent. Rate dependent models can be found in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF]. In [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], Miehe mentions that applying viscous parameters can stabilize the numerical treatment even for a rate-independent model. However, based on the tests in [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], this stabilization by an artificial viscosity might not be necessary to get converged results. Therefore, in the case of the loading-rate independent model, we do not need to apply viscosity.

Eq. (4.42) is rewritten below:

g ′ (d ) ϕ 0 + G c c 0 l c α ′ (d ) -∇ • 2G c l c c 0 ∇d = 0.
In order to simplify the formulations and help the reader pick up the phase-field key param-

Model α (d ) α ′ (d ) c 0 g (d ) l c Crack bandwidth
Miehe [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF]

d 2 2d 2 (1 -d ) 2 l c = 27 256 E 0 G c σ 2 c ∞ Pham [200] d 1 8/3 (1 -d ) 2 l c = 3 8 E 0 G c σ 2 c 4l c
Wu -Linear [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF] 2d

-d 2 2 -2d π (1-d ) 2 (1-d ) 2 +a 1 d 1-1 2 d Numerical parameter πl c
Table 4.3: The parameters of different type of phase-field methods.

eters, a summary of all the parameters of three phase-fields models is shown in Table 4.3. As for the energy ϕ 0 , it is an independent parameter that can be chosen separately.

Accordingly, the evolution equation of different phase-fields can be written as follows:

Miehe's model -2 (1 -d ) ϕ 0 + G c l c d -∇ • (G c l c ∇d ) = 0, (4.99)
Pham's model

-2 (1 -d ) ϕ 0 + 3 8 G c l c -∇ • 3G c l c 4 ∇d = 0, (4.100)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

Wu's model linear softening law

- 4a 1 (1 -d ) 2a 1 d -4d -a 1 d 2 + 2d 2 + 2 2 ϕ 0 + 2 π G c l c (1 -d ) -∇ • 2G c l c π ∇d = 0. (4.101)

Chapter conclusions

In this chapter, a review of existing damage modeling has been accomplished. First, we have phenomenological failure criteria widely used in industrial engineering, although they are not suitable for complex analysis. Then we have LEFM thanks to the pioneering work of Griffith and Irwin, which is a good tool for studying crack propagation. However, the crack tip stress calculated by LEFM is infinite hence leading to instantaneous crack propagation. In addition, as the displacement field in the standard FEM is continuous, some special techniques have to be implemented to introduce the displacement discontinuity due to the crack. The 3D implementation of these techniques for problems with multi-cracks branching and merging is still a tough task.

The cohesive zone model has been proposed to solve the infinite stress at the crack tip in LEFM. A process zone is introduced around the crack tip, and crack opening is interpreted as atomic or molecular bonds breaking between two surfaces. This approach is suitable for studying the interfacial damage problem, but it only deals with the separation of the crack surface, not its propagation. Extra approaches such as LEFM are necessary to determine the propagation of the crack. Considering the exigence of a pre-defined crack for LEFM, from a fundamental point of view, LEFM and CZM are not self-contained, requiring additional criteria to determine when/where a crack initiates and grows and in which direction.

In parallel with fracture mechanics, the Continuum Damage Mechanics have been developed to deal with the degradation of material properties due to an accumulation and evolution of microcracks. In CDM, the damaged regions can naturally grow, branch, and coalesce without any numerical difficulty. To resolve the mesh dependency, non-local damage models have been proposed, such as the integration and gradient-based models. Then, we have the phase-field method that can be categorized as a gradient-based non-local model. Phasefield uses the damage variable to describe, on the one hand, the degradation of the material, and on the other hand, the cracks, which are regularized as diffusive crack bands. By utilizing the variational approach to brittle fracture, which minimizes the total potential energy, it is possible to solve the displacement fields and cracks simultaneously.

Three different phase-field models, Wu, Pham, and Miehe, have been presented in detail in this chapter. In phase-field methods, the most important parameter is l c , which controls the width of the diffusive crack. For Miehe's and Pham's model, the length parameter l c should be considered as a material parameter [START_REF] Pham | Gradient damage models and their use to approximate brittle fracture[END_REF][START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF][START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF], like in table 4.3. Even though these formulations come from the analytical solution of a 1D bar with homogeneous crack phasefield [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF] or from the experimentation and inverse analysis [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF], they can give good results in some cases, like in [START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF][START_REF] Mesgarnejad | Validation simulations for the variational approach to fracture[END_REF]. However, if l c is considered a material parameter, its value should be fixed for a given material. These estimations are often too large for the problem dimensions, which lead to some bizarre phenomena [START_REF] Mandal | Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture[END_REF]. Furthermore, the diffusive crack bandwidths of Miehe's and Pham's models are quite large, which results in crack merging

Chapter overview

As stated in the previous chapter, the characteristic length l c is a significant parameter in Miehe's phase-field model. However, according to the state of the art [START_REF] Wu | Phase-field modeling of fracture[END_REF][START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF], researchers have often studied its influence on homogeneous materials while the case of heterogeneous materials has rarely been investigated. This is the key part of this chapter. A review of Chen's implementation of Miehe's model will be presented in Section 5.2 as well as the problem statement in Section 5.3. The analysis of l c based on the homogeneous formulation (Eq. (4.40)) will be shown in Section 5.4, and that based on the heterogeneous formulation (Eq. (4.42)) will be presented in Section 5.5.

Review of Chen's implementation of Miehe's model

Inspired from the work of Bourdin [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF], Miehe's phase-field model [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF] is probably the most widely used phase-field model nowadays. However, one of this model's drawbacks is the high computational cost requirement. A sufficiently refined mesh is necessary to accurately describe the gradient term, and the convergence of the conventional alternating minimization algorithm is relatively slow.

Therefore, one of the solutions is using parallel implementation to benefit from the modern computers. Fortunately, as presented in Section 1.5, the FFT method is an intrinsically parallelized solver. Furthermore, some complex mathematical operations in the real space are much simpler in the frequency domain, like the gradient term in FFT is a simple multiplication. Therefore, a combination between the phase-field method and the FFT solver could overcome some of the drawbacks of the phase-field approach. Such a combination and a weak coupling solving algorithm have been developed by Chen in [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], the analysis in this chapter is mainly based on these developments. Because Miehe's method is an l c sensitive approach, where l c should be considered as a material parameter, l c should be unique for each material. In the case of a heterogeneous model, G c and l c are functions of x in the model, so the divergence of G c l c is not equal to zero. On the contrary, in a homogeneous material model, they are constant and can be treated as simple factors. That is the difference between Eqs. (4.42) and (4.40), where these two equations are recalled in the following.

The generic homogeneous formulation (Eq. (4.40)) is

g ′ (d ) ϕ 0 + G c c 0 l c α ′ (d ) - 2G c l c c 0 • ∆d = 0. ( 5.1)
The generic heterogeneous formulation (Eq. (4.42)) is

g ′ (d ) ϕ 0 + G c c 0 l c α ′ (d ) -∇ • 2G c l c c 0 ∇d = 0. (5.2) FIBROUS MATERIAL ON MICRO-LEVEL
However, in most applications in the literature, the heterogeneity in the mathematical equation has been omitted in the case of heterogeneous media, like [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], and like [START_REF] Espadas-Escalante | Mesoscale analysis of the transverse cracking kinetics in woven composite laminates using a phase-field fracture theory[END_REF][START_REF] Sargado | High-accuracy phasefield models for brittle fracture based on a new family of degradation functions[END_REF][START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF][START_REF] Cao | A novel fft-based phase field model for damage and cracking behavior of heterogeneous materials[END_REF] for other phase-fields. Recent studies [START_REF] Ma | FFT-based homogenization of hypoelastic plasticity at finite strains[END_REF][START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF] propose a multi-phase-field method but only for solving anisotropic problems, and [START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF][START_REF] Li | Modeling microfracture evolution in heterogeneous composites: A coupled cohesive phase-field model[END_REF][START_REF] Paggi | Revisiting the problem of a crack impinging on an interface:a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF][START_REF] Guillén-Hernández | A micromechanical analysis of inter-fiber failure in long reinforced composites based on the phase field approach of fracture combined with the cohesive zone model[END_REF] study the interface debonding and the crack propagation without considering this heterogeneity. To summarize, in these papers, Eq. (4.40), which is called homogeneous formulation, is applied instead of Eq. (4.42) even in the presence of a heterogeneous medium.

In this chapter, the homogeneous formulation of Miehe's model that omits the heterogeneity is applied for the first part of the analysis (Section 5.4), and is then compared to the corrected formulation which is introduced in the second part (Section 5.5). Accordingly, at first, Eq. (4.99) can re-written as:

-2 (1 -d ) H + G c l c d -G c l c ∆d = 0 (5.3)
which is the damage evolution equation. Besides, the mechanical local governing equations are

       div (σ (u, d )) = 0, σ = g (d ) λ 0 〈t r (ε n )〉 + I + 2µ 0 ε + + λ 0 〈t r (ε n )〉 -I + 2µ 0 ε -, ε = 1 2 ∇u + ∇ T u .
(

The energy here is always split both for the mechanical and damage part with the strain decomposition method presented in Section 4.3.6. The so-called maximum positive elastic energy history H is defined as

H (ε (x, t )) = max t ∈[0,t ] ϕ + 0 (ε (x, t )) . (5.5)
Additionally, in [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], a viscous regularization is implemented in order to improve the stability of the numerical solution. However, based on the tests of [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], adding the viscous parameter was not very efficient, thus, in the following, this parameter is not taken into account.

In [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], Miehe's phase-field model is solved by a fixed-point algorithm. Eq. ( 5.3) can be simplified into the form:

A t n (x) d (x) -∆d (x) = B t n (x) (5.6)
which can be written as:

1 l 2 c (x) + 2H t n (x) G c (x) l c (x) d (x) -∆d (x) = 2H t n (x) G c (x) l c (x) (5.7) with A t n (x) = 1 l 2 c (x) + 2H t n (x) G c (x) l c (x) , B t n (x) = 2H t n (x) G c (x) l c (x) . ( 5.8)
Because H t n (x) is calculated at each time step t n , A t n and B t n only need to be updated at t n . By introducing a reference term, Eq. (5.6) can be rewritten as:

A t n 0 d (x) -∆d (x) = τ (x) , (5.9)
with τ (x) the polarized term of phase-field: 0 can be given as:

τ (x) = B t n (x) -A t n -A t n 0 , (5.10)
A t n 0 = min A t n + max A t n 2 , (5.11)
which is an homogeneous constant at each time step t n . In Fourier space, the derivation can be easily transformed into a multiplication with frequency vector ξ, where ∆d (x) becomes

F (∆d ) (ξ) = -(ξ • ξ) F (d ) (ξ) (5.12)
with F (-) the Fourier transform. Therefore, a fixed-point algorithm can be easily executed by assuming the polarization term is already known at the i th iteration, which gives:

d i +1 (ξ) = τ i (ξ) A t n 0 + ξ • ξ (5.13)
until the following convergence criterion is verified

e = ||d i +1 -d i || 2 ≤ 1.10 -6 (5.14)
where ||-|| 2 represents the L 2 norm over the model. Fig. 5.1 represents the flow chart implemented in [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] for solving the Miehe's phase-field.

Figure 5.1: The flow chart of Miehe's phase-field method based on the fixed-point algorithm and FFT solver proposed in [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF] (remark: CV is convergence check, ACV is convergence acceleration algorithm).

In this scheme, the damage field is calculated based on the energy calculation at the previous time increment which requires a very fine loading time step to get accurate results. As shown in Fig. 5.2 in [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF], the results are not stable anymore when δε is larger than 5.10 -7 . This time step has been checked in our cases but further verifications are needed for the general case. In summary, for all our analyses done with Chen's implementation of Miehe's model, the time step will be 5.0 . 10 -7 .

Besides, although it has been concluded in Chapter 2 that Willot's rotated scheme is preferable among the improved schemes, we did not succeed in implementing this scheme for the phase-field solver. To avoid using two different schemes for the mechanical and phase-field parts, Moulinec and Suquet's original collocation scheme has been used for both solvers. FIBROUS MATERIAL ON MICRO-LEVEL

Problem statement

Papers [START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF][START_REF] Mandal | Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture[END_REF] show clearly the l c sensitivity in Miehe's phase-field model, which is one of the most critical issues that limit its application. However, this kind of analysis is mainly based on a homogeneous model with only one phase. The impact of l c with a heterogeneous model with several phases like in composites (fiber and matrix system) or polycrystals has rarely been reported, even if the phase-field has been successfully applied on these materials for several years.

Furthermore, as illustrated in the previous chapter, Miehe's phase-field has no initial elastic stage. That means the damage in one material is strictly positive no matter how hard it is to be damaged. Therefore, it is interesting to study the damage field evolution in a heterogeneous material such as composites composed of rigid fibers and a soft matrix.

In the FFT method, every voxel needs to be assigned certain properties, so that even models with pre-cracks like single or double notch models become heterogeneous medium models. Indeed, in those models, artificial elastic properties are usually assigned in the pre-crack or void [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF].

The objective of this chapter is to investigate the influence of Miehe's model's l c for heterogeneous materials. As shown in [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF], this influence can be related to the inter-fiber distance for fiber/matrix materials like composites. A similar but deeper analysis is proposed in this chapter, for the homogeneous formulation in Eq. (4.40) and then the corrected formulation for heterogeneous materials in Eq. (4.42).

In the following, the models are subjected to a macroscopic strain along the x-direction (transverse direction) with a time step δ〈ε xx 〉 = 5 . 10 -7 until the final failure while stress-free conditions are imposed in the other directions (〈σ y y 〉 = 〈σ zz 〉 = 〈σ x y 〉 = 〈σ xz 〉 = 〈σ y z 〉 = 0). These mixed-type loading conditions are applied using the method presented in [START_REF] Kabel | Mixed boundary conditions for FFT-based homogenization at finite strains[END_REF]. It is reminded that periodic boundary conditions are intrinsic to FFT methods. Due to stress-free boundaries, Poisson's effect is active during loading.

Analysis of the impact of l c with homogeneous formulation

First, in this part, the homogeneous formulation presented in Eq. 5.3 is applied for all test cases.

Unit-cell model with two half fibers on the borders (l f constant)

Inspired from the paper [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF], first of all, a unit-cell model with two half fibers on the borders is analyzed (as shown in Fig. 5.3). The model size is 20 µm (0.02mm) with fiber volume fraction V f = 0.55, and the inter-fiber distance (l f ) is 3.29 . 10 -3 mm while one voxel size is used in the thickness. The resolution of the model is N = 225. The blue part is the fiber for which the material is E-glass with properties: E f = 74 000 MPa, ν f = 0.2, and G f c = 9.0 . 10 10 N/mm to make sure that the fiber is unbreakable. The red part is the matrix for which the material is Epoxy MY750 with properties: E m = 4650 MPa, ν m = 0.35, G m c = 9.6 . 10 -4 N/mm which are taken from the paper [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF]. In this series of simulations, the value of l c varies from 2.7 . 10 -4 mm to 8.2 . 10 -3 mm while the inter-fiber distance is fixed. The minimum value of l c satisfies Miehe's criterion which states that l mi n c /h = 3 > 2, where h denotes the voxel size. As shown in Fig. 5.4, the first observation is that the maximum stress decreases with increasing l c . When l c increases from 2.7 . 10 -4 mm to 1.6 . 10 -3 mm, the peak stress almost drops to half its initial value, which confirms the conclusion that l c should be considered as a material parameter.

Furthermore, an interesting fact, which has not been presented in other papers, is that when l f is much larger than l c , it can be seen that the final failure strain will converge to a stable value. However, when the l f /l c ratio is less than 2, the crack bandwidth becomes larger than the inter-fiber distance and the final failure strain increases with l c . This effect cannot be observed for homogeneous models like in Figure 19 of [START_REF] Wu | Phase-field modeling of fracture[END_REF] or Figure 6 of [START_REF] Zhang | Numerical evaluation of the phasefield model for brittle fracture with emphasis on the length scale[END_REF]. It seems FIBROUS MATERIAL ON MICRO-LEVEL that, in our cases, the presence of the fibers delays the damage evolution when the l f /l c ratio is sufficiently small, while the fibers have little effect on matrix damage evolution when this ratio is sufficiently large. In the Table 5.1, the maximum damage values in the fibers at the final failure strain are presented for different l f /l c ratios. It is surprising that the damage in the fiber is not null even if G f c = 9.0 . 10 10 N/mm, which should be high enough to prevent fiber damage. Therefore, a deeper analysis needs to be carried out to investigate the non-physical damage in the fiber. in the fiber at the final failure strain for different l f /l c ratio.

l f /
The local damage fields for different l f l c ratios are presented in Fig. 5.5. When the ratio is equal to 12.2, there is a clear change of the damage field between matrix and fibers. For a ratio of 2.5, the interfaces become smeared and they cannot be distinguished anymore for ratios of 1 and 0.5. Fig. 5.5 implies that the damage in the matrix is diffused into the fiber, and Fig. 5.6 confirms this. As shown in the figure, whatever the l f /l c ratio, the damage profiles are so smooth that we cannot even distinguish the transition between matrix and fibers. Considering that G f c is high enough to prevent fiber failure, the only explanation is that the damage in the fibers is due to the diffusion from the matrix and is not physical.

To confirm this, the value of G f c is changed to 9.6 . 10 4 and 9.0 . 10 20 N/mm. As shown in Fig. 5.7, whatever the value, the macroscopic stress-strain curve does not change, and the damage profile also remains stable. This confirms that the damage in the fiber is artificially caused by diffusion of matrix damage.

Unit-cell model with two fibers (l c constant)

In this section, the model has been changed to the two-fibers unit-cell model presented in Fig. 5.8 with a total fiber volume fraction of 0.15 and one voxel size is used in the thickness.

The resolution of the model is N = 225. Considering that Miehe's model is an l c sensitive approach, herein, we will make the analysis of macroscopic mechanical response and local crack patterns by considering l c as a constant, and varying the inter-fiber instance l f from 2.75 . 10 -3 mm to 2.67 . 10 -4 mm.

The blue part is fiber for which the material is E-glass with properties: E f = 74000 MPa, ν f = 0.2. Meanwhile, the red part is the matrix for which the material is Epoxy MY750 with properties: E m = 4650 MPa, ν m = 0. Furthermore, the l f herein is set as much smaller than twice the distance from the fibermatrix interface to the border, and since periodic boundary conditions are used, the strain concentration point is expected to be the center zone of the model. The closer the two fibers are, the more concentrated should the strain be. Examples of pure elastic tests with different l f values (as shown in Fig. 5.9) can well illustrate this phenomenon. Figs. 5.10 and 5.11 show the local damage fields and the macroscopic stress-strain curves for different l f /l c ratios. Fig. 5.10 demonstrates that when l c is fixed, for both two tests, if l f is too small, the crack will jump to the borders. However, normally, the smaller l f is, the closer two fibers are, and there should be more strain concentration at the center of the model, as shown in Fig. 5.9.

Furthermore, if we look at the stress-strain curves (Fig. 5.11), if we decrease l f from the highest value, there is a first stage where the final failure is progressively brought forward due to an increase of strain concentration. If we continue decreasing l f until a certain value of 9.78 . 10 -4 mm for Test 1 and 6.2 . 10 -4 mm for Test 2, however, there is a second stage where failure is delayed, which is not physical anymore. The only explanation is that l c is too large compared to l f . With the same l c values, the non-physical damage distribution from the matrix to the fiber can be observed for both tests because the damage variation from the matrix to the fiber is continuous. The smaller the l f is, the higher damage of the fibers have. Secondly, the crack bandwidth appears to be slightly affected by the inter-fiber distance. For both tests, if we decrease l f , the crack bandwidth decreases slightly. Besides, the crack bandwidth of Test 2 is smaller than Test 1 because Test 2 has smaller l c .

Unit-cell model with two fibers (l f constant)

Meanwhile, the two-fibers unit-cell model is taken to analyze the behavior in the case of l f being constant but varying the l c values. The model is shown in Fig. 5.13, where the l f is fixed at 2.75 . 10 -3 mm. However, if we look at the local damage fields in Fig. 5.15, on the one hand, for a l f l c ratio from 12 to 3.5, the damage profile becomes more and more diffusive. Note that the minimum value of l c satisfies Miehe's criterion as l mi n c /h = 2.6 > 2, where h denotes the voxel size. On the other hand, if l c is large, as in Fig. 5.15(f ), the crack jumps from the model center to the border, which is obviously non-physical. Indeed, as mentioned before, the inter-fiber distance is much lower than the fiber-border distance., the strain concentration should be at the center, and the damage initiation and the crack propagation should also be at the center as observed in Fig. 5.15(a-e).

In Fig. 5.16, the damage profile along the center cross-section line is plotted. In the figure, the non-physical damage in the fiber can be clearly observed. The case l f l c = 1.7 is shown separately to emphasize that the cracks are on the borders.

Another interesting point is that, in Fig. 5.16, the borders are actually regions of damage concentration for all l f l c ratios, which may explain why the crack will grow on the border if l c is sufficiently large.

Unit-cell model with two fibers (decreasing fiber l c )

As shown before, there is non-physical damage in the fiber due to the matrix damage diffusion. In the phase-field model, l c controls the crack bandwidth, so a smaller l c means thinner damage diffusion bandwidth. Therefore, the idea is to downsize the l c value of fiber to form a thinner damage diffusion bandwidth and reduce the non-physical damage state in the fiber. In this section, the model is taken from Section 5.4.3, where the inter-fiber distance is fixed to 2.75 . 10 -3 mm. The elastic properties and the damage properties of the materials are the same as for Test 1 in section 5.4.2. The characteristic length of the matrix is fixed to l m c = 8.8.10 -4 mm while for the fiber the two values in As shown in Fig. 5.17, the damage in the fibers is reduced from 0.3 to 0.08. In Fig. 5.18, however, we see that the crack in the matrix jumps from the center (l It can be seen that the behavior for Test 1 is normal and physical, as the energy at the center is always higher than at the border until final failure due to the strain concentration. For Test 2, however, the strain energy at all locations perfectly coincides with Test 1 in a first stage but at some point the strain energy at the border grows much faster. It is hard to explain why the crack would initiate at the border from a physical point of view. As a conclusion, the impact of l c for heterogeneous materials can be considered as twofold.

On the one hand, the choice of l m c impacts the behavior of the fiber. On the other hand, the choice of l f c also impacts the matrix behavior. In both cases, this influence due to the formulation is artificial and cannot be justified from a physical point of view.

RVE multi-fiber model

In this part, the RVE shown in Fig. 5.21 is tested for different values of l c . The same value is used for fiber and matrix to avoid the non-physical behavior discussed above. This RVE is reproduced from the paper [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF] as well as the material properties: the blue part is fiber for which the material is E-glass with properties: E f = 74000MPa, ν f = 0.2, and G f c = 9.0 . 10 10 N/mm to make sure that the fiber cannot fail. On the other side, the red part is matrix for which the material is Epoxy MY750 with properties: E m = 4650MPa, ν m = 0.35, G m c = 9.6 . 10 -4 N/mm. In Fig. 5.22, two phenomena can also be observed. The first is that the crack patterns are influenced by l c . The second is that the damage in the matrix is diffused into the fibers. Fig. 5.23 shows also that the peak maximum stress decreases when l c is increased.

Analysis of the impact of l c with heterogeneous formulations

As aforementioned, in heterogeneous media, parameters G c and l c are not constant over the domain. If we use Eq. ( 5.3) in a heterogeneous model, as analyzed in the previous section, the damage in one material diffuses into the others without any physical meaning, and varying the damage properties does not reduce this diffusion.

Accordingly, in this part, the heterogeneity correction (Eq. (4.41)) is added into Eq. ( 5.3) to form a heterogeneous formulation that can be written as:

-2 (1 -d ) H + G c l c d -∇ • (G c l c ∇d ) = 0 (5.16)
without considering the rate-dependent case. In order to implement Eq. 5.16 into the fixed point solver, inspired from the paper [START_REF] Jeulin | Towards crack paths simulations in media with a random fracture energy[END_REF], we introduce:

Q (x) = G c (x) l c (x) .
(5.17)

Introducing a fluctuation term Q ′ (x), we have:

Q (x) = Q 0 + Q ′ (x) (5.18)
where

Q 0 = max (Q (x)) + min (Q (x)) 2 . ( 5.19)
Therefore, Eq. (5.16) can be rewritten as:

G c (x) l c (x) d (x) -∇ • (Q 0 (x) ∇ (d (x))) -∇ • Q ′ (x) ∇d (x) = 2 (1 -d ) H (5.20)
where H, the history, is expressed in Eq. (5.5). Because Q 0 is constant over the domain, Eq. (5.20) gives

G c (x) l c (x) + 2H t n (x) d (x) -Q 0 (x) ∆d (x) = 2H t n (x) + ∇ • Q ′ (x) ∇d (x) .
(5.21)

Setting A t n = G c (x) l c (x) + 2H t n (x), B t n = 2H t n (x), D = ∇ • Q ′ (x) ∇d (x)
, we have:

A t n d (x) -Q 0 (x) ∆d (x) = B t n + D. ( 5

.22)

Polarizing A with A

t n 0 = max(A)+min(A) 2
, Eq. 5.22 gives:

A t n 0 d (x) -Q 0 (x) ∆d (x) = τ (x) (5.23) with τ (x) = B t n -A t n -A t n 0 d (x) + D. (5.24)
At a given time step t n+1 , and iteration i , a new damage field can be solved from Eq. (5.23): While convergence criterion is not met:

d t n +1 i +1 = τ i A t n 0 + Q 0 (ξ • ξ) (5.
(a 0 ) d 0 (x) = 0 (b 0 ) ∇ • Q ′ (x) ∇d 0 (x) = 0 (c 0 ) Calculate Q ′ (x)
(a 2 ) τ i (x) = F[τ i (x)] i t h iteration (b 2 ) d t n+1 i +1 (ξ) = τ i (x) A t n 0 +Q 0 .(ξ•ξ) (c 2 ) d t n+1 i +1 (x) = F -1 [ d t n+1 i +1 ] (x) (c 2 )
Calculate D i +1 (x), and τ i +1 (x)

(d 2 )
Convergence test (5.26) The calculation of D (x) can also be performed in Fourier space, because the derivation operation in real space is a simple multiplication in Fourier space, which can be given as:

D (x) = F -1 J • ξ • F Q ′ (x) • F -1 J • ξ • d (x) (5.27)
where J = -1. Compared to the uncorrected formulation, the correction adds one more Fourier transform, and two more inverse transforms. Besides, the convergence test is the same as Eq. (5.14)

Unit-cell model with two fibers

In this section, the model geometry is the same as in Section 5.4.3, as well as the elastic properties and loading conditions. Regarding the damage properties, the G f c is given as 60 N/mm and G m c = 0.0030 N/mm. Meanwhile, two l c values, summarized in Table 5.4, are taken for this study. Fig. 5.26 presents the macroscopic responses for two l c values after the heterogeneity correction. It implies that this correction has a non-negligible effect on macro-level σε relations. On the contrary to the statement in [START_REF] Chen | Fft phase-field model combined with cohesive composite voxels for fracture of composite materials with interfaces[END_REF], this phenomenon states that Eq. (4.41) can never be omitted for heterogeneous materials. Then damage profiles along the centerline are plotted in Fig. 5.27. Regarding this figure, the damage in fiber is suppressed, which suggests that the non-physical damage diffusion between phases can be stopped by adding this heterogeneity correction. Table 5.5: The total iteration numbers of homogeneous and heterogeneous formula.

In this study, the Anderson acceleration scheme has been applied for both homogeneous and heterogeneous formulation to speed up convergence. This acceleration scheme has been detailed in [START_REF] Chen | Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography[END_REF] and studied for homogeneous formulation in [START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF]. Table 5.5 implies that the heterogeneous formulation will add more than 200 times of iterations despite the accelerating algorithm. Besides, our following study will show that the heterogeneous formula is sensible to the contrast of damage properties, like G c . The higher the contrast is, the more iterations are required. Thus, the application of the heterogeneous formula in high contrast scenarios (e.g., composite with voids) should be further investigated.

Single notch model

In order to better study the impact of the Anderson acceleration scheme and contrast of damage properties, a simpler model (single-notch model) has been studied. As shown in Fig. 5.29, it is a square plate of length 1mm. A straight horizontal notch of 0.5 mm is located at the mid-height of the left edge. The width of the notch is one voxel size. Because the FFT solver is intrinsically periodic, we need to add lateral bands, denoted as M2 (purple material), to avoid spurious effects due to the boundary conditions. Material M2 is also used for the notch. Note that, in this study, the resolution is given as N = 251 (mesh size h = 4.00 . 10 -3 mm) and l c = 0.02 mm. 5.7 shows the total iteration numbers of the fixed-point solver, where several conclusions can be obtained. Compared to the homogeneous formulation, the heterogeneous formula is sensible to the contrast of damage properties. The higher the contrast is, the more iterations are required. On the other hand, the Anderson acceleration scheme can significantly reduce the iteration numbers. They, nevertheless, still have a similar order of magnitude to those without acceleration. Hence, may a more efficient scheme can be investigated in the future.

Chapter conclusions

Thanks to the different tests and analyses conducted in this chapter, several conclusions can be raised.

At first, many works have reported that the l c should be treated as a material parameter because the choice of l c impacts material behaviors. Based on the above studies, we would like to say the second part of this theory «choice of l c impacts material behaviors» has been reconfirmed. However, its first part should be discussed more prudently for heterogeneous materials.

In the case of a heterogeneous material like a composite composed of a soft matrix and rigid inclusions/fibers, the relation between the inter-fiber distance l f and the characteristic length l c has a huge impact on both the crack propagation patterns and the macroscopic stress-strain response. From a mechanical point of view, if two fibers are closer than others, there is more concentration of strain in the matrix zone between these two fibers. Therefore, damage should initiate in this zone. If two fibers are placed closer and closer, damage initiation should occur earlier.

With the application of Miehe's phase-field, this is not verified. On the one hand, if the l f /l c ratio is not large enough (the l f is too small or the l c too large), although there is much more concentration of strain in this inter-fiber zone, the crack jumps to somewhere else. On the other hand, if the l f /l c ratio allows for an inter-fiber crack, but if it is not sufficiently large, CHAPTER 5. SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL damage initiation is delayed and the peak macroscopic stress increases when the inter-fiber distance decreases, which is not physical. Consequently, an appropriate l f /l c ratio is necessary for heterogeneous materials to get a proper crack initiation and propagation pattern and macroscopic stress-strain relation. However, the safe l f /l c ratio varies a lot for different combinations of fiber/matrix properties. Based on our experience, l f /l c ≥ 4 is recommended.

Thirdly, in Miehe's model without heterogeneity correction in the mathematical equation, damage can diffuse from one phase to the other. Based on our analysis of a fibrous material, this can lead to two issues:

• For a composite material under transverse loading, the matrix damage should be dominant. However, the damage in the matrix can be diffused into the fiber, even if the fiber fracture toughness G f c is very large. This diffusion strongly depends on l c . The larger the l c , the stronger the diffusion.

• The interactions of damage behavior between matrix and fiber are twofold. The damage in the matrix can affect the damage in the fiber, and reciprocally. For example, if we only reduce l c in fibers to reduce the fiber damage, at the same time, the damage behavior in the matrix also changes.

Therefore, setting up different l c values between phases by, for example, giving a smaller l c value for a more rigid material to decrease the diffusion, seems difficult. Based on our tests, this idea can help reduce the damage in the fiber. It nevertheless leads to damage initiation and cracks propagation in the matrix in a way that cannot be explained.

Then, a correction of the implementation of Miehe's phase-field model in the FFT solver for heterogeneous materials has been assessed. With this correction, the non-local term becomes ∇ • (G c l c ∇d ), and two changes can be observed:

• The macro-level mechanical responses are changed, which implies that the heterogeneity in the mathematical formulation cannot be omitted.

• The damage diffusion between phases can be stopped at the interfaces.

Hence, the heterogeneous formulation should be applied for heterogeneous materials. Nevertheless, this correction is costly, raising hundred times of complexity based on our studies. On the other hand, this increase of CPU time is sensible to the contrast of material properties. The higher the contrast is, the more iterations are required for convergence. Thus, efficient accelerating algorithms and improved schemes should be investigated in the future.

For a general conclusion, in order to get a regular damage field, the characteristic length l c should be as small as possible. However, in Miehe's method, it is related to material behavior. In our experience, a ratio l f /l c > 4 should be used, which widely limits the application of this phase-field model to heterogeneous materials.

One of the promising solutions for heterogeneous materials may be the implementation of a model which would suppress the sensitivity on l c . As introduced in Chapter 4, one of the candidates is Wu's model [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF]. Thus, Wu's phase-field model will be introduced in the FFT solver in the next chapter, where the heterogeneous formulation is naturally implemented. In this method, l c is considered as a purely numerical parameter that does not affect the mechanical behavior as long as it is chosen small enough.

Chapter overview

Based on the investigations presented in the previous chapter, it is clear that the characteristic length l c has a significant impact on the damage behavior, especially for heterogeneous materials. Crack initiation and propagation are influenced, and non-physical damage diffusion occurs between phases. Some solutions to improve the phase-field model have been attempted. The main idea was to assign l c values as small as possible. In Miehe's model, unfortunately, l c is defined as a material parameter, and a too-small value of l c leads to a nearly linear behavior, which is not physical. We need a phase-field method that would define l c as a numerical parameter. In particular, we would like the results to converge in terms of crack initiation and propagation, dissipated energy and overall strength for a sufficiently small l c . Wu claimed that his method possesses these features in [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF]. We choose to investigate this method in this chapter.

A reminder of Wu's phase-field model is presented in Section 6.2. Sections 6.3 and 6.4 will present some trials to implement it in the FFT solver, and the results of testing cases will be in Section 6.5.

Reminder of Wu's phase-field model

Based on Eqs. (4.56) and (4.57), it can be clearly seen that when there is an elastic stage, the displacement because of the opening of the crack is function of σ c and d . At final failure, when d = 1, the final crack opening displacement w c only depends on σ c and not on l c . Furthermore, the peak stress is limited by σ c . Note that the numerical stability condition a 1 ≥ 3 2 requires l c to be as small as possible. Therefore, any phase-field model that has a purely elastic stage can theoretically be implemented as an l c insensitive method, like [START_REF] Geelen | A phase-field formulation for dynamic cohesive fracture[END_REF]. However, due to the irreversibility of the damage field ḋ ≥ 0, the initial damage localization bandwidth D 0 needs to be inferior or equal to the damage localization bandwidth at final failure D u (D 0 ≤ D u ). Based on Fig. 4.10, in order to integrate the linear softening law, Wu's phase-field is an optimal choice. For Wu's phase-field with the linear softening law, the phase-field evolution equation can be expressed as:

2 π G c l c (1 -d ) -∇ • ( 2G c l c π ∇d ) + g ′ (d )H max = 0, ( 6.1)
where the generic evolution equation of phase-field is written in Eq. (4.42), and the different parameters are written in Table . 4.3. Besides, because g (d ) is a rational function, its derivation is straight forward and we have

g ′ (d ) = - 4a 1 (1 -d ) (2a 1 d -4d -a 1 d 2 + 2d 2 + 2) 2 . (6.2)
The history variable is denoted H and is defined as

H = < σ 1 > 2 2E 0 (6.3)
where the Rankine criterion is used with the stress split, < ->= max{-, 0} being the positive part operator. Stress σ 1 represents the maximum principal stress. Because Wu's phase-field CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER has an elastic stage, if the evolution is used as is, the damage calculated in the elastic stage will be negative. In order to correct this issue, and enforce the irreversibility of the damage field, an extra condition needs to be satisfied:

H max (x, t n ) := max 0≤t ≤t n H(x, t ), 1 2
σ 2 c E 0 . (6.4)
This conditions ensures d = 0 during the elastic stage. If d > 1 in the calculation, then the condition d = 1 is applied.

The difficulty of applying Fixed-point algorithm

Following the procedure used for dealing with Miehe's model with the heterogeneity correction, the evolution equation of Wu's phase-field can be written as:

A 0 d (x) + Q 0 ∆d (x) = τ(x), (6.5) with τ(x) = B -(A -A 0 )d (x) -D. (6.6)
The coefficients are given by 

      A = 2 π G c l c B = g ′ (d )H max + 2 π G c l c D = ∇ Q ′ (x) • ∇d (x) (6.7) with A 0 = max(A)+min(A) 2 , Q 0 = max(Q)+min(Q) 2
, and Q

(x) = Q ′ (x) + Q 0 = 2G c l c π .
The damage field d i +1 at iteration i + 1 at time increment t n can be calculated as.

d i +1 t n = τ i A 0 -Q 0 • (ξ • ξ) . (6.8)
There is an inconsistency with this fixed-point algorithm that can be proven for a simple 1D case discretized with only one element. Indeed, in this case the damage would be distributed in a homogeneous way (∇d = 0), so the non-local term would be equal to zero (D = 0). Thus, it would be possible to rewrite Eq. (6.1) as

2 π G c l c (1 -d ) = 4a 1 (1 -d ) (2a 1 d -4d -a 1 d 2 + 2d 2 + 2) 2 H max (6.9)
and obtain the analytical solution

d hom = 1 - Ma 1 -2Ma 1 H max M(a 1 -2) (6.10)
with M = G c πl c . At the last time increment t n of the elastic stage, where the peak stress would be reached in the element we would have d (x) = 0 and H max = H = 1 2 σ 2 c E 0 . As a result, we would have E 0 . At the first iteration of the fixedpoint algorithm, we would need to calculate τ(x), in which d (x) would be taken from the last iteration of time increment t n , in our case d (x) = 0. We would be facing the problem that H max > 1 2 σ 2 c E 0 would lead to B < 0, which would result in a negative damage field (d i =1 t n+1 ). One might think that getting a negative damage field at some iteration of the fixed-point algorithm would not be a major issue as the algorithm might still converge to a physically acceptable value. Fig. 6.1, however, clearly shows that this would be unlikely. Indeed, this figure shows that once d becomes negative, g ′ (d ) decreases drastically down to very large but negative values, leading to B ≪ 0. This issue is not met for Miehe's model, where the fixed-point algorithm always computes positive damage values. In conclusion, the same algorithm cannot be used for Wu's phase-field model.

g ′ (0) = -a 1 = - 4E 0 G c πσ 2 c l c ⇒ B = g ′ (0)H max + 2 π G c l c = - 4E 0 G c πσ 2 c l c 1 2 σ 2 c E 0 + 2 π G c l c = 0. ( 6

Applying Wu's phase-field model by Newton-Krylov algorithm in FFT solver

Therefore, Newton-Raphson's method, one of the most traditional method for solving nonlinear equations, is applied in this chapter. This approach is based on the iteration:

d i +1 = d i - R(d ) R ′ (d ) , ( 6.12)
where the residual function R(d ) is

R(d ) = g ′ (d )H max + 2 π G c l c (1 -d ) -∇ • ( 2G c l c π ∇d ) (6.
      g ′′ (d ) = a 1 -3(1-d ) 2 +a 1 (1.5d 2 -3d +2) (1-d ) 2 +a 1 d (1-0.5d ) 3 R ′ (d )δd = g ′′ (d )H max δd -2 π G c l c δd -∇ • ( 2G c l c π • ∇δd ). (6.16)
This linear problem is solved using a Conjugate Gradient (CG) solver, which does not require to introduce a reference material or a polarization term. The non-local term is directly computed using

∇ • ( 2G c l c π ∇d ) = F -1 J • k • F Q ′ (x) • F -1 (J • k • d (x)) (6.17)
with J = -1. The FFT and its inverse only need to be computed for the calculation of this non-local term. Note that we did not either be able to implement Willot's rotated scheme for the phase-field solver for this chapter. To avoid using two different schemes for the mechanical and phase-field parts, we used Moulinec and Suquet's original collocation scheme for both solvers. This implementation of the heterogeneous formulation deals correctly with complex numbers both in real and Fourier space in the intermediary steps. This is in contrast with the code used in the previous chapter for Miehe's model.

Regarding the CG solver, we need to introduce LP(x) and RP(x) as follows:

LP(d ) = -2 π G c l c δd + g ′′ (d )H max δd -∇ • ( 2G c l c π • ∇δd ), RP(d ) = -g ′ (d )H max -2 π G c l c (1 -d ) + ∇ • ( 2G c l c π ∇d ). (6.18)
The pseudo-code of the CG solver for Wu's phase-field is presented in (6. [START_REF] Hopkins | Unique set of micromechanics equations for high-temperature metal matrix composites[END_REF]. It is a standard CG solver except for operation f 1 which is a special treatment to remove the accumulation of floating point error. The convergence criterion of the CG solver is set as 1.0 . 10 -6 .

The result of the CG solver is δd , thus, d new = d ol d + δd , and convergence is verified for the new damage field. As opposed to the conventional convergence check in Eq. 5.14, in our implementation, the convergence criterion for Newton-Raphson's algorithm is based on the residual RP(d ) of the phase-field evolution equation. The algorithm is summarized in pseudo-code in (6.20).

In order to simplify the calculation, in this implementation, the hybrid formulation is applied, where there is a split in the history variable used for damage evolution (Eqs. (6.3) and (6.4)) but no tension/compression split in the mechanical part. Thus, the equations for the mechanical part can be written as:

       div(σ(u, d )) = 0, σ = g (d ) ∂ϕ 0 ∂ε , ε = ∇u+∇ T u 2 . (6.21)
This formulation could be solved using the same weak coupling used for Miehe's model (as shown in Fig. 5.1), where we chose to solve the damage and displacement fields only once per time increment. This approach required to use a sufficiently fine time step to obtain converged results. In this chapter, in order to use larger time steps, we preferred a stronger coupling and switched to an alternating minimization algorithm. The flow chart of alternating minimization is given in Fig. 6.2. As shown in this figure, first we have the mechanical solve, and then the phase-field solve, for which we included the details of the Newton-Raphson algorithm. Outside these two loops, there is another loop, called the alternating minimization loop, to make sure that all damage and mechanical fields are well updated in the current time increment. This alternating minimization loop is controlled by the alternating minimization iteration number «al t » and the convergence check «AltCV», which is a based on the damage change:

i = i + 1 (b 1 ) Calculate LP(a, x) (c 1 ) Define q = LP(a, x) (d 1 ) Calculate α = d el new a T .q (e 1 ) Calculate δ i d = δ i -1 d + α * a ( f 1 ) if i is divisible
(b 2 ) d t n +1 i +1 = d t n +1 i + δd (c 2 ) Re-calculate RP(d t n +1 i +1 , x) (d 2 ) if RP < than tolerance (1..10 -6 ), convergence is met (6.20)
eral times until convergence of the alternating minimization scheme, thus ensuring a strong coupling between the displacement field and the phase-field. The advantage of this algorithm is that we do not need a very refined time step to get accurate results. Based on our tests, δ〈ε〉 ≤ 1.00 . 10 -4 is acceptable (results are shown in Section 6.5.1), which is 200 times larger than the value used with the weak coupling.

Numerical examples

The proposed algorithm is designed to solve general 3D problems, but to limit the computation time in our tests we used models with only one voxel thickness in the third dimension. Due to the periodic boundary conditions, the out-of-plane strain components are constant in the third direction, which is equivalent to the generalized plane strain condition.

CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER

As for loading conditions, in these analyses, the notch models are subjected to a macroscopic strain along the y-direction with a loading rate δ〈ε〉 y y = 1.10 -4 (thanks to alternating minimization) until final failure while stress-free conditions are imposed in the other directions (〈σ〉 xx = 〈σ〉 zz = 〈σ〉 x y = 〈σ〉 xz = 〈σ〉 y z = 0). These mixed-type loading conditions are applied using the method presented in [START_REF] Kabel | Mixed boundary conditions for FFT-based homogenization at finite strains[END_REF]. Due to stress-free boundaries, Poisson's effect is active during loading. Regarding the two-fibers composite models, the loading direction is changed to x with δ〈ε〉 xx = 1.10 -4 , and 〈σ〉 y y = 0, while others are unchanged.

Single notch model

Inspired from the literature, the first model we tested is the single notch plate. This example is a popular benchmark test to verify phase-field models for brittle fracture [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Chen | A fft solver for variational phase-field modeling of brittle fracture[END_REF][START_REF] Ambati | A review on phase-field models of brittle fracture and a new fast hybrid formulation[END_REF]. As shown in Fig. 6.3, it is a square plate of length 1mm. A straight horizontal notch of 0.5 mm is located at the mid-height of the left edge. The width of the notch is one voxel size.

Because the FFT solver is intrinsically periodic, we need to add lateral bands, denoted as M2 (purple material) to avoid spurious effects due to the boundary conditions. Material M2 is also used for the notch. In addition, we need to assign properties in the notch as well as the lateral bands that are low enough to model voids but not too low to cause numerical instabilities. Fig. 6.3(a) represents the model for voxel size convergence tests, for which we ensured that the volume fraction of the lateral bands is maintained at 10% for all resolutions. Apart from these voxel size convergence tests, we used bands of a thickness as small as possible to minimize the computational cost. We found that a thickness of 5 voxels, as shown in 6.3(b), is the minimum value to get convergence in the FFT solver and break the periodicity. As shown in Fig. 6.6(a), the macroscopic response changes with the variation of l c , which differs with the conclusion declared in paper [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF]. We observed the same phenomenon using an FEM implementation of Wu's method, as shown in 6.6(b).

To investigate the cause of this sensitivity, we simulated the single notch tests using the FEM solver again but without lateral bands and with a zero-thickness notch created by duplicating nodes (monophase model without void-like material 2). Results with a very fine resolution are plotted as red curves in Fig. 6.6(a) and (b). These curves feature an earlier failure than the first model, which might impute the l c sensitivity to the presence of a void-like material.

From these two figures, it can also be seen that the bigger l c , the higher the peak stress and the later the final failure. The opposite effect is observed with the homogeneous formulation of Miehe's model, as shown in Fig. 6.6(c). The fact that the two models lead to two opposite tendencies suggests that there is a particular effect of the heterogeneity due to material 2 for Wu's phase-field model.

In order to verify this hypothesis, three additional tests are simulated by increasing the damage properties of M2 to G e c = 2.7. The local damage field for resolution N = 251 is shown in Fig. 6.8. It can be clearly seen that the crack bandwidth becomes larger and larger when increasing l c .

In Section 5.5. Besides, the time step convergence studies are also carried out herein, where three time steps are also adopted, which are δ〈ε〉 = 1.00 . 10 -4 , 5.00 . 10 -5 and 2.00 . 10 -5 . We first show the macro-level stress-strain relations in Fig. 6.9. As shown in the figure, the macro-level stressstrain curves are coincident for all three contrasts even if decreasing the time step by five times. Thus, a time step 1.00 . 10 -4 can be acceptable.

Asymmetric double notch model

The next model that we propose to analyze is the asymmetric double notch model, which is presented in Fig. 6.10. It is a square plate of length 1mm with a voxel size thickness. Two notches of 0.2 mm are located at two edges of the model, the one is at 0.55mm of height and the other is at 0.45mm of height. The height of the notch is kept at 1 voxel size as previously and the width of each lateral band at 5 voxel sizes as well. notch model as from the single notch model. As for damage fields, it can be seen that damage diffusion is larger when increasing l c . However, if l c is too large (in our case, l c = 0.03mm, Fig. 6.12(e)), the two cracks connect and merge. Although this merge is not shown in Miehe's model because of the convergence problem, this tendency can also be seen in the local damage fields. The macroscopic responses are shown in Fig. 6.14. The observations and conclusions are similar as previously: when increasing l c , the peak stress becomes larger, and the damage propagation is delayed. Besides, we can also see clearly two stages of crack propagation. The first stage is the crack propagation from the crack tip to the fiber, and the second stage is the propagation along the fiber interface.

Meanwhile, the damage field is presented in Fig. 6.15. The images on the top show the first stage of crack propagation, while the second stage is shown on the bottom. As shown in these figures, for l c = 0.01mm, and 0.02mm, the crack will first propagate from the notch tip to the fiber, then propagate along the fiber interface until final failure. However, for l c = 0.03mm, it seems that crack propagation along the fiber interface occurs independently at the top and bottom, with no connection to the crack propagating from the notch tip in the first stage. This is another proof that the characteristic length l c affects the local damage pattern when it is too large.

CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER

Unit-cell model with two fibers with l c constant

In order to analyze the matrix damage diffusion into the fibers, the two-fibers unit cell model is tested with a fixed l c = 5.0 . 10 -4 mm. All the model details and material parameters are the same as in the model with with l f fixed, except that l f now varies from 2.75 . 10 -3 mm to 2.67 . 10 -4 mm.

The macroscopic response is shown in Fig. 6. [START_REF] Chamis | Mechanics of Composite Materials: Past, Present, and Future[END_REF], where we see a similar tendency as for Miehe's model. In a first stage, when we decrease l f damage initiation and failure occurs earlier. Nevertheless, if we continue decreasing l f , despite the increasing concentration of strain, failure is delayed, and we get a higher peak stress.

This phenomenon can be explained by considering the crack bandwidth as discussed in the previous section. When l f reduces from 1.87 . 10 -3 mm to 1.53 . 10 -3 mm, the l f /l c ratio reduces from 3.74 to 3.06. Since π ≈ 3.14, the fibers end up out of the crack bandwidth, and the response becomes physical (damage initiation and failure occur earlier). However, if we reduce l f from 1.53 . 10 -3 mm to 9.78 . 10 -4 mm, the l f /l c ratio decreases from 3.06 to 1.96, so parts of the fibers end up within the crack bandwidth. The peak stress is higher because of the fibers in the crack propagation path, and the final failure strain becomes larger. Therefore, it is interesting to see if there is damage diffusion in Wu's model. The local damage fields are presented in Fig. 6.19. When l f > 6.2 . 10 -4 mm, the crack passes through the model center, while when l f is too small, the crack propagates on the other side of the fiber-matrix interface rather than through the center. It can also be seen that the damage in the fiber is much lower than in the matrix because the crack at the model center features a necking effect, which means matrix damage diffusion into fibers is successfully prevented.

The damage profiles along the center cross-section line shown in Fig. 6.20 re-confirm these conclusions. For Wu's method, the matrix damage diffusion into the fibers is totally suppressed for all l f values except for the smallest one. On the contrary, we can see a significant and non-physical damage diffusion in Miehe's method.

CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER

Conclusions

In Wu's phase-field model, degradation function g (d ) is nonlinear, which makes it quite difficult to solve Wu's phase-field evolution equation using the FFT solver with a fixed-point algorithm. Thus, we developed a new phase-field solver based on the FFT method. We took advantage of this new implementation to consider a strong coupling of the phase-field with the displacement field to alleviate restrictions on the time step.

This new algorithm involves three nested loops to compute the displacement field and phasefield at each increment. The top loop is the alternating minimization, in which we alternately solve the displacement solver where the degradation function is updated at each iteration and the phase-field solver where the history variable is updated at each iteration. Due to its nonlinearity, the phase-field evolution equation is solved using a Newton-Raphson loop, of which each iteration involves a linear system of equations that is solved using a CG solver. The last loop is the CG solver itself, where at each iteration, the FFT and its inverse are used to compute the non-local term with the correct formulation for heterogeneous materials.

Wu's model is l c insensitive in the literature, but only for homogeneous models [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF]. Based on our tests, for a heterogeneous model, if the crack initiates at the interface between different materials, l c impacts the response. The larger l c , the bigger the impact is. On the contrary, in our models inspired from the composites, if the crack diffusion bandwidth does not include a second material, the method becomes l c insensitive. In Wu's phase-field, the crack bandwidth is well defined and has a size of πl c . In the case of composite materials, if l f ≥ πl c , the model response is surely l c insensitive. Furthermore, the crack patterns can also be impacted if l f /l c ratios are too small even with Wu's phase-field model. Fortunately, in Wu's model, l c is considered as a purely numerical parameter and can be chosen small enough to get converged predictions, as long as it remains large enough compared to the voxel size.

Chapter 7

General conclusions and perspectives

Conclusions

In this Ph.D. thesis report, three years of research on elastic and damage modeling of heterogeneous materials based on the Fast Fourier Transform (FFT) and phase-field models have been summarized. The FFT and phase-field methods have been applied diversely in different domains in the literature. This work has been the opportunity to review a large part of the existing literature. In the beginning, our objective was to use FFT and phase-field to study the textile composite, where the voxel-based mesh is preferable. However, when profoundly studying these methods from the most uncomplicated cases, we found three major issues: The first is the presence of strong oscillations in FFT methods. The second is the non-physical damage diffusion between phases in phase-field methods, and the third is the sensitivity of l c of Miehe's phase-field model on crack patterns. Thus, we decided first to resolve these issues, which led to our several original contributions. These reviews and contributions have been presented in six chapters, with the scope of the three first chapters restricted to elastic behavior and that of the three last chapters widened to model damage and fracture.

In the first chapter, textile composites have been introduced, and a review of existing modeling techniques for these materials has been presented. Different analytical approaches for the microscopic scale as well as the mesoscopic scale have been reviewed and compared to evaluate their differences, advantages, and drawbacks. Due to the demand for complex structure and behavior modeling, it has been shown that numerical approaches have become more and more popular, and in particular the Finite Element Method (FEM).

In the same chapter, we have also discussed different geometry modeling techniques for the microscopic and mesoscopic scales to form realistic microstructures and mesostructures of composite materials. Due to the complexity of these geometries, automatically generating conformal meshes is a tough task, with multiple issues such as non-physical interpenetration and voids. Although correction techniques exist to deal with those issues, our review has led us to conclude that they can involve high computational costs.

Consequently, we have turned our investigation to voxel-based meshes, which offered the possibility to overcome some drawbacks of conformal meshes, and also create models directly from images. In this chapter, we have discussed the steps for generating the voxel From the fourth chapter, we have turned our attention to damage modeling in composite materials. In this chapter, a review of existing damage modeling has been presented. Among damage modeling techniques, the phase-field model has gained an increased attention because it allows dropping the requirement of re-meshing or front tracking methods to follow the evolution of the discontinuities. Indeed, the evolution of the phase-field variable (d ) can completely describe the crack and the damaged regions of the material by regularizing the sharp cracks as a diffusive crack band. Meanwhile, the displacement field and crack propagation are computed simultaneously by minimizing the total energy with a variational approach to brittle fracture. Therefore, in the phase-field model, we found that crack initiation, growth, and coalescence can be automatically determined, and multi-cracks merging and branching problems can be solved without additional efforts.

In this chapter, we emphasized that the phase-field method imposed severe restrictions regarding the mesh size, which could only be satisfied either by an adaptive meshing technique or a uniformly refined mesh. This is why we came up with the idea of combining the phasefield and FFT methods. Our vision was that the efficient and simple parallelization of the FFT and its inverse for voxel-based meshes made the FFT method a serious candidate for phase-field modeling of damage and fracture in composite materials.

We focused on three phase-field models, namely Wu's, Pham's, and Miehe's, which we presented in details in this chapter. In phase-field methods, l c is an essential parameter, which controls the width of the diffusive crack. In Miehe's and Pham's models, the length parameter l c appeared in our review to be considered as a material parameter, while it appeared as a purely numerical one in Wu's phase-field. Most studies, however, discussed the influence of l c only for homogeneous materials.

In the fifth chapter, we analyzed the impact of l c in heterogeneous media with Miehe's phasefield model. Our main concern was that in Miehe's model, l c should be fixed, as it is considered a material parameter, but in heterogeneous materials like composites, the fibers can be arranged randomly, and we can find a wide range of inter-fiber distances l f . Therefore, we investigated the case of a composite composed of a soft matrix and unbreakable rigid inclusions/fibers. We identified that when l f was small compared to l c (typically l f /l c < 4), non-physical damage diffusion occurred from the matrix to the fibers, which corrupted the results both in terms of crack patterns and overall mechanical response.

In summary, on the one side, we need an appropriate l c value to obtain a reasonable macroscopic mechanical response. On the other side, a strong limitation is set to this choice because of the local crack patterns. Therefore, the application of Miehe's phase-field model for modeling heterogeneous materials is strongly limited.

Finally, we found during our analysis that our implementation of Miehe's phase-field model in the FFT solver was not correct for heterogeneous materials. We presented a correction in this chapter as well as its implementation. Although we could not successfully implement this correction in the FFT solver used in this chapter, our preliminary results showed that this correction could reduce the non-physical damage diffusion from matrix to fibers but led to non-physical crack patterns.

In the last chapter, we chose to study Wu's length insensitive phase-field model. Our idea was that since in composite materials, the characteristic length l c should be small enough compared to the inter-fiber distance l f , we needed a phase-field model that defined l c as a purely numerical parameter which could take arbitrarily small values.

In this chapter, we first presented the implementation of Wu's model in the FFT solver by applying the Newton-Krylov method. In addition, we used an alternating minimization to alleviate restrictions on the time step. We also took advantage of our new implementation to include the correction of the formulation for heterogeneous materials.

Although Wu's model has been proved to be l c insensitive for homogeneous materials in the literature, our tests on a heterogeneous material led us to differ. When the crack initiated near an interface between phases with different damage properties (e.g. G c , and σ c ), we found that l c had an impact on the behavior. Only when the second material was not inside the crack diffusion bandwidth did the method become l c insensitive.

In the case of composite materials, small l f /l c values also led to a non-physical macroscopic response and crack pattern as observed with Miehe's model in the fifth chapter. Fortunately, in Wu's model, l c is considered as a purely numerical parameter. Accordingly, we could choose a characteristic length l c as small as possible as long as we also respected Miehe's criterion regarding the mesh size (l c /h > 2, h being the mesh size). For a sufficiently small value, we observed that the non-physical damage diffusion from matrix to fibers was reduced and that the predictions became length-insensitive both in terms of crack pattern and overall response.

To summarize, the key original contributions of this Ph.D. thesis are as follows:

• Developed the method to reduce the numerical artifacts in the FFT solver.

-We have quantitatively compared and ranked different causes of oscillations in FFT methods when modeling the elastic response of heterogeneous materials. It shows that the non-smooth interface due to the voxel-type discretization is the most critical cause.

-Numerical experiments show that Willot's rotated scheme is a good compromise between the oscillation reduction effect and CPU time. However, it cannot suppress the oscillations due to the non-smooth interface, requiring an interface smoothing technique.

-The neighbor voxel average method has been re-proposed in the FFT solver. Meanwhile, an enhanced composite voxel method has been proposed to facilitate its application on image-based models, where the non-parametric representation of interfaces is often the case. The latter has a better oscillation reduction effect, which is more preferred.

• Developed the methods to improve the FFT solver's phase-field accuracy and generality.

-On the contrary to the literature, we have proven that the homogeneous formulation cannot be omitted, which leads to damage diffusion for heterogeneous materials. Thus, we implemented this heterogeneity correction for Miehe's and Wu's model in FFT.

-A thorough analysis of Miehe's phase-field has shown that choosing an appropriate l c is a challenging task, Wu's phase-field model was then developed in FFT solver.

-Numerical results show that Wu's model is partially l c insensitive for heterogeneous materials. If one phase is outside another phase's fracture zone, the model is l c insensitive. Otherwise, the model becomes l c sensitive. To reduce or avoid this effect, l c in Wu's model should be chosen as small as possible.

Besides, these contributions have been disseminated in the form of:

• a paper in a peer-reviewed international journal:

-X.Ma,M. Shakoor, D. Vasiukov, S. V. Lomov, and C. H. Park, "Numerical artifacts of fast fourier transform solvers for elastic problems of multi-phase materials: their causes and reduction methods," Computational Mechanics 2021 67:6, pp. 1661-1683, Apr. 2021,

• two oral presentations at international conferences (the author of this thesis being the presenter for the first one):

-X. Ma, M. Shakoor , D. Vasiukov, S.V. Lomov, and C.H. Park , "Numerical artifacts in FFT solver for elastic problems of multi-phase materials and their reduction with a level-set based composite voxel method," in ACEX conference, Malta, July 4-8, 2021 -D. Vasiukov, S.V. Lomov, X. Ma and M. Mehdikhani, "In situ optical observation of the transverse crack nucleation and growth in cross-ply laminates", in HyFiSyn Conference, Leuven(Belgium), Sept. 15-16, 2021

• and a poster presentation at an international work-shop:

-X. Ma, K-K. Parvathaneni, S.V. Lomov, D.Vasiukov, M.Shakoor and C-H.Park, "Quantitative comparison between fast fourier transform and finite element method for micromechanical modeling of composite," in FiBre-MOD conference, Leuven(Belgium), Dec. 9-12, 2019

Perspectives

In addition to original contributions, we also identified several new possibilities that this research work has opened and should be considered in future work. Several flaws of the current work that we noted should also be filled.

Future works on our enhanced composite voxel method

In our actual enhanced composite method, we need to calculate the signed distance function from each voxel center to the interface. However, the calculation is based on the non-smooth interface that generates a non-smooth signed distance function. We plot the signed distance function of the first layer voxel of the matrix along the interface (model C in Chapters 2 and 3 is used) in Fig. 7.1(a) that clearly presents this phenomenon. A possible idea is to apply a low-pass filtering process to reduce the oscillation state of the signed distance function, like in Fig. 7.1(b). May this improve the oscillation reduction of the present enhanced composite voxel method again. where Π ph i describes the strain energy for each phase. This means the energy stored in sharp and diffusive interfaces should be the same.

Acceleration schemes for phase-field models with heterogeneous formulation

Chapters 5 and 6 show that the heterogeneous formulation, also named complete formulation, is very important for heterogeneous materials. However, we found their limits in terms of iteration numbers. In particular, it is sensitive to the contrast of damage properties, which will cause issues in the presence of voids in composite materials. This increase of computational complexity occurs for both two phase-fields.

In Miehe's model, the Anderson acceleration scheme has been applied. It can reduce the iteration numbers but with limited effect. As for Wu's model, none of the accelerating algorithms has been added. Thus, a more efficient scheme can be proposed for these two phase-fields in the future.

Multi-phase-fields

As discussed in the sixth chapter, Wu's phase-field model becomes l c sensitive when the crack initiates at the interface. One idea proposed in this thesis is choosing an l c as small as possible. Another idea is to use multi-phase-fields [START_REF] Ernesti | Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures[END_REF], where each phase has its own damage variable (d ph i ) that is solved independently. One of the drawbacks is that calculations become much more complex for the multi-phase system. On the other hand, when using multi-phase-fields, each phase-field only deals with one phase, where the homogeneous formulation can be applied naturally. Thus, it is interesting to compare the CPU time of these two methods (multi-phase-fields vs. heterogeneous formula) for multi-phase models.

interface and interphase modeling

On the one hand, in this work, we have focused on damage and fracture modeling mainly in the matrix phase, although we could set up different properties for each phase, including fibers. Interface modeling, however, was completely discarded. Several works have been carried out in the literature recently to access this point, where there are two main directions.

One is coupling the phase-field model with the cohesive zone model like in [START_REF] Paggi | Revisiting the problem of a crack impinging on an interface:a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model[END_REF][START_REF] Chen | Fft phase-field model combined with cohesive composite voxels for fracture of composite materials with interfaces[END_REF]. Another is to smear the sharp interface by a diffusive interface band with another phase-field like in [START_REF] Zhang | Modelling progressive failure in multi-phase materials using a phase field method[END_REF][START_REF] Li | Modeling microfracture evolution in heterogeneous composites: A coupled cohesive phase-field model[END_REF]. Inspired by the latter, an interesting and innovative point is to replace the diffusive interface (phase-field) with our enhanced composite voxel method, which creates a diffusive interface band between two phases.

Thus, we actually have three phases for a two-phase composite system: matrix, fiber, and interface formed of the composite voxels. The interface herein has a thickness hence can also be considered as interphase. Interface or interphase? Two directions can be followed: One introduces energy conservation regularization functions to recover the sharp interface as aforementioned. Another is considering this "interface" as true interphase because [START_REF] Wang | Effect of heterogeneous interphase on the mechanical properties of unidirectional fiber composites studied by fft-based method[END_REF][START_REF] Wu | Imaging the interphase of carbon fiber composites using transmission electron microscopy: Preparations by focused ion beam, ion beam etching, and ultramicrotomy[END_REF] have implied that a transition domain between matrix and fiber for composite materials truly exists.

As for the damage model, Wu's phase-field model can be adopted because l c can be chosen as small as possible to reduce the l c effects.

More complex models

As mentioned in Chapters 5 and 6, we did not succeed in implementing Willot's rotated scheme for phase-field solvers. Filling this gap is also an important part of future work.

After finishing all the above works, our parallelized and efficient implementation enables the modeling of complex structures with a reasonable computational cost. More complex models can be achieved. The first part can be the study of the 3D crack propagation for UD composites, where these experiments are accomplished with 3P bending tests, shown in Fig. 7.3. In the figure, first, we can see a UD composite material under 3P bending loading. Then, a high-speed camera is set to capture the crack propagation along the transverse direction (Fig. 7.3(b)). Meanwhile, another camera is used for longitudinal direction crack propagation (Fig. 7.3(a)). With this setup, we can observe the 3D crack propagation frame by frame. Our objective is to create a real microstructure model and observe the real-time crack propagation in virtual experiments using FFT and phase-fields. The second part can go back to the original background of this thesis: Multi-scale modeling of damage behavior of textile composites with detailed microstructures. Even though these studies have been massively carried out for years, they were often limited in one scale. An interesting idea can be found in [START_REF] Fang | A coupled FE-FFT multiscale method for progressive damage analysis of 3D braided composite beam under bending load[END_REF], in which they coupled the macro-level and meso-level, considering yarns as homogeneous materials.

Therefore, further work can couple the meso-level and micro-level that involve the fibers and the matrix on the microscopic scale in the yarns. Because the FFT methods can enable the creation of millions of voxels, one of the scenarios can be replacing the homogeneous yarns with the real microstructure of fibers and matrix, which will create a huge model with high accuracy. If this idea is far beyond our capability, a second option is to re-utilize the logic in [START_REF] Fang | A coupled FE-FFT multiscale method for progressive damage analysis of 3D braided composite beam under bending load[END_REF], which associates each yarn voxel to a microscopic RVE, whose behavior can be obtained by an upscaling procedure from the results of the microscopic scale.

Besides, considering that some types of polymers in composite materials often show loading rate-dependent nonlinear plastic behavior, further investigation can be carried out by taking it into account. For modeling loading rate-dependent behavior, a viscous parameter can be simply introduced like in [START_REF] Miehe | Thermodynamically consistent phasefield models of fracture: Variational principles and multi-field fe implementations[END_REF][START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF]. The papers [START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF] and [START_REF] Gélébart | Non-linear extension of FFT-based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials[END_REF] can also help us for achieving this goal, as elasto-viscoplastic behavior of polycrystals has been modeled using an FFT solver in the first one, and an accelerated version of the FFT solver for non-linear cases has been proposed in the second one.

List of acronyms

Abstract

This PhD thesis addresses numerical modeling of the fracture of heterogeneous materials using a Fast Fourier Transform (FFT) based method and a phase-field model. FFT-based numerical methods relying on the voxel type mesh, show higher computational efficiency and similar accuracy as the finite element method for the same mesh type. These methods, however, are well-known to lead to numerical artifacts (spurious oscillations). The first part focuses on these artifacts and their causes. A neighbor voxel average and an improved composite voxel method relying on a signed distance function are proposed to reduce those oscillations. The second part of the thesis focuses on damage modeling using a length-insensitive phase-field model. The correct implementation of this model for heterogeneous materials within the FFT solver is presented, and it is shown that this model suppresses the influence of the characteristic length on local and global responses as compared to classical models.

Résumé

Cette thèse de doctorat aborde la modélisation de la rupture de matériaux hétérogènes en utilisant une méthode numérique basée sur la Transformée de Fourier Rapide (TFR) et une approche champ de phase. Les méthodes basées sur la TFR ont démontré une haute efficacité et une précision similaire à celle des éléments finis pour des maillages de type voxel. Cependant, elles génèrent des artefacts (oscillations) numériques. La première partie de cette thèse se concentre sur l'étude de ces artefacts et de leurs causes, et sur les traitements numériques permettant de les atténuer. La deuxième partie de cette thèse concerne la modélisation de la rupture à l'aide d'une approche champ de phase insensible à la longueur caractéristique. La mise en oeuvre correcte de ce modèle pour des matériaux hétérogènes dans le solveur TFR est présentée, et il est montré que l'influence de la longueur caractéristique est supprimée par rapport aux modèles classiques.

xxv

Contents 1 . 1 1 . 2 1 . 4 1 . 5 Fast

 11121415 Chapter overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Material introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Composite materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Hierarchy of textile composite material . . . . . . . . . . . . . . . . . . . 1.3 Analytical modeling techniques for composites . . . . . . . . . . . . . . . . . 1.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Micromechanics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Mesoscopic scale homogenization techniques . . . . . . . . . . . . . . . 1.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical modeling techniques for composites . . . . . . . . . . . . . . . . 1.4.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Microscale geometry modeling . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Mesoscale geometry modeling . . . . . . . . . . . . . . . . . . . . . . . . 1.4.4 Meshing technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fourier Transform method (FFT) . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Mathematical description of basic scheme . . . . . . . . . . . . . . . . . 1.5.3 Advantages and drawbacks of FFT . . . . . . . . . . . . . . . . . . . . . . 1.5.4 Improvements of FFT methods . . . . . . . . . . . . . . . . . . . . . . . 1.5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . CHAPTER 1. INTRODUCTION OF MATERIALS AND MODELING TECHNIQUES

Figure 1 . 1 :

 11 Figure 1.1: Schematic of three models proposed by Ishikawa and Chou [36]: (a) mosaic model, (b) undulation model, and (c) bridging model.

CHAPTER 1 .Figure 1 . 2 :

 112 Figure 1.2: The principle of the 3SHM model [44].

Figure 1 . 5 :

 15 Figure 1.5: Multi-scale modeling strategy for textile composites.

Figure 1 . 6 :

 16 Figure 1.6: Several structured periodic fiber arrangements of fiber volume fraction V f = 0.475 [55].

CHAPTER 1 .

 1 Figure 1.7: Modeling of a binder yarn by the DEM [78].

CHAPTER 1 .

 1 Figure 1.8: The illustration flow chart of realistic modeling process based on a microscopic image [80].

CHAPTER 1 .Figure 1 .

 11 Figure 1.10: General algorithm to generate a conformal mesh along the contact zone between yarns [89].

CHAPTER 1 .Figure 1 . 11 :

 1111 Figure 1.11: Assignment of voxels to yarns for different voxel sizes [93].

CHAPTER 1 .Figure 1 . 13 :

 1113 Figure 1.13: Scale separation for composite materials.

CHAPTER 1 .Figure 1 . 14 :

 1114 Figure 1.14: Comparisons between FFT and FEM on (a) macro-level σϵ relation and (b) CPU time (8 cores parallelization) (reproduced from [121]).

Figure 1 . 15 :

 115 Figure 1.15: Some drawbacks of the FFT basic scheme algorithm (reproduced from Lionel Gelebart -FFT-based solvers to evaluate stress distributions in RPV steels -SOTERIA Training School 2018).

Figure 1 .

 1 Figure 1.16: Number of iterations at convergence as a function of the (a) elastic contrast (Young moduli ratio) and as a function of the (b) reference material for an elastic contrast of 10 3 (α is multiplier of basic scheme's optimum reference [126]).

Figure 2 . 2 :

 22 Figure 2.2: Illustration of hourglass effect: (a) 2D finite element with reduced integration under a volume-conserving shear strain, (b) Deformed element with no strain energy.

CHAPTER 2 .

 2 Figure 2.3: Comparison of local stress fields between voxel mesh and conformal (consistent) mesh of textile composite for FEM [93].

Figure 2 . 4 :

 24 Figure 2.4: Illustration of Willot's rotated scheme: A pixel with edges parallel to the original Cartesian basis (e 1 , e 2 ) and the rotated basis with 45 • ( f 1 , f 2 ). The strain and stress fields are evaluated at the pixel center x. The displacement and the divergence of the stress field are output along the pixel corners.

Figure 2 . 5 :

 25 Figure 2.5: Comparison of local stress fields between different Green operators: (a) Green operator of Moulinec-Suquet ( Γ M-S ), (b) Green operator of Willot with centered scheme ( Γ WI-C ), (c) Green operator of Willot with forward-and-backward scheme ( Γ WI-W ), (d) Green operator of Willot with rotated scheme ( Γ WI-R ) [120] (remark: (b), (c) and (d) are the schemes based on finite differences).

  (a) Nodes numbering of the trilinear hexahedral reference element. (b) The displacements are at the element corners, and the strains are at the Gauss points.

Figure 2 . 6 :

 26 Figure 2.6: Nodal numbering and Gauss points for trilinear hexahedral finite elements.

Figure 2 .

 2 Figure 2.7: Illustration of spurious oscillations (check-board patterns) of local stress field (σ xx ) present in FFT with rotated scheme proposed in [120] (left half part) and FEM (right half part) under tension in x-direction with the contrast between the fiber and matrix Young's moduli set to E f /E m = 22.

CHAPTER 2 .Figure 2 . 8 :

 228 Figure 2.8: Illustration of the three models: (a) model A, (b) model B and (c) model C (remark: the size along y-direction of model A is set as one voxel size).

CHAPTER 2 .Figure 2 . 9 :

 229 Figure 2.9: The variation of elastic properties (E x and G x y ) of different FFT algorithms with different resolutions: (a) and (b) are for model A, (c) and (d) are for model B (remark: in all four figures, the curves of W-S and FEMR are superimposed, and the curves colors are the same for (a)-(d)).

CHAPTER 2 .

 2 Figure 2.10: The variation of fiber volume fraction of model C with different resolutions with a target V C f = 0.55.

  225, 315, and 405. The resolutions chosen for model C are N C =225, 255, 305, 345, 375 and 455. The fiber volume fractions in this part are kept as V A f = 0.6 for model A, V B f = 0.36 for model B and V C f = 0.55 for model C.

CHAPTER 2 .Figure 2 . 11 :

 2211 Figure 2.11: The interface normal and tangential strain (ε r r and ε θθ ) of model C, as well as the angle θ.

Figure 2 .

 2 Figure 2.12: (a) The average relative errors for FFT results compared with FEMR and FEMF of the model A of different resolutions, (b) The relative errors along the model A for FFT results compared with FEMR and FEMF of resolutions 315.

CHAPTER 2 .

 2 INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT Regarding model B, Figs. 2.13(a,b) show the maximum error along the interface and Figs. 2.13(c,d) shows the error on the maximum strain along the interface. To assess the hourglass effect, FEMR is used as reference for Figs. 2.13(a,c) while FEMF is used for Figs. 2.13(b,d).

Figure 2 . 13 :

 213 Figure 2.13: Model B: maximum relative absolute error along the interface (a) with FEMR and (b) with FEMF, relative absolute error on the maximum strain along the interface (c) with FEMR and (d) with FEMF (remark: the curves colors are the same for all plots).

Fig. 2 .

 2 Fig.2.14 shows the average interface (orange contour) relative errors (a,b), the average relative errors at the fiber center (green point) (c,d) and those of whole volume (e,f ). Similarly as in Fig.2.13, the left parts of Fig.2.14 are the relative errors with FEMR and the right parts are those with FEMF. These two figures re-confirm that W-S is always more accurate when compared with FEMR, while SF is more accurate when compared with FEMF, except for the average relative errors at the fiber center.

CHAPTER 2 .Figure 2 . 14 :

 2214 Figure 2.14: Model B: The average errors for FFT results along the interface compared to: (a) FEMR and (b) FEMF. The average errors for FFT results at the fiber center compared to: (c) FEMR and (d) FEMF; The average errors of FFT results over the whole volume compared to: (e) FEMR and (f) FEMF (remark: the curves colors are the same for all plots).

Figure 2 . 15 :Figure 2 .

 2152 Figure 2.15: Model C: maximum relative absolute error along the interface (a) with FEMR and (b) with FEMF, relative absolute error on the maximum strain along the interface (c) with FEMR and (d) with FEMF (remark: the curves colors are the same for both all plots).

  Figure 2.17: Illustration scheme of strain variation along a line in a direction (reference curve VS oscillated curve).

CHAPTER 2 .Figure 2 . 18 :

 2218 Figure 2.18: Model B: The D values of horizontal interface line with resolution 135 (a) and resolution 405 (b). The D values of vertical interface line with resolution 135 (c) and resolution 405 (d).

Figure 2 . 19 :

 219 Figure 2.19: Illustrations of different MD values : (a) low MD values (weak oscillated curve), (b) high MD values (strong oscillated curve).

CHAPTER 2 .

 2 INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

Figs. 2 .

 2 Figs. 2.20, 2.22 and 2.25 present the average of MD j over interface lines (denoted MD) for model A, B and C respectively, which can be used to evaluate the impact of each cause of oscillations. In Figs. 2.20 and 2.22, the local strain component ε xx is chosen to study the oscillations while the radial strain component ε r r and the tangential strain component ε θθ are chosen instead in Fig. 2.25.

Figure 2 . 20 :

 220 Figure 2.20: The value of MD for model A using the different methods.

  2.22(a) shows the value of MD along the horizontal interface line. Fig. 2.22(b) shows the MD value along the vertical interface line. Note that the strain values for calculating MD are taken on the matrix side, which are more critical during loading (orange line in Fig. 2.22).

Figure 2 .

 2 Figure 2.21: (a) describes the ε xx component of model B along the horizontal interface line with N B = 405; (b) shows the ε r r component of model C along the interface with N C = 375 (remark: the green line in (b) is the curve of FEM of full integration with conformal mesh, and the curves colors are the same for both (a) and (b)).

Figure 2 . 22 :

 222 Figure 2.22: The values of MD for model B along (a) the horizontal interface line and (b) those of the vertical interface line. (remark: the curves colors are the same for both (a) and (b)).

Figure 2 .

 2 Figure 2.23: Model B: (a) the values of function D (see Eq. (2.9)) of the center line with resolution 135 and (b) resolution 405 (remark: The D values is based on the local strain component ε xx , and the curves colors are the same for both (a) and (b)).

CHAPTER 2 .

 2 Figure 2.24: The oscillations of radial strain along the interface of model C of different methods.

Figure 2 .

 2 Figure 2.25: The MD values for model C (a) of the interface radial strain component and (b) of the interface tangential strain component (remark: the curves colors are the same for both (a) and (b)).

CHAPTER 2 .

 2 INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT

Figure 3 . 1 :

 31 Figure 3.1: Illustration of neighbour voxel average technique and layers used for averaging with different window sizes in 2D: (b) N L = 1, (c) N L = 2 and (d) N L = 3 on a bi-phase model with 81 voxels (remark: the red and blue voxels represent different phases).

Fig. 3

 3 Fig. 3.1(a) presents a bi-phase model. The red and blue voxels represent each phase, and the orange line describes the real interface. The orange voxel is the target voxel of the matrix phase that is being treated. Figs. 3.1(b), (c), and (d) show the averaging for different layer counts. As shown in these figures, only voxels of the same phase are considered for averaging (the voxels marked by a cross).

CHAPTER 3 .

 3 INTERFACE SMOOTHING TECHNIQUE AND OUR IMPROVEMENTSshown in Eq. (3.2).

Figure 3 . 2 :Figure 3

 323 Figure 3.2: Local stress (MPa) fields of model C under transverse tensile loading along x direction (〈ε xx 〉 = 1.10 -3 ): (a) FEMF Conformal mesh, (b) FFT results without average, (c) FFT results with average of N L = 2. From left to right, they are σ xx , σ y y , σ x y . (remark: the models in the figure are of exactly the same size, although the images, which are captured using different software, could suggest otherwise)

Figure 3 . 4 :

 34 Figure 3.4: Illustration of a model without composite voxel (a) and with composite voxel (b). The orange line is the real interface [163].

Figure 3 . 5 :

 35 Figure 3.5: Illustration of different mixing rules: (a) Voigt, (b) Reuss and (c) Laminate.

Figure 3 . 6 :

 36 Figure 3.6: Comparison of relative error for different mixing rules [163].

  .7, the first step is dividing the composite voxel into N sub sub-voxels. The second step is calculating the center position of each sub-voxel (x and determining the constituent (ph i ) at each position with the help of the real interface. This enables to assign the properties of the chosen constituent (C ph i ) to the corresponding sub-voxel. Finally, the number of subvoxels with the same constituent (N ph i sub ) divided by the total number of sub-voxels inside the composite voxel is the volume fraction of this constituent in the composite voxel (V ph i f = N ph i sub /N sub ).

Figure 3 . 7 :

 37 Figure 3.7: Description of the calculation of the volume fraction of each constituent inside a composite voxel following the method proposed by Kabel: the composite voxel consists of 25 sub-voxels, each is assigned to a phase based on its centroid (the real interface is shown as an orange curve and the composites voxels are shaded).

Figure 3 . 8 :

 38 Figure 3.8: Illustration of the normal vector calculation: (a) the red phase is dominant, and (b) The blue phase is dominant.

Figure 3 . 9 :

 39 Figure 3.9: Illustration of the improved composite voxel method: (a) identification of the interface, (b) calculation of the minimum distance d ph i , and (c) selection of l and determination of the composite voxels.

Figure 3

 3 Figure 3.10: The interface 2nd normal strain (a) and tangential strain (b) with N C = 375 and using FEM conformal mesh (reference result) and W-S alone or combined with the conventional and enhanced composite voxel methods. (remark: the curves colors are the same for all plots).

Fig. 3 .

 3 Fig. 3.11 clearly shows the effect of the reduction of oscillations based on the OCV and the TCV methods.

Figure 3 . 11 :

 311 Figure 3.11: The MD value of interface 2nd normal (a) and tangential (b) strains using W-S and different composite voxels approaches (remark: the curves colors are the same for all plots).

Figure 3 . 12 :

 312 Figure 3.12: The MD values of interface 2nd normal strain component (a) and of the tangential strain component (b). (remark: AVE 1L means N L = 1 and so on, and the curves colors are the same for all plots).

Figure 3 . 13 :Figure 3 . 14 :

 313314 Figure 3.13: The interface 2nd normal strain (ε r r ) component of different approaches of model C with N C = 375 (remark: AVE 1L means N L = 1 and so on).

Figure 3 . 15 :

 315 Figure 3.15: Relative errors with respect to FEMR along the layer of interface 2nd using W-S and different smoothing techniques: (a) average relative error, (b) maximum relative errors; (c) errors on the maximum strain (remark: AVE 2L means N L = 2, and the curves colors are the same for all plots).

Figure 4 . 1 :

 41 Figure 4.1: Schematic representation of the cohesive zone. The green arrows represent the distribution of tractions over the process zone [175].

Figure 4 . 2 :

 42 Figure 4.2: Bilinear traction-separation law (a) under mode I, (b) under mode II and III (compression can not lead to a cohesive zone failure).

1 (Traction), δ f 2 (

 12 .6) Meanwhile, the final displacements for the state of complete debonding, δ f Shear), δ f 3 (Shear) are obtained as:

  eq , (i = Lt ; Lc; Tt ; Tc; Zt ; Zc; Mt ; Mc) (4.16) where Lt means longitudinal tension, Lc means longitudinal compression, Tt is transversal tension, Tc is transversal compression, Zt is z-direction tension, Zc is z-direction compression, Mt denotes the matrix tension, and Mc denotes the matrix compression. This is illustrated in Fig. 4.3.

Figure 4 . 3 :

 43 Figure 4.3: The schematic of a strand yarn with directions [186].

CHAPTER 4 .

 4 REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD gradient-enhanced damage model phase-field model PDE equation rl 2 c ∆r = r dl 2 c ∆d = -g ′ (d ) Hl c /G c Table

4 . 5 .Figure 4 . 5 :

 4545 Figure 4.5: A solid body Ω with the crack set Γ: (a) sharp cracks with its modeling and (b) approximated diffuse crack bands.

Figure 4 . 6 :

 46 Figure 4.6: Illustration of diffusive damage profile of a crack at x = 0 for various length parameters l c (Miehe's method).

  can be seen that they all satisfy d = 1 when x = 0, and that d vanishes away from x = 0. This is also shown in Fig.4.7.

Figure 4 . 7 :

 47 Figure 4.7: The damage profile of three diffusive crack functions.

Fig. 4 .

 4 8 presents the profile of α (d ) with different ζ values. Since α (d ) is not monotonic in the ranges ζ < 0 and ζ > 2, they are forbidden to avoid non-physical damage diffusion with the increase of d .

CHAPTER 4 .Figure 4 . 8 :

 448 Figure 4.8: The crack geometry function α (d ) for different values of parameter ζ.

  d ) dx is the apparent displacement jump due to the opening of the crack. At the same time, with Eqs. (4.45) and (4.40), we can get:

. 48 )

 48 It is reminded that the generic geometric crack function is α (d ) = ζd + (1 -ζ) d 2 , where ζ ∈ [0, 2] to avoid non-physical damage (explained in Section 4.3.2)

. 49 )

 49 When the stress exceeds the critical value σ > σ c , the diffusive crack will localize into a localization band with its size controlled by the length scale l c . As mentioned above, the crack is limited into a zone [-D, D]. Building upon the above setting, multiplying Eq. (4.46) by d ′ and integrating with respect to x ∈ [-D, D], we get σ 2 ω (d ) -A 0 α (d )l 2 c (∇d ) 2 = 0. (4.50)

. 56 )

 56 From Eqs. (4.55) and (4.56), we can discuss three different cases: • Miehe's model (ζ = 0): In this case, σ c and the initial elastic stage do not exist. Therefore, the jump distance w(d * ) can only be described by Eq. (4.55), where the phasefield response always depends on l c . • Pham's model (ζ = 1, σ c varies): In this case, σ c and initial elastic stage exists. Based on Eq. (4.49), σ c is a function of a 1 and l c . However, a 1 = 2 in this case. Thus, the value of σ c and the phase-field response depends on the choice of l c . • Wu's model (ζ = 2, σ c is fixed): In this case, σ c and initial elastic stage exists, and σ c is fixed. The displacement jump (Eq. (4.56)) will become independent to l c . When d * = 1, the opening crack displacement reaches its maximum value, we have:

Figure 4

 4 Figure 4.9: (a) The relation between the stress σ and the opening crack displacement w for the linear softening law [208], (b) Energetic degradation function for the linear softening law with different values of a 1 [201].

(4. 71 )

 71 which is also a function of ζ. Because of the irreversibility condition, the crack bandwidth cannot shrink with increasing d which yields the following condition:D 0 (ζ) ≤ D u (ζ) . (4.72)Therefore, a relation between the initial and final semi-crack bandwidth depending on the variation of ζ is shown in Fig.4.10. As shown in the figure, for Pham's diffusive crack (ζ = 1), three softening laws can be used except for the linear case that is the most common law.

Figure 4 .

 4 Figure 4.10: Initial and ultimate semi-crack bandwidth with the variation of ζ [201].

CHAPTER 5 .

 5 SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL

CHAPTER 5 .

 5 SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL where A t n

Figure 5 . 2 :

 52 Figure 5.2: Stress-strain relations with different time increment δε.

CHAPTER 5 .

 5 SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL

Figure 5 .

 5 Figure 5.3: The composite model with two half fibers on the borders (blue part: fiber; red part: matrix).

Figure 5 . 4 :

 54 Figure 5.4: The macroscopic stress-strain curves of the unit-cell model with two half fibers on the borders with different values of l c .

CHAPTER 5 .Figure 5 . 5 :

 555 Figure 5.5: The local damage fields at final failure for different values of l f l c : (a) 12.2; (b) 2.53; (c) 1; (d) 0.5 (remark: the semi-circles with dash lines are the fiber contours).

Figure 5 . 6 :Figure 5

 565 Figure 5.6: The damage profile of the model along the center line with different l f /l c ratios.

Figure 5 .

 5 Figure 5.8: Two-fibers unit-cell composite model with l f constant (blue part: fiber and red part: matrix).

Figure 5 . 9 :

 59 Figure 5.9: The strain component ε xx of the two fibers unit cell model without the phase-field (purely elastic with 〈ε〉 = 0.005): (a) l f = 6.2.10 -4 mm, and (b) l f = 2.67.10 -4 mm (remark: all conditions are the same as the tests with phase-field).

CHAPTER 5 .Figure 5 .

 55 Figure 5.10: The local damage fields at final failure of unit-cell two fiber models with l c constant, and l f = Test 1 :(a) 6.2 . 10 -4 mm, (b) 2.67 . 10 -4 mm; Test 2: (c) 6.2 . 10 -4 mm, (d) 2.67 . 10 -4 mm (remark: the circles with dash lines are the fiber contours).

Figure 5 . 11 :

 511 Figure 5.11: The macroscopic stress-strain curves of the two fibers unit-cell model (l c fixed) shown with different values of l c : (a) Test 1 and (b) Test 2.

Figure 5 . 12 :

 512 Figure 5.12: The damage profile of the model along the center line with different l f /l c ratios: (a) Test 1 and (b) Test 2.

Figure 5 . 13 :Fig. 5 .

 5135 Figure 5.13: Two-fibers unit-cell composite model with l f fixed blue part: fiber and red part: matrix.

Figure 5 . 14 :

 514 Figure 5.14: The macroscopic stress-strain curves of the model shown in Fig. 5.8 for different values of l c .

CHAPTER 5 .Figure 5 . 15 :

 5515 Figure 5.15: The local damage fields at final failure with different l c values, and l f l c =: (a) 12; (b) 9.2; (c) 6.9; (d) 5.5; (e) 3.5; (f) 1.7 (remark: the circles with dash lines are the fiber contours).

Figure 5 .

 5 Figure 5.16: The damage profile of the model along the center line with different l f /l c ratios.

Figure 5 .

 5 Figure 5.17: The damage field in the fiber at the final failure: (a) l f c = 8.8.10 -4 mm, (b) l f c = 2.3.10 -4 mm.

Figure 5 . 18 :

 518 Figure 5.18: The damage field of the model at the final failure: (a) l f c = 8.8.10 -4 mm, (b) l f c = 2.3.10 -4 mm (Remark: the circles with dash lines are the fiber contours).

  f c = 8.8.10 -4 mm) to the border (l f c = 2.3.10 -4 mm). Meanwhile, Fig. 5.19 shows that Test 2 features a larger peak stress and final failure strain.

Figure 5 . 19 :

 519 Figure 5.19: The macroscopic stress-strain curves for two different l f c values.

Fig. 5 .

 5 Fig.5.20(a) shows the values of the history variable at different locations during loading. It can be seen that the behavior for Test 1 is normal and physical, as the energy at the center is always higher than at the border until final failure due to the strain concentration. For Test 2, however, the strain energy at all locations perfectly coincides with Test 1 in a first stage but at some point the strain energy at the border grows much faster. Fig.5.20(b) which shows the damage evolution at the same locations during loading leads to a similar conclusion. The behavior for Test 1 is also more physical while Test 2 behavior is incomprehensible.

Fig. 5 .

 5 Fig.5.20(a) shows the values of the history variable at different locations during loading. It can be seen that the behavior for Test 1 is normal and physical, as the energy at the center is always higher than at the border until final failure due to the strain concentration. For Test 2, however, the strain energy at all locations perfectly coincides with Test 1 in a first stage but at some point the strain energy at the border grows much faster. Fig.5.20(b) which shows the damage evolution at the same locations during loading leads to a similar conclusion. The behavior for Test 1 is also more physical while Test 2 behavior is incomprehensible.

Figure 5 . 20 :

 520 Figure 5.20: The history (a) and the damage (b) values at the border and the center of the model with two different l f c values (remark: the curves begin at 〈ε〉 xx = 5 . 10 -4 ).

Figure 5 . 21 :

 521 Figure 5.21: The model RVE with multi-fibers inspired from [238].

CHAPTER 5 .Figure 5 . 22 :

 5522 Figure 5.22: The local damage fields at final failure for the RVE model (resolution 225) with different values of l c =: (a) 2.9 . 10 -4 mm, (b) 4.8 . 10 -4 mm, (c) 9.5 . 10 -4 mm and (d) 1.9 . 10 -3 mm.

Figure 5 . 23 :CHAPTER 5 .

 5235 Figure 5.23: The macroscopic stress-strain curves of the RVE model (resolution 225) with different l c values.

CHAPTER 5 . 4 Table 5 . 4 :

 5454 SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL Test 1 Test 2 l c (mm) 7.00 . 10 -4 5.00 . 10 -The damage properties of different tests.

Figure 5 . 26 :

 526 Figure 5.26: The macroscopic stress-strain relations with/without heterogeneity correction.

Figure 5 . 27 :Figure 5 . 28 :

 527528 Figure 5.27: The damage profiles along the centerline of different l c values: (a) l c = 7.00 . 10 -4 mm; (b) l c = 5.00 . 10 -4 mm.

CHAPTER 5 .

 5 SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL

Figure 5 . 29 :

 529 Figure 5.29: The geometry of single notch plate (Material 1 (M1) is in yellow and Material 2 (M2) in purple).

1 6. 2

 12 Chapter overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Reminder of Wu's phase-field model . . . . . . . . . . . . . . . . . . . . . . . 6.3 The difficulty of applying Fixed-point algorithm . . . . . . . . . . . . . . . . 6.4 Applying Wu's phase-field model by Newton-Krylov algorithm in FFT solver 6.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.1 Single notch model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5.2 Asymmetric double notch model . . . . . . . . . . . . . . . . . . . . . . . 6.5.3 Single notch with a reinforced fiber . . . . . . . . . . . . . . . . . . . . . 6.5.4 Unit-cell model with two fibers with l f constant . . . . . . . . . . . . . 6.5.5 Unit-cell model with two fibers with l c constant . . . . . . . . . . . . . . 6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 . 1 :

 61 Figure 6.1: Illustration plot of the values of g ′ (d ) when d < 0.

CHAPTER 6 .

 6 DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER Initialization: d (x) is already known for previous iteration (a 0 ) Define CG solver tolerance (tol) (b 0 ) i = 0 (Initialize iteration of CG solver), and define i max (c 0 ) δd i =0 = 0 (d 0 ) Calculate RP(d , x) and LP(δd i =0 , x) (e 0 ) Calculate r = RP -LP ( f 0 ) Define a = r (g 0 ) Define d el new = r T • r While (d el new ≥ tol 2 ): (a 1 ) Calculate

  by 10 Calculate LP(δ i d , x) Calculate r = RP -LP else Calculate r = r -α.q end (g 1 ) Define d el ol d = d el new (h 1 ) Calculate d el new = r T .r (i 1 ) Define β = d el new d el ol d ( j 1 ) Calculate a = r + β.a End (6.19)

Figure 6 . 2 :

 62 Figure 6.2: The flow chart of Wu's phase-field method based on the Newton-Krylov algorithm in FFT solver (remark: CV is damage field convergence check. AltCV is alternating minimization convergence check. The part in the red rectangular is the main algorithm for solving the damage field).

Figure 6 . 3 :

 63 Figure 6.3: The geometry of single notch plate: (a) Lateral bands thickness fixed for testing mesh size convergence. (b) Lateral bands number of voxels fixed for other tests. Material 1 (M1) is in yellow and Material 2 (M2) in purple.

CHAPTER 6 .Figure 6 . 6 :

 666 Fig.6.3(b). As mentioned before, the resolution chosen here is N = 251, with a voxel size h = 4.00 . 10 -3 mm. In these tests, the smallest value used for the characteristic length is l c = 0.01mm, which results in l c /h = 2.51 and hence satisfies Miehe's criterion.

  10 4 N/mm, and σ e c = 2445.42 MPa. These tests of Wu's CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER model are accomplished using the FEM solver with a resolution N = 63 (63 × 73), l c = 0.04 and 0.05mm. The results are shown in Fig. 6.7. As it can be seen, increasing the damage properties of M2 affects significantly the macroscopic response, even though M2 has a low volume fraction and should play a little role in damage evolution as compared to M1.

Figure 6 . 7 :FFT 1 :

 671 Figure 6.7: The macroscopic stress-strain curves of M1 with G e c = 2.7 . 10 4 N/mm, and σ c = 2445.42MPa, and model resolution N = 63 using the FEM solver.

Figure 6

 6 Figure 6.9: Macro-level stress-strain relation curves of (a) contrast 1, (b) contrast 10 and (c) contrast 1000 (remark: the curve thickness of different time steps are different, which is for distinguishing the curves between each other).

Figure 6 .Figure 6 . 11 :

 6611 Figure 6.10: Illustration of Asymmetric double notch model.

Figure 6 . 12 :

 612 Figure 6.12: The damage fields of single-notch model(N = 251) under y-direction loading using different methods: Wu's model: (a) l c = 0.01 mm, (b) l c = 0.02 mm, (c) l c = 0.03 mm. Miehe's model: (d) l c = 0.01 mm, (e) l c = 0.02 mm, (f) l c = 0.03 mm.

Figure 6 . 13 :

 613 Figure 6.13: Illustration of the single notch model with fiber reinforcement.

Figure 6 . 18 :

 618 Figure 6.18: The macroscopic stress-strain curves of two fiber unit-cell models (l c = 5.10 -4 mm) with different inter fiber distance l f .

Figure 7 . 1 :

 71 Figure 7.1: The signed distance function of the first matrix voxel layer along the interface: (a) raw values and (b) values after low-pass filtering ( f p s describes the bandwidth of the FFT low-pass filter).

Figure 7 . 2 :

 72 Figure 7.2: Illustration scheme of the model with diffusive interface powered by phase-field [241].

Figure 7 . 3 :

 73 Figure 7.3: 3P bending tests of UD composites.

  1 + e Jξ 2 1 + e Jξ 3 .

	CHAPTER 2. INSIGHTFUL AND QUANTITATIVE ANALYSIS OF THE CAUSES OF
	NUMERICAL ARTIFACTS AND IMPROVED SCHEMES FOR FFT
	(2.6)

Table 2 .

 2

2 presents the mesh convergence of model C. The Relative Errors of Properties (REP) are calculated as follows

Table 2 .

 2

	.7)

2: The elastic properties (E x and G x y (MPa)) of model C of different FFT algorithms with different resolutions and their relative error compared to the FEM.

Chapter 3 Interface smoothing technique and our improvements Contents 3

.1 Chapter overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Neighbour voxels average method . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Overview and description . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Test cases and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Composite voxel method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Conventional composite voxel method . . . . . . . . . . . . . . . . . . . 3.3.3 Improved composite voxel method . . . . . . . . . . . . . . . . . . . . . 3.4 Comparisons of interface smoothing technique . . . . . . . . . . . . . . . . 3.5 Chapter conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 4 .

 4 2: The parameters of different types of diffusive cracks.

		0	d 2	2	exp -|x| l c
	Pham [200]	1	d	8/3	l c 1 -|x|	2
	Wu [207, 208]	2 2d -d 2	π	1 -sin |x| l c

  .51) CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD Integrating Eqs. (4.43), (4.44) and (4.29) into Eq. (4.51), the stress σ (d * ) can be re-written as:

Table 5 .

 5

	l c	12.2	6.1	4.1	2.5	2.1	1.5	1	0.61 0.5	0.4
	d max f	0.085 0.15 0.23 0.38 0.45 0.56 0.67 0.8 0.84 0.88

1: The maximum damage d max f

Table 5 .

 5 35. Regarding the damage properties of the constituents, CHAPTER 5. SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL two series of damage properties are tested, as shown in table 5.2. Note that it may seem that these two tests are irrelevant and not comparable because the damage properties are totally different. However, it is interesting to see that they both lead to the same conclusion. Besides, the choices of l c also satisfy the criterion, where l c /h > 2. 2: The damage properties of different tests.

		Test 1	Test 2
	G c (N/mm) f	9 . 10 10	60
	G m c (N/mm) 9.6 . 10 -4	0.003
	l c (mm)	8.8 . 10 -4 7.0 . 10 -4

Table 5

 5 .3 are tested. They both satisfy Miehe's criterion l f /h > 2, where h is the voxel size.

		Test 1	Test 2
	l	f c (mm) 8.8.10 -4 2.3.10 -4

Table 5 .

 5 3: The fiber characteristic length values for different test cases.

  [START_REF] Loos | Chapter 2 -composites[END_REF] where a fixed point algorithm can still be established. It can be summarized as follows: CHAPTER 5. SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL Initialization, if t n = 0 :

Table 5 .

 5 7: The total iteration numbers of homogeneous and heterogeneous formula.

  .11) CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY FFT SOLVER At the next time increment t n+1 , we would have H max > 1

	σ 2 c
	2

  +1d i . Because most terms in R(d ) are linear, we only have to compute g ′′ (d ).

	CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY
	FFT SOLVER
	with δd = d i We get
	
					13)
	and its derivative	R ′ (d ) =	∂R(d ) ∂d	.	(6.14)
	Eq. 6.12 can be re-written as:
		R ′ (d )δd = -R(d )	(6.15)

  2, we have studied the iteration numbers under different contrasts of damage properties. Similar tests are also presented here. Note that the model's resolution is N = 251, with a voxel size h = 4.00 . 10 -3 mm and l c = 0.02 mm. Moreover, three contrast values are presented in table 6.2.

	G e c (N/mm) 2.70 27 2700
	Contrast G e c /G c	1	10 1000
	Table 6.2: The damage properties of different tests.

  ).Next, the total iteration numbers of the CG solver are presented in Table6.3. Like Miehe's model with heterogeneous formulation, Wu's model is also sensitive to the contrast of damage properties. Thus, in the case of high contrast models, an acceleration scheme should be considered.Table6.3: The total iteration numbers until the final failure of different cases.

	CHAPTER 6. DAMAGE MODELING OF MATERIAL WITH WU'S PHASE-FIELD MODEL BY
	FFT SOLVER
			Time step
		1.00 . 10 -4 5.00 . 10 -5 2.00 . 10 -5
	Contrast 1	1.48 . 10 5	2.08 . 10 5	2.92 . 10 5
	Contrast 10	3.34 . 10 5	4.58 . 10 5	7.96 . 10 5
	Contrast 1000 3.22 . 10 6	4.17 . 10 6	6.44 . 10 6

  FFTW Fastest Fourier Transform in the West: https://www.fftw.org/. 16, 19

	3SHM Three Stages Homogenization Method. 5
	CDM Continuum Damage Mechanics. 76
	CG Conjuagate Gradient method. 144
	CGO Continuous Green Operator. 16
	CLT Classical Laminate Theory. 5
	CZM Cohesive Zone Model. 73
	DEM Digital Element Method. 10
	DFT Discrete Fourier Transform. 18
	DGO Discrete Green Operator. 29, 30, 52
	DST Domain Superposition Technique. 13
	FEM Finite Element Method. 8
	FFT Fast Fourier Transform. 2, 16
	GED Gradient-Enhanced Damage model. 79
	LEFM Linear Elastic Fracture Mechanics. 72
	MD Molecular Dynamic algorithm. 10
	micro-CT micro-Computed Tomography. 8
	MOC Method Of Cell. 5, 7
	M-T Mori-Tanaka model. 5

NNA Nearest Neighbor Algorithm. 10 ODE Ordinary Differential Equation.

[START_REF] Espadas-Escalante | Mesoscale analysis of the transverse cracking kinetics in woven composite laminates using a phase-field fracture theory[END_REF]

As shown in Fig.2.5, the Green operator of Moulinec-Suquet and Willot with centered finite difference show strong oscillations that are reduced with the green operator of Willot with the forward-and-backward scheme and the rotated scheme. Considering the forwardand-backward finite difference scheme may break the symmetry. Consequently, the Γ WI-R is chosen for subsequent analysis. (a) Γ M-S . (b) Γ WI-C . (c) Γ WI-W . (d) Γ WI-R .

(a) (b) (c) (d)

(a) (b)

CHAPTER 4. REVIEW OF DAMAGE MODELING OF MATERIAL AND PHASE-FIELD METHOD

when two cracks are close to each other [START_REF] Martínez-Pañeda | A phase field formulation for hydrogen assisted cracking[END_REF] or impinging on the hole if the crack is close [START_REF] Molnár | 2d and 3d abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF].

On the contrary, Wu's scheme is free of this dilemma as l c can always be treated as a numerical parameter and be chosen as a small positive number [START_REF] Wu | A length scale insensitive phase-field damage model for brittle fracture[END_REF][START_REF] Mandal | Length scale and mesh bias sensitivity of phase-field models for brittle and cohesive fracture[END_REF]. Furthermore, the bandwidth of Wu's model is much more narrow than the two others. However, compared to the two other phase-fields, the evolution equation of Wu's phase-field becomes nonlinear, resulting in increased complexity in solving the equation and higher computational costs.

In the last two chapters, Wu's and Miehe's models are chosen to study the impact of l c for heterogeneous materials.

CHAPTER 5. SENSITIVITY OF MIEHE'S PHASE-FIELD TO THE CHOICE OF THE L C FOR FIBROUS MATERIAL ON MICRO-LEVEL

We find the properties in [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF] which reports a σ c for the matrix of σ m c = 80Mpa, which gives a l c = 7.3 . 10 -5mm. In order to integrate this l c into the model, an adaptation is necessary to meet Miehe's criterion, so the resolution is increased to 875. Besides, l c = 3.0 . 10 -4 mm is also tested, which is the value used in [START_REF] Espadas-Escalante | A phase-field model for strength and fracture analyses of fiber-reinforced composites[END_REF].

The damage fields are presented in Fig. 5.24, and the macroscopic stress-strain curves are shown in Fig. 5.25. In all single notch model tests, the following material parameters are adopted for M1: Young's modulus: E 0 = 2.1.10 5 MPa, Poisson's ratio ν = 0.3, failure strength σ c = 2445.42 MPa, and fracture critical energy release rate G c = 2.7 N/mm. These properties result in an internal length l ch = E 0 G c σ 2 c = 0.095 mm. Considering that a 1 ≥ 3 2 is an exigence for Wu's phase-field, the condition l c ≤ 0.08mm has to be satisfied for this material. Miehe's criterion regarding the ratio between l c and the voxel size h still has to be verified for Wu's model.

Unit-cell model with two fibers with l f constant

In this part, the model utilized is the same as in section 5.4.3, and is shown in Fig. 6.16. The model size is 0.02 × 0.02mm with two fibers of diameter 6.2.10 -3 mm located at the midheight. The inter fiber distance is fixed to 2.76.10 -3 mm.

The blue part is fiber, for which the material is E-glass of properties: is much smaller than twice the distance from the fiber-matrix interface to the border. Considering the periodic boundary conditions, the strain concentration point should be at the center of model.

In Fig. 6.17, the macroscopic response of Wu's and Miehe's models, with the same model and material parameters, are presented. Viewing the macroscopic response, it can be seen that Miehe's model is l c sensitive, as increasing l c reduces the maximum stress, and moves the crack propagation forward. Wu's method, on the contrary, shows no l c -sensitivity. Based on table 4.3, it can be seen that the crack bandwidth of Wu's model is πl c , which means that damage vanishes completely outside this bandwidth interval. From l c = 2.3 . 10 -4 mm to l c = 8.0 . 10 -4 mm, the ratio l f /l c decreases from 12 to 3.5, which means the fiber is always out of the crack bandwidth, which could explain the l c insensitivity. However, for the models with one or two notches, the crack initiates close to an interface between two materials and can be affected by the second material. CHAPTER 7. GENERAL CONCLUSIONS AND PERSPECTIVES mesh, the drawbacks of the voxel mesh, and the existing solutions to overcome some of those drawbacks.

A numerical method tailored for voxel-based meshes, the Fast Fourier Transform (FFT) based method initially proposed by Moulinec-Suquet, has been presented in the chapter. We have detailed the mathematical aspects of the method to demonstrate why it has a higher computational efficiency than conventional FEM. At the end of this chapter, some drawbacks and improved schemes have been presented.

In the second chapter, we have thoroughly investigated a specific drawback of the FFT method and existing solutions to overcome it. We have focused on the presence of numerical artifacts in the form of spurious oscillations in the results obtained with the FFT method.

Among existing solutions to deal with this issue, we have focused on Discrete Green Operators (DGOs) based on finite-difference schemes or FEM, as proposed respectively by Willot and Schneider. Quantitative comparisons between FFT methods and FEM with reduced and full integration have led us to conclude that Schneider's scheme with full integration is similar to the full integration formulation of FEM, and has a better accuracy than reduced integration formulations like Willot's scheme. The full integration scheme, however, requires eight times more memory and a much greater computation time. Thus, we have concluded that Willot's scheme could be a better choice in many cases, as it had a good balance between computational accuracy and speed.

Although the different causes of the numerical artifacts had already been identified in the literature, our quantitative analysis has enabled us to rank them in terms of impact on the predictions. In general, the non-smooth zig-zag pattern of the interface caused by the voxelbased discretization appeared to be the most significant cause of oscillations that contributed more than five times than others.

Interface smoothing techniques to reduce oscillations caused by the non-smooth interface, namely the composite voxel method and the neighbor voxels average method, have been presented in the third chapter.

The neighbor voxels average method (AVE) has been shown to reduce the oscillations significantly, with peak stresses much closer to FEM conformal mesh results, and a much cleaner local field.

An improved composite voxel method (OCV ) has also been proposed in the chapter. Compared with the conventional composite voxel method TCV , OCV has been enhanced to be applicable directly to interfaces for which there is no parametric representations, such as those found in digital images. The most important aspect that we found is a better oscillations reduction effect with OCV as compared to TCV .

In general, we observed that both AVE and OCV had a good oscillation reduction effect. We felt problematic, however, that AVE is only a post-processing operation, which does not ensure that all results still satisfy equilibrium and compatibility equations. Even though AVE could be considered as easier to implement, we concluded that it should be used only when a non-intrusive approach is needed. When technically possible, our analysis led us to recommend OCV as it showed a good relative error and oscillation reduction.