Last but not least, I would like to thank my family: my

Multiscale modeling and simulation of materials with mesostructure gradients

by Minh Vuong LE Understanding mechanisms at small scales in composite structures is crucial for the design of many applications in engineering, including structures for aircraft engines which allow a reduction of CO 2 emissions. The objective of this Ph.D. is to develop numerical approaches based on the finite element method to accurately describe the thermomechanical fields in these complex structures at different scales, at reasonable computations costs.

One difficulty in such problems is the lack of scale separation, corresponding to a weak difference between the characteristic dimensions of the heterogeneities and the fluctuation length of the applied loads. Therefore, the classical homogenization methods are rather inaccurate near stress concentrations or when global gradient deformation modes such as bending occur.

In a first part, I investigate the use of a previously developed method, the filter-based computational homogenization method, in the context of woven composite structures. The accuracy of the method is tested and a comparison with classical first-order homogenization methods, including with enhanced re-localization process, is conducted.

In a second part, an original method, developed during this Ph.D., is presented. It uses advantages of both computational homogenization and domain decomposition techniques. The approach, called CMCM (Coarse Mesh Condensation Multiscale) method, is based on several ingredients. The structure is first meshed with a fine mesh, at the scale of heterogeneities. The obtained mesh can be extremely large and lead to intractable computational costs with classical Finite Element solvers. The structure is then decomposed in several subdomains. On each subdomain, preliminary calculations allow constructing a basis of local solutions, which are then condensed on unknown degrees of freedom of a coarse mesh covering the whole structure. The full resolution in the fine mesh is then substituted by solving a problem on the coarse mesh only, allowing computational times reductions by several orders of magnitude. The technique is first applied to a linear elastic structure. An error analysis is conducted, and several applications are proposed, including an industrial case whose microstructural mesh features more than a billion finite elements, or a calculation on a microstructure whose description is provided by micro tomography. CMCM is extended to unstructured and non conforming meshes. The technique is finally extended to other behaviors, such as thermo-elasticity and general nonlinear regimes. It is shown that the method allows obtaining a satisfying approximation of the local solution on the fine mesh together with a large computational cost reduction.

Résumé

Modélisation et simulation des matériaux à gradients de mésostructure par Minh Vuong LE La compréhension des mécanismes liés aux petites échelles dans les structures composites est un enjeu crucial pour développer de nombreuses applications en ingénierie, notamment les structures entrant dans la composition des moteurs d'avion permettant une réduction des émissions de CO 2 . L'objectif de cette thèse est de développer des méthodes numériques fondées sur la méthode des éléments finis pour permettre de décrire finement les champs thermomécaniques dans ces structures complexes à différentes échelles et à des coûts de calculs raisonnables.

L'une des difficultés est la non séparation des échelles, correspondant à une faible différence entre longueur caractéristique des hétérogénéités et longueurs de fluctuations des champs appliqués, rendant les méthodes d'homogénéisation classiques peut précises près des concentrations de contraintes ou en présence de gradients de déformations macroscopiques comme de la flexion.

Dans une première partie, nous évaluons les performances d'une méthode développée précédemment, la méthode d'homogénéisation numérique par filtres, pour des structures modèles en composites tissés. Une quantification des erreurs obtenues par la méthode et une comparaison avec des techniques d'homogénéisation du premier ordre avec relocalisation améliorée sont menées.

Dans une deuxième partie, une méthode originale est développée, tirant parti à la fois des techniques d'homogénéisation numériques et des méthodes de décomposition de domaine. La technique, nommée méthode CMCM (Coarse Mesh Condensation Multiscale), repose sur plusieurs étapes. La structure est tout d'abord complètement maillée à l'échelle des hétérogénéités. Les maillages obtenus peuvent alors être extrêmement larges et impraticables pour des solveurs éléments finis classiques. La structure est ensuite décomposée en sous-domaines. Sur chaque sous-domaine, des calculs préliminaires permettent de construire une base de solutions locales. Les solutions locales sont ensuite condensées sur les inconnues nodales d'un maillage grossier. La résolution de la structure complètement maillée est alors substituée par une résolution sur le maillage grossier, permettant des gains de taille de calculs de plusieurs ordres de grandeurs. La technique est d'abord testée dans un cas élastique linéaire. Une analyse quantitative des erreurs est menée, et des applications sont proposées, comme un cas industriel d'une structure composite tissée maillée avec plus d'un milliard d'éléments ou un calcul sur une image de micro structure obtenue par micro tomographie. La méthode CMCM est étendue à des maillages grossiers réguliers ou non, et pouvant être non conformes avec le maillage micro sous-jacent. La méthode est ensuite étendue et appliquée à des problèmes de structures thermoélastiques et des comportements locaux non linéaires. Il est montré que la méthode permet d'obtenir des qualités de solution satisfaisantes y compris dans des cas ne présentant pas de séparation d'échelles, avec des réductions de coûts de calcul importantes.

First and foremost, I would like to express my sincere gratitude to my Ph.D. advisor, Prof. Julien Yvonnet, for his continuous support throughout the three years of my Ph.D.; for his patience, motivation and immerse scientific knowledge. His guidance has helped me in all the time of research and writing of this thesis.

I would also like to thank my two scientific and industrial supervisors, Fabrice Detrez and Nicolas Feld. Thank you for the precious scientific discussions, they have helped me solving a lot of difficulties during my three years of Ph.D.. Thank you Fabrice for giving me a lot of your time at the lab, helping me with any kind of problem, even the smallest one. Thank you Nicolas for guiding, helping me since my very first day at Safran Tech. I really appreciate your advice in helping me improving my skills. I also want to thank my supervisors for their time and effort in proofreading this thesis.

I would like to thank my fellow labmates in the MSME lab for the interesting discussions throughout my Ph.D.. Thank you for all the amusing moments at the international conferences at Caltech, Glasgow, Giens. Thank you Daicong for including me in your article, even though I just contribute a small part.

My thanks also go to all my colleges at the CSM team at Safran Tech. Thank you for the warm welcoming and a lot of constructive comments about my work. Thank you Tonya for helping me proofreading my first article.

I also want to thank all my friends in Vietnam and France, for all the interesting conversations, relaxing after-work moments, which were very helpful after the intense research works. [START_REF] Zeman | A finite element perspective on nonlinear FFT-based micromechanical simulations[END_REF] . . 1.2 Comparison between Cosserat, Cauchy continuum and the reference solution of the y-displacement along a horizontal line of a bending structure [START_REF] Forest | Cosserat overall modelling of heterogeneous materials[END_REF]) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Second-order computational homogenization scheme [START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF]) 1.4 Variation of the component C 2222 of the sixth-order elasticity tensor in function of the volume fraction τ [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF]) . . . . . . . . . . . . . . . . 1.5 Illustration of one-dimensional basis functions [START_REF] Efendiev | Multiscale finite element methods: theory and applications[END_REF] . 1.6 Schematic description of the HMM method (Abdulle et al. (2012)) . . . . . 1.7 Illustration of different scales used in the approach [START_REF] Wangermez | Interface coupling method for the global-local analysis of heterogeneous models: A second-order homogenization-based strategy[END_REF] 1.8 Structure decompsoed into 2 subdomains Ω (1) and Ω (2) and their interface Υ [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF] Vectors and second order tensors, as well as matrices, are denoted by bold letters A. Third order tensors are denoted by calligraphic uppercase letters G, fourth-order, fifthorder and sixth-order tensors are denoted by double case letters A. Double contraction of indices for second order tensors A and B is denoted by A : B = A ij B ij , inner product for two vectors a and b by a • b = a i b i , and simple contraction of indices for a second order tensor A and a vector b is denoted by (Ab) i = A ij b j . The gradient operator is denoted by ∇(.) and the divergence operator by ∇ • (.).

List of Figures

Notations

Let u be defined as the displacement vector and x as a material coordinate, we define the linearized strain as:

ε ij = 1 2 ∂u i ∂x j + ∂u j ∂x i (1)
the second gradient displacement tensor by:

G ijk = ∂ 2 u i ∂x j ∂x k (2)
and the strain gradient tensor as:

∇ε ijk = 1 2

Introduction

With the steady increase of air traffic, airline companies have to face serious problems related to greenhouse gas emission reduction. Indeed, according to research studies from EUROCONTROL forecast, the number of flights in Europe by 2035 should be increased by 150% as compared to 2012. Despite all the developments in aerospace industry, the impact of CO 2 , toxic gases and noise in the environment is still growing. The "Advisory Council for Aviation Research and innovation in Europe" (i.e. ACARE) has set ambitious targets through Clean Sky 1 & 2 programs to reduce the environmental and societal impact of air traffic in Europe and in the world. These include 75% of CO 2 emission reduction, 90% of NO x reduction and 65% of noise reduction by 2050 as compared to 2000. The efforts range from 30% to 40% on engine manufacturers and from 40% to 50% on aircraft manufacturers. Consequently, these companies, among which Safran Group is an international leader on aircraft engines and equipment, are enforced to innovate on new aircraft structure designs to reduce harmful emissions while maintaining highest safety standards. One of the most efficient and virtuous ways of achieving these requirements is through the design of lighter engines, landing gear and nacelles.

One possible solution for structure mass reduction is the use of new materials which are more efficient than existing ones. Composite materials, and more specifically 3D woven composites, for which Safran has gained an expertise for several years, constitute an effective response to these ambitious objectives. Their properties, combining strength and lightness help creating stronger and lighter parts, for example fan blades, fan cases or landing gear rods (see e.g. Fig. 2), than their metallic predecessors. This technology plays a key role in the manufacturing of the new LEAP (Leading Edge Aviation Propulsion) commercial airplane engine, for which a 15% fuel reduction has been obtained as compared to the previous CFM56 model. There are several types of woven composites, including e.g. organic matrix-composites, which consist of ceramic fibers embedded in organic matrix (OMC), or ceramic matrix composite which are composed of ceramic fibers embedded in ceramic matrix (CMC).

Introduction

When studying 3D woven composites, three possible scales can be considered (Fig. 3). The finest scale is that of discrete carbon. The order of magnitude of heterogeneities at this scale typically lies within 1 to 100 micrometers. The intermediate scale is the meso scale where the fibers are packed into supposedly homogeneous yarns, which are woven into specific patterns. The typical scale of heterogeneities ranges here between 0.1 and 10 mm. The largest scale of study for 3D woven composites is the macro scale, which is the scale of the complete, homogeneous part. The order of magnitude of this scale should theoretically exceed a few centimeters. In this thesis, I will mainly focus on the transition between the meso and macro scales. With the constant increase in computer performance, numerical methods and simulations play a crucial role in the study and design of heterogeneous structures, and particularly those made of composites. They allow the analysis of these materials without performing expensive laboratory experiments. However, when studying 3D woven composites, many challenges must be tackled in order to obtain a representative and predictive modeling by the Finite Element Method (FEM). One reason is the strong influence of the manufacturing process on the mesostructure, and hence on the overall behavior, that characterizes them. This is the framework of the Mécanique des Composites project, a Collaborative Research Project (PRC) funded by the Direction Générale de l'Aviation Civile (DGAC), in which this Ph.D. belongs. Indeed, the process of manufacturing parts in 3D woven composites makes it possible to modify their mesostructure according to local aerodynamics and mechanical requirements. The same flexibility leads to a difficulty in the constitutive modeling of these materials, since many elementary principles underlying the formulation of homogeneous constitutive models are violated, like the separation of scales assumption between meso, macrostructure and gradients of loading, as well as the repeatability of the mesostructure that would characterize a "Representative Volume Element". None of these conditions being met, there is a strong need for taking into Introduction account the influence of said mesostructure on the overall response of the components using non-standard approaches. Several complementary approaches are developed as part of the A3 work package of the PRC: "Multi-scale dimensioning of composite parts taking into account the mesostructure and the process-induced effects".

Micro scale

Meso scale

Macro scale

Homogenization methods are powerful tools to study heterogeneous materials while maintaining the local analysis limited to a Representative Volume Element (RVE), and then avoiding the costly calculations in a fully detailed structure. Homogenization methods and multiscale approaches will be presented in more details in Chapter 1. Among the main challenges to model aerospace 3D woven composites structures, the lack of scale separation assumption mentioned above induces a severe difficulty. Indeed, under the hypothesis of scale separation, the characteristic dimensions of the mesoscopic heterogeneities (here the yarn) should be much smaller than the fluctuation length of local strain fields. This statement can be violated in 3D woven CMCs or OMCs, as the dimensions of the yarn, or at least that of the weaving pattern, might be quite large compared to the scale of the part. This effect is even more pronounced around stress/strain concentrations. If this assumption is not met, classical homogenization methods may not be valid anymore. On the other hand, solving directly a composite structure at the scale of yarns with the Finite Element Method may involve solving huge problems with more than billions of unknowns.

The main objective of this work is to develop a new numerical methodology able to:

• Obtain the whole full-field (stress and strain) solution in heterogeneous structures involving several level of scales, with or without periodic repeatability, for arbitrary scale separation assumption;

• Obtain the solution at reasonable computational costs;

• Handle arbitrary boundary conditions on the structure;

• Handle a large variety of mechanical behaviors, including linear and nonlinear ones.

The organization of this thesis is as follows. In Chapter 1, a brief summary of existing multiscale approaches for problems without clear scale separation is presented. The advantages and drawbacks of each method are pointed out to give a clear overview of the state of the art. In Chapter 2, a Filter-based computational homogenization approach, which was investigated at the beginning of the Ph.D., is presented and tested on problems which lack scale separation. The pros and cons of this method are also discussed to justify the need for another approach, which constitutes the main original contribution of this Ph.D.. In Chapter 3, the basics of this so-called Coarse Mesh Condensation Multiscale (CMCM) method are presented in the context of heterogeneous linear periodic structures. The method is described step by step, from local problems on subdomains to the construction of a relation between local and global scales, and finally the resolution of the global problem. The numerical difficulties related to the inconsistency between local and global meshes are also investigated in this chapter. Benchmark numerical examples on 2D and 3D heterogeneous structures are presented in order to demonstrate the performance of the proposed method in solving full-field linear elastic problems. The CMCM method is extended in Chapter 4 to thermoelasticity and in Chapter 5 to nonlinear problems. A detailed description of the technical implementation of the proposed method in the Python programming language is presented in Appendix B, to show the extensibility of the method to High Performance Computing (HPC) architectures.

Chapter 1

State of the art

In this chapter, I present a brief review of computational methods used to simulate heterogeneous structures such as composite parts. Starting from classical homogenization methods, I then focus on techniques that can be used in problems where the scale separation assumption is not necessarily met. This assumption states that the typical size of heterogeneities should be much smaller than the fluctuation lengths associated with macroscopic loads, i.e. at the engineering scale. These techniques include higher-order homogenization methods, multiscale methods and High Performance Computing techniques. I present the advantages and drawbacks of these approaches to justify the introduction of a new method in the following chapters.

First-order computational homogenization approaches

In this section, first-order computational homogenization approaches are presented. I will restrict this survey to computational approaches, and will not mention the classical analytical techniques which have been described e.g. in [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF]; [START_REF] Milton | Theory of composites. Cambridge monographs on applied and computational mathematics[END_REF]; [START_REF] Buryachenko | Micromechanics of heterogeneous materials[END_REF]; [START_REF] Bornert | Homogenization in mechanics of materials[END_REF]; [START_REF] Li | Introduction to micromechanics and nanomechanics[END_REF]; [START_REF] Auriault | Homogénéisation de phénomènes couplés en milieux hétérogènes[END_REF]; [START_REF] Dvorak | Micromechanics of composite materials[END_REF]; [START_REF] Suquet | Continuum micromechanics[END_REF].

FEM linear approaches

As mentioned in the Introduction, computational homogenization methods are nowadays widely used in both academic and industrial domains in order to study the behavior of composite materials. The first work related to computational homogenization was introduced in [START_REF] Adams | Transverse normal loading of a unidirectional composite[END_REF], where the solutions at the micro scale and the macro properties were calculated through numerical simulations. The main principles of the homogenization problem using finite elements were introduced in [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics, Homogenization Techniques for Composite Media[END_REF]. Since then, many research works have been done on the development of the so called first-order computational homogenization method, see for example references in [START_REF] Geers | Computational homogenization of structures and materials[END_REF][START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF]; [START_REF] Geers | Multiscale modeling of microstructure-property relations[END_REF]. First-order homogenization is not only applicable to elasticity but also to many other problems like thermoelaticity [START_REF] Temizer | Homogenization in finite thermoelasticity[END_REF][START_REF] Temizer | On the asymptotic expansion treatment of two-scale finite thermoelasticity[END_REF]), piezoelectricity Tich ỳ et al. (2010); [START_REF] Yang | Special topics in the theory of piezoelectricity[END_REF]; [START_REF] Brenner | Numerical computation of the response of piezoelectric composites using Fourier transform[END_REF]; [START_REF] Pettermann | A comprehensive unit cell model: a study of coupled effects in piezoelectric 1-3 composites[END_REF]; [START_REF] Berger | Unit cell models of piezoelectric fiber composites for numerical and analytical calculation of effective properties[END_REF]), poroelasticity Ly (2015) or linear viscoelasticity [START_REF] Tran | A simple computational homogenization method for structures made of linear heterogeneous viscoelastic materials[END_REF], among many others. The first-order method has been successfully applied to evaluate the effective elastic properties of random media in [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF]. Several methods for solving linear computational homogenization problems are reviewed in [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]. Even though local fields can be reconstructed in first-order homogenization through the so-called "localization" process, first-order homogenization methods lack in representing boundary effects and localized strain/stress Chapter 1. State of the art fields [START_REF] Kruch | Computation of coarse grain structures using a homogeneous equivalent medium[END_REF]. For this reason, other techniques mentioned in the following of this chapter will be introduced.

FEM nonlinear homogenization approaches

Nonlinear homogenization introduces critical additional difficulties, mainly due to the fact that the superposition principle no more holds. A first class of approaches can be found in the literature under names such as "Concurrent Multiscale Methods", "Multilevel Finite Element", and "Computational Homogenization" (see e.g., [START_REF] Smit | Prediction of the large-strain mechanical response of heterogeneous polymer systems: local and global deformation behaviour of a representative volume element of voided polycarbonate[END_REF]; [START_REF] Feyel | Multiscale FE 2 elastoviscoplastic analysis of composite structure[END_REF]; [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF]; [START_REF] Terada | A class of general algorithms for multi-scale analyses of heterogeneous media[END_REF]; [START_REF] Kouznetsova | An approach to micro-macro modeling of heterogeneous materials[END_REF]; [START_REF] Ghosh | A multilevel computational model for multiscale damage analysis in composite and porous media[END_REF]; [START_REF] Feyel | A multilevel finite element method (FE 2 ) to describe the response of highly non-linear structures using generalized continua[END_REF]; [START_REF] Geers | Multiscale first-order and second-order computational homogenization of microstructures towards continua[END_REF]; [START_REF] Yvonnet | The Reduced Model Multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains[END_REF]; [START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF], and popularized under the name of FE 2 as introduced in [START_REF] Feyel | Multiscale FE 2 elastoviscoplastic analysis of composite structure[END_REF]. In this method, the main idea is to solve a heterogeneous nonlinear structure problem by means of a macroscopic mesh, whose strain state defines boundary conditions on a representative volume elements (RVE) associated with a lower scale, which after solving the local problem provides the macro stress through averaging. One appealing feature of FE 2 is that the local (fine scale) problems can be solved in parallel, as being independent of each others. FE 2 has been widely applied and extended to many problems including second-order homogenization, transient conduction problems, model reduction, damage localization or topology optimization, among many others [START_REF] Kouznetsova | Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF]; [START_REF] Yvonnet | The Reduced Model Multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains[END_REF]; [START_REF] Ozdemir | Computational homogenization for heat conduction in heterogeneous solids[END_REF]; [START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF]; [START_REF] Coenen | A multi-scale approach to bridge microscale damage and macroscale failure: a nested computational homogenization-localization framework[END_REF]; [START_REF] Xia | Concurrent topology optimization design of material and structure within FE 2 nonlinear multiscale analysis framework[END_REF]. Other studies on nonlinear behavior can be found in [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF]; [START_REF] Matouš | A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials[END_REF].

Fourier Transform approaches

Fourier solvers have become efficient tools in solving heterogeneous microstructures. A very well known method that can be listed in this category is the Fast Fourier Transformbased approach (FFT-based homogenization). The very first idea of FFT-based homogenization has been introduced in Moulinec andSuquet (1994, 1995), using an iterative algorithm to compute the local and overall response of heterogeneous structures from their microstructural images. The approach allows the expansion of elastic solution into Neumann series, whose convergence depends not only on the mechanical properties of each component, but also on the choice of the reference medium [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF]. As an illustration, considering the balance equation in linear elasticity at small strains:

∇ • (C(x) : ε(x)) = 0 (1.1)
Expressing the elasticity tensor into a constant reference medium C 0 and a complement as C(x) = C 0 + C(x) -C 0 , and expanding the strain tensor into a constant part ε and a fluctuation ε as ε(x) = ε + ε(x) (1.1) yields:

∇ • C 0 : ε(x) = -∇ • τ (ε(x)) (1.2) with τ(x) = C(x) -C 0 : ε(x)
is the polarization stress tensor. Solution of (1.2) can be obtained in the form of the so-called Lippmann-Schwinger equation as:

ε(x) = ε -Γ 0 * τ (ε(x)) , (1.3) 1.2

. Gradient and nonlocal homogenization approaches

where Γ 0 is the Green function associated with the linear operator in (1.2). This equation can be solved iteratively, by alternatively expressing variables in the Fourier space and in the real space domain, using a fixed point iterative algorithm does not require solving a matrix system. The conditions that ensure the convergence of the Neumann series have been introduced in [START_REF] Michel | A computational scheme for linear and non-linear composites with arbitrary phase contrast[END_REF]; [START_REF] Milton | Theory of composites. Cambridge monographs on applied and computational mathematics[END_REF]. There have been many research works exploiting the FFT-based homogenization method, for example [START_REF] Monchiet | Numerical homogenization of nonlinear composites with a polarization-based FFT iterative scheme[END_REF]; [START_REF] Kabel | Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations[END_REF]; [START_REF] Monchiet | Combining FFT methods and standard variational principles to compute bounds and estimates for the properties of elastic composites[END_REF]; [START_REF] Schneider | FFT-based homogenization for microstructures discretized by linear hexahedral elements[END_REF]. Recently in [START_REF] Brisard | FFT-based methods for the mechanics of composites: A general variational framework[END_REF], a novel FFT-based scheme based on the Hashin-Shtrikman variational principle [START_REF] Hashin | On some variational principles in anisotropic and nonhomogeneous elasticity[END_REF] has been proposed using the conjugate gradient method for the resolution of the problem. Despite its efficiency in determining the overall properties of heterogeneous structures, the FFT-based homogenization method is only applicable to periodic structures over regular grids (Fig. 1.1), therefore it is usually applied to microstructures obtained from X-ray tomography images and with distinct scale separation. The method has been verified by solving a two-phase laminate with inelastic ratedependent phases using a Newton scheme and conjugate gradient algorithm. The approach has recently been extended to finite strain problems in De Geus et al. (2017).

Gradient and nonlocal homogenization approaches

In this section, homogenization approaches taking into account the gradient at the macro scale are presented. These methods are often used for highly heterogeneous structures, where the separation of scales cannot be considered.

Cosserat and Generalized continuum media

Many research works have been proposed in the past to solve problems of heterogeneous structures without scale separation. For instance, a Cosserat-type generalized continuum media approach has been proposed in [START_REF] Forest | Cosserat overall modelling of heterogeneous materials[END_REF] and later applied in [START_REF] Feyel | A multilevel finite element method (FE 2 ) to describe the response of highly non-linear structures using generalized continua[END_REF] for nonlinear cases, to replace the classical homogenization method. The principle Chapter 1. State of the art consists in including additional degrees of freedom or higher order gradients of displacement, taking into account non-separated scale effects at the micro scale. The boundary conditions to be prescribed on a unit cell are derived from the macro displacement field. The heterogeneous Cauchy medium is then replaced by a homogeneous Cosserat continuum. The transition between micro and macro scales is constructed by minimizing the distance of kinematics between two scales, such distance L is described as the following on a considered unit cell Ω (Forest and Sab (1998)):

L = Ω (u(x) -u(x) -Φ(x) × (x -x)) 2 dΩ (1.4)
where x and x are the position vectors at micro and macro scale, respectively; u(x) is the displacement field at micro scale; u(x) and Φ(x) are the displacement and rotation fields at macro scale, respectively. From the results obtained by a numerical example presented in Forest and Sab (1998) (see Fig. 1.2), it can be seen that this approach did increase the precision of the structural response. However, taking into account more localized loads or nonperiodic heterogeneous structures is challenging for such techniques.

Second-order computational homogenization

Another well known approach for non-separated problems is the second-order homogenization initially introduced in [START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF][START_REF] Kouznetsova | Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy[END_REF]Kouznetsova et al. ( , 2004a)). This is an extension of the first-order (classical) homogenization method, whose scheme is depicted in Fig 1 .3. The main idea of the approach is to impose the macroscopic deformation tensor and its gradient on a representative volume element (RVE). In order to construct a boundary value problem on an RVE, the starting point is to build the transformation of a material vector x at the micro scale from the macro scale x using a Taylor series expansion truncated at second-order term:

x = F M • x + 1 2 x • 3 G M • x + x (1.5)
where F M is the deformation gradient tensor, 3 G M is the second-order deformation gradient tensor; x is the extra term taking into account the fluctuations of the micro field. After solving the local problem on a RVE, the effective behavior of the structure is determined using the Hill-Mandel energy condition Suquet (1985a); [START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF]. As a result, the macro scale is a full-gradient higher-order continuum, whose macroscopic stress tensor is calculated by averaging the local stress tensors at the micro scale, and the higher-order stress tensor is defined as follows:

3 Q M = 1 2V 0 V 0 (P c m x + xP m ) dV 0 (1.6)
where P m is the first Piola-Kirchhoff stress tensor at the micro scale.

In contrast to the first-order linear homogenization approach, at second-order, the size of the RVE does have an effect on the determination of the macroscopic behavior. The role of the size of RVE has been investigated in [START_REF] Kouznetsova | Size of a representative volume element in a second-order computational homogenization framework[END_REF]. In addition, a more general procedure for the enforcement of the boundary conditions on the RVE problem can be found in [START_REF] Kaczmarczyk | Scale transition and enforcement of RVE boundary conditions in second-order computational homogenization[END_REF].

In Yvonnet et al. (2020a); [START_REF] Monchiet | Strain-gradient homogenization: a bridge between asymptotic expansion and quadratic boundary condition methods[END_REF], new boundary conditions for the second-order homogenization method, consistent with asymptotic expansion approximations (see next section), have been proposed to evaluate the fully anisotropic fifth and sixth-order tensors in linear second-order homogenization by finite elements. Finally, an extension of strain gradient computational homogenization coupled to electrical conduction phenomena (called as flexoelectricity) has been recently proposed in [START_REF] Yvonnet | Apparent Flexoelectricity due to Heterogeneous Piezoelectricity[END_REF].

It can be seen that the second-order homogenization is able to properly assess the second-order constitutive equations emerging from a lower scale, and to naturally pass deformation gradients to a RVE, which allows to take into account the effects of heterogeneities when scales are not separated [START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF]. As an application, in Yvonnet et al. (2020a), the response of anisotropic materials taking into account the effective strain gradient behavior has been studied using the second-order homogenization approach. Other research works on second-order micromorphic theory or stress-gradient theory can be found in [START_REF] Forest | Finite-deformation second-order micromorphic theory and its relations to strain and stress gradient models[END_REF]; [START_REF] Hütter | Kinematics and constitutive relations in the stress-gradient theory: interpretation by homogenization[END_REF].

Chapter 1. State of the art

Asymptotic methods

In [START_REF] Tran | A micromechanics-based approach for the derivation of constitutive elastic coefficients of strain-gradient media[END_REF], the authors have proposed an asymptotic homogenization approach for linear elasticity taking into account the strain gradient effects at the macro scale. The method consists in using the series expansion at the local scale while keeping higherorder terms to account for the microstructural effects due to the non-separation of scales in certain situations. The total displacement field at the local scale reads:

u(x, y) = U(x) + X 1 (y) : E(x) + 2 X 2 (y) . . .G(x) + 3 X 3 (y) :: D(x)... (1.7)
where U is the macroscopic displacement field, whereas E(x), G(x) and D(x) are the macroscopic strain, strain gradient and double strain gradient, respectively; is the scale factor, which is the ratio between the characteristic length of the micro scale and the one of macro scale under loading; the X i (y) tensors can be obtained by solving the local problems on the RVE. The micro-macro transition is determined using the elastic energy density, where the higher-order terms are kept to take into account the non-separated scale effects. The efficiency of this approach has been demonstrated through a periodic 2D reinforced fiber composite (see e.g. Fig. 1.4). However, the method has only been applied to academic periodic structures.

Nonlocal homogenization

Another category of methodology is based on nonlocal elasticity theories. Usually, these methods use a nonlocal operator to define the behavior, which depends on a nonlocal definition of the strains. For example, the Eringen model [START_REF] Eringen | On nonlocal elasticity[END_REF] (see also [START_REF] Gao | Modeling plasticity at the micrometer scale[END_REF]) for a more recent asymmetric theory of nonlocal elasticity) defines the effective stress as the spatial convolution of the strain with a nonlocal elastic operator:

σ(x) = Ω C(x -x ) : ε(x )dx .
(1.8)

Multiscale methods

In (1.8), C(xx ) is a nonlocal elasticity tensor, e.g. the product of the local elasticity tensor with a monotonically decreasing function such as a Gaussian or bell-shaped function defined in a domain Ω.

These types of nonlocal elastic models are however empirical and not founded on a microstructural analysis. Furthermore, taking into account holes or cracks in such models constitutes a tough issue [START_REF] Polizzotto | Nonlocal elasticity and related variational principles[END_REF]. Non-local homogenization has been used for damage problems [START_REF] Fish | Computational damage mechanics for composite materials based on mathematical homogenization[END_REF] or wave propagation Fish et al. (2002b,a); [START_REF] Hui | A nonlocal homogenization model for wave dispersion in dissipative composite materials[END_REF] in heterogeneous materials.

In Yvonnet and Bonnet (2014a,b); [START_REF] Tognevi | A multi-scale modeling method for heterogeneous structures without scale separation using filter-based homogenization scheme[END_REF], the authors proposed a nonlocal homogenization method based on numerical filters such as Gaussian filters, where the averaging operators are replaced by nonlocal operators. This method allows defining a homogenized model for a given characteristic wavelength, independently of scale separation. This technique will be further described and applied in chapter 2.

Multiscale methods

In this section, I present a non exhaustive list of multiscale approaches, which are strongly based on the coupling between the macro scale and a lower scale (micro or meso).

Multiscale Finite Element Method

The multiscale finite element method (MsFEM) is also known to be a very effective approach when dealing with heterogeneous structures [START_REF] Efendiev | The multiscale finite element method (MsFEM) and its applications[END_REF]; [START_REF] Efendiev | Multiscale finite element methods: theory and applications[END_REF]. The objective of the approach is to solve partial differential equations with multiscale solutions Hou andWu (1997, 1999); Hou et al. (1999); [START_REF] Efendiev | The multiscale finite element method (MsFEM) and its applications[END_REF][START_REF] Efendiev | Accurate multiscale finite element methods for two-phase flow simulations[END_REF]; [START_REF] Chen | A mixed multiscale finite element method for elliptic problems with oscillating coefficients[END_REF]; [START_REF] Aarnes | An efficient domain decomposition preconditioner for multiscale elliptic problems with high aspect ratios[END_REF]; [START_REF] Wu | Analysis of upscaling absolute permeability[END_REF]. The main idea of MsFEM is to capture the information at the micro scale through multiscale basis functions and then those basis functions are coupled using a global numerical formulation. The choice of a local basis function can be arbitrary (see e.g. Fig. 1.5). Chapter 1. State of the art Despite the good performance of the approach, the treatments to reduce discontinuities across interfaces induce discontinuities in the displacement field in the reconstructed solution at the fine scale.

Heterogeneous Multiscale Method

Another multiscale method, namely the Heterogeneous Multiscale Method (HMM) has been proposed by Abdulle et al. Abdulle et al. (2012). The solution of the method involves a set of macroscopic variables U and a set of microscopic variables u. The schematic description of the method is shown in Fig. 1.6. In this framework, D is the sought part of the model depending on the problem type. For example, if it is a mechanic problem, D can be the constitutive relation between the stress and strain; or if it is a thermal problem, D might be the effective conductivity tensor. HMM estimates the sought data by using the microscale model at each point where the data is needed. For this purpose, the microscale model has to be constrained so that its macro state is similar to that which we are interested in:

f (u, d(U)) = 0 (1.9)
where d(U) is the constraint of the microscale model. In practical terms, HMM requires two main components:

• A macroscopic solver, for example a finite element solver

• An algorithm to determine the sought data D using the micromodel, involving the constrained microscale simulation and the post-processing of the micro solutions.

The HMM method has shown great promise in solving heterogeneous structures [START_REF] Engquist | The heterogeneous multi-scale method[END_REF]; [START_REF] Ming | Analysis of the heterogeneous multiscale method for elliptic homogenization problems[END_REF]; [START_REF] Weinan | Heterogeneous multiscale method: a general methodology for multiscale modeling[END_REF]; [START_REF] Abdulle | The finite element heterogeneous multiscale method: a computational strategy for multiscale PDEs[END_REF]; [START_REF] Zhang | Analysis of the heterogeneous multiscale method for parabolic homogenization problems[END_REF]; [START_REF] Ren | Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics[END_REF]; [START_REF] Abdulle | The finite element heterogeneous multiscale method: a computational strategy for multiscale PDEs[END_REF]Bai (2013, 2012). However, several difficulties still need to be addressed, for example the fluctuations of the microstructure when the scales are not separated.

Multiscale methods

Global-local analysis

Many research works have been proposed in the literature to deal with heterogeneous structures using global-local techniques. For example, in [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF], the authors have proposed a sub structuring technique, namely the Arlequin method. The Arlequin method offers a framework of mixing and gluing different models; the main principles of the method are:

• A superposition of mechanical states in a subzone S inside an entire domain Ω

• An energy distribution between the states in the subzone S using weight functions for the conservation of local energies

• A weak and compatible gluing of these states in S

The method allows to locally take into account defects and low-scale mechanical phenomena of heterogeneous structures. The Arlequin approach can also be used for contact problems [START_REF] Dhia | Contact in the Arlequin framework[END_REF], non-linear problems [START_REF] Hu | Multi-scale nonlinear modelling of sandwich structures using the Arlequin method[END_REF], or many other applications [START_REF] Bauman | On the application of the Arlequin method to the coupling of particle and continuum models[END_REF].

Another well known method in global-local analysis is the bridging domain method for coupling continua with molecular dynamics [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF]. In this approach, the continuum and the molecular domains are overlapped in a bridging subdomain. The reason that this method is considered as a multiscale method is that the spectra and resolution of the continuum model have much smaller cutoff frequencies than the molecular dynamic model. Many applications of this approach has been done in the literature, for instance [START_REF] Xu | Conservation properties of the bridging domain method for coupled molecular/continuum dynamics[END_REF] The method requires solutions of the local problems on subdomains or on RVE, by prescribing boundary conditions similar to the second-order homogenization [START_REF] Kouznetsova | Multi-scale constitutive modeling of heterogeneous materials with gradient enhanced computational homogenization scheme[END_REF]. The boundary displacement W γ of the RVE is expressed as:

W γ = u * G + (F M -I d ) .dX + G M : ∆X G ⊗ dX + 1 2 G M : dX ⊗ dX + ω W M γ (1.10) with ∆X G = X G -X C
where X C and X G are the centroid of the RVE and the right interface of a rectangular RVE.

The results have shown that the method exhibits good performance compared to Mortar coupling technique for compatible models. However, the approach requires the periodicity of the RVE. Theoretically, the method is able to solve the problem at any zone of the microstructure. However, solving the problem at each zone requires an iteration scheme with non-negligible computational cost. Therefore it is preferable for solving only certain zones of interest of the microstructure.

Others

Recently, a Wavelet-Reduced Order Model for homogenization of heterogeneous structures has been proposed in [START_REF] Van Tuijl | Multi-dimensional wavelet reduction for the homogenisation of microstructures[END_REF]. This approach consists of 2 main steps: i) construction of a reduced order model in order to reduce the dimensionality of the microstructural model ii) construction of a wavelet-reduced order model in order to reduce the cost of the numerical integration. One of the advantages of this approach is the direct control of the internal force balance approximation, which is defined by the imposed wavelet tolerance. The method has been applied so far on academic periodic structures with a reduction of 10% to 30% of the stored history parameters and material models.

Other multiscale techniques for heterogeneous structures can be found in various research works. For instance, in [START_REF] Hautefeuille | A multi-scale approach to model localized failure with softening[END_REF], the authors have proposed a multiscale method for modelling localized failure with softening; in [START_REF] Daghia | A micro-meso computational strategy for the prediction of the damage and failure of laminates[END_REF], a coupling micro-meso technique has been presented for the prediction and failure of laminates. Other research works related to multiscale techniques can be found in [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF]; [START_REF] Kelley | Mesh requirements for the analysis of a stress concentration by the specified boundary displacement method[END_REF].

High Performance Computing approaches

In this section, High Performance Computing (HPC) approaches for large heterogeneous structures are presented. Such approaches are not necessarily related to multiscale calculations, but they are often used to solve complex heterogeneous structure containing a large number of unknowns, which is one of the objectives of this Ph.D..

Domain decomposition methods

One of the HPC approaches that are well known for solving large heterogeneous structures is the Balanced Domain Decomposition (BDD) method [START_REF] Tallec | Domain decomposition methods for large linearly elliptic three-dimensional problems[END_REF]. The method consists in decomposing the structure into non-overlapping, unstructured subdomains of arbitrary shape and solving problems at interfaces using the associated trace operator and a preconditioned conjugate gradient approach. Supposing that we have two subdomains Ω (1) and Ω (2) that share one interface (see e.g. Fig. 1.8). The main 1.4. High Performance Computing approaches idea of the method is to solve finite element problems on each subdomain satisfying the following condition:

u 1 = u 2 on Υ (1.11)
where u 1 and u 2 are the displacement solutions on subdomains Ω (1) and Ω (2) , respectively.

FIGURE 1.8: Structure decompsoed into 2 subdomains Ω (1) and Ω (2) and their interface Υ (Gosselet and Rey ( 2006))

In order to solve the local problem on each subdomain, local condensed operators can be used, for example the primal Schur complement, the dual Schur complement or the hybrid Schur complement [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF].

The Finite Element Tearing and Interconnecting approach is also a non-overlapping domain decomposition method. However, unlike the BDD approach, it consists in using interface tractions as the main unknown [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]:

σ (1) n (1) + σ (2) n (2) = 0 on Υ (1.12)
The BDD and FETI methods were initially developed for homogeneous structures, and later extended to homogeneous structures using preconditioners Rixen and Farhat (1999) and initialization [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF]. However when dealing with structures containing large heterogeneities and high contrast of phase properties, these methods often exhibit poor performance. More specifically, when the subdomain interface intersects the heterogeneities, the interface problem lacks an efficient preconditioner. Several techniques have been proposed to alleviate these issues such as FETI-Geneo (FETIgeneralized eigenvalues in the overlaps) [START_REF] Spillane | Solving generalized eigenvalue problems on the interfaces to build a robust two level FETI method[END_REF], multipreconditioned FETI method (MPFETI) [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF] or an Adaptive MultiPreconditioned algorithm (AMPFETI) Spillane (2016); [START_REF] Bovet | Adaptive multipreconditioned FETI: scalability results and robustness assessment[END_REF]. BDD and FETI-based methods have been widely applied for solving large scale problems, for example [START_REF] Klawonn | Dual-primal FETI methods for linear elasticity[END_REF]; [START_REF] Li | FETI-DP, BDDC, and block Cholesky methods[END_REF]; [START_REF] Klawonn | Robust FETI-DP methods for heterogeneous three dimensional elasticity problems[END_REF]; [START_REF] Pechstein | Analysis of FETI methods for multiscale PDEs[END_REF]; Rixen et al. (1999); [START_REF] Farhat | FETI-DPH: a dual-primal domain decomposition method for acoustic scattering[END_REF]; Farhat and Mandel (1998); Farhat et al. (1998); [START_REF] Li | A dual-primal FETI method for incompressible stokes equations[END_REF]; [START_REF] Klawonn | FETI-DP methods with an adaptive coarse space[END_REF]; [START_REF] Prakash | A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics[END_REF].

Algebraic multigrid method

Other parallel strategies have been proposed, such as the algebraic multigrid (AMG) method [START_REF] Stüben | Algebraic multigrid (AMG): experiences and comparisons[END_REF]; Ruge andStüben (1984, 1987); [START_REF] Stüben | A review of algebraic multigrid[END_REF], developed to solve efficiently large linear systems of equations. Such an approach is a purely algebraic Chapter 1. State of the art matrix-based approach used to deal with large sparse linear systems and can be considered as a solver to various finite element discretization problems without any geometrical or physical background. The main principle of AMG lies in the coarsening of a given linear system of equations using a coarse grid in order to reduce the problem size.

In a classical AMG algorithm, a coarse space is a subspace defined by an interpolation matrix, whose dimension is a fraction of that of the finer space [START_REF] Xu | Algebraic multigrid methods[END_REF]. There are many variations of multigrid algorithms beside the classical AMG, for example aggregation-based AMG [START_REF] Notay | Aggregation-based algebraic multilevel preconditioning[END_REF], bootstrap AMG Brandt et al. (2011) or adaptive reduction-based AMG [START_REF] Maclachlan | Adaptive reduction-based AMG[END_REF]. However, they always share the following important steps: [START_REF] Arbenz | A comparison of eigensolvers for large-scale 3D modal analysis using AMG-preconditioned iterative methods[END_REF]. However, this is a purely algebraic approach and often viewed as a solver more than a physical approach.

Domain decomposition method based on material regularization

Another domain decomposition method for heterogeneous structures has been proposed by Zhodi et al. [START_REF] Zohdi | A domain decomposition method for bodies with heterogeneous microstructure basedon material regularization[END_REF]; [START_REF] Zohdi | A method of substructuring large-scale computational micromechanical problems[END_REF]. The main idea of the approach is to solve boundary value problems on non-overlapping subdomains whose union forms the entire microstructure. The boundary conditions of the subdomain problems can be displacement-type or traction-type, which are obtained by solving an inexpensive auxiliary boundary value problem on a regularized microstructure.

The errors between reference and approximated solutions are calculated by:

σ -σL,IS 2 E -1 (Ω) = N ξ ∑ ζ=1 Γ ξ (error in traction) • (displacement jump)ds (1.13) u -ũR,ID 2 E -1 (Ω) = N ξ ∑ ζ=1 Γ ξ (traction jump) • (error in displacement)ds (1.14)
where Γ ξ is an interior subdomain interface and N ξ is the total number of subdomain interfaces. The difference between this approach and the BDD approach is that the boundary value problems on subdomains are computed using classical homogenization boundary conditions. Another technique, the Large Time Increment (LATIN) method [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF], is a non-incremental iterative computational strategy where local (nonlinear) problems solutions in subdomains are updated through a global linear problem through interfaces conditions and an appropriate search direction to ensure the global convergence.

Conclusion

Conclusion

In this chapter, I have presented a state-of-the-art of computational approaches to model heterogeneous structures. First-order homogenization methods fail to accurately reproduce boundary and size effects in finite structures.

As shown above, High Performance Computing approaches such as FETI, BDD or LaTIn have high potential in solving heterogeneous structures, both with or without scale separation. However, these approaches are often viewed as solvers for direct calculations of the structure. In addition, HPC methods require iterative solvers and can be very expensive, especially when the structure is highly heterogeneous.

Other techniques are proposed in the literature which make use of multiscale principles, such as multiscale methods or gradient and nonlocal homogenization approaches. Gradient homogenization is an extension of the first-order homogenization, where the higher order strain and stress at the macro scale are taken into account to capture the effects of non-separated scales and macroscopic load gradients. However, these approaches are restricted to periodic structures and their performance has not been demonstrated for industrial composite structures. On the other hand, multiscale approaches directly use the information at the lower scale to solve the problem at the scale of the structure. However, they still do not fully take into account the non-separation of scales (HMM method) or produce discontinuities in the macro scale solution (MsFEM). The global-local approach has the best potential in capturing micro scale effects, but this method is preferable for solving a limited zone of interest inside the structure.

The objective of this Ph.D. is to analyze complex heterogeneous structures without assumptions on scale separation and allowing to take into account strong localized loads and size effects, at reasonable computational costs. To avoid limitations of homogenization methods, and with the aim to develop a technique that still manages to differentiate macroscopic and mesoscopic responses when relevant, a new method will be developed in the next chapters.

Chapter 2

Filter-based computational homogenization

In this chapter, a first method, initially proposed in Yvonnet and Bonnet (2014a,b); [START_REF] Tognevi | A multi-scale modeling method for heterogeneous structures without scale separation using filter-based homogenization scheme[END_REF], is investigated to analyze fiber-reinforced composites in the context of nonseparated scales. The method can be viewed as a "non-local" homogenization approach (see Chapter 1). Even though the method was not originally proposed in this Ph.D., I have first investigated it as a candidate for our study. More specifically, several original tests using the Filter-based computational homogenization are conducted here to quantitatively analyze the accuracy of the method in conditions of non-separated scales. It is also compared with first order homogenization without and with enhanced re-localization techniques [START_REF] Kruch | Homogenized and relocalized mechanical fields[END_REF]. I first describe the method and then apply it to composite structures.

General description

The main idea of the filter-based homogenization method is to construct a non-local model taking into account the influence of the heterogeneity at the microscopic scale. For this purpose, the strain and stress fields at the micro and macro scales are related through a numerical filter based on finite element shape functions, used to cut off the micro-fluctuations. The averaging operator used in the first-order homogenization is replaced by filtering operators related to each node of the coarse mesh. The filter-based homogenization method consists of 4 main steps as follows:

(i) Resolution of the local problems on a RVE using an iterative scheme (ii) Determination of non-local macro behavior of the structure using filtering operators (iii) Resolution of the global problem on a coarse mesh (iv) Relocalization of the local strain and stress fields

Definition of numerical filter and filtered macroscopic fields

First, I define the linear operator to filter the microscopic fields. It is based on the leastsquare minimization, of the deviation of said fields with respect to a solution defined over the finite elements shape functions.

Chapter 2. Filter-based computational homogenization Micro mesh Macro mesh FIGURE 2.1: Fine and coarse mesh of a RVE

Definition of numerical filters

Let us note the micro strain field ε(x), defined on a fine mesh composed of N micro nodes. The main objective of the filter is to determine the nodal values ε(x), on a coarse mesh composed of N macro nodes (Fig. 2.1), minimizing the distance between the given micro field and the sought macro field. Such distance is defined as the following:

U[ε p ] = N micro ∑ m=1 N macro ∑ p=1 M p (x m ) ε p -ε (x m ) 2 (2.1)
where ε p is the nodal strain associated to node p of the coarse mesh.

Minimizing (2.1) requires:

dU dε q ij = 0; q = 1, 2, ..., N macro (2.2)
Substituting (2.1) in (2.2) we have:

2 N micro ∑ m=1 M q (x m ) N macro ∑ p=1 M p (x m )ε p ij -ε ij (x m ) = 0 q = 1, 2, ..., N macro (2.3) ⇒ N macro ∑ p=1 N micro ∑ m=1 M p (x m )M q (x m ) ε p ij = N micro ∑ m=1 M q (x m )ε ij (x m ) (2.4)
Then the values of the strain field at the macro scale (the coefficients ε ij ) are determined by solving the following linear system:

Au = b (2.5)
where:

A pq = N micro ∑ m=1 M p (x m )M q (x m ); b q = N micro ∑ m=1 M q (x m )ε ij (x m ); u = ε 1 ij , ε 2 ij , ..., ε N macro ij (2.6)
with M p (x m ) and M q (x m ) the values of the macro mesh shape functions at the nodal coordinates of the micro mesh.

2.2. Definition of numerical filter and filtered macroscopic fields

Definition of macroscopic fields

We recall the strain and stress fields at the microscopic scale: ε(x) and σ(x), while ε(x) and σ(x) are fields at the macroscopic scale. The following relations connect these fields:

σ(x) = F h [σ(x)] (2.7) and ε(x) = F h [ε(x)] (2.8)
where F(.) is the linear filter operator that has properties related to the characteristic wavelength h of the macroscopic scale: Fig. 2.3 illustrates the filtered fields in 1D with different contrasts between micro and macro scale lengths, using the shape-function-based filter presented in the previous section. It can be seen that, when the scale contrast decreases, i.e. the macro mesh size is Chapter 2. Filter-based computational homogenization closer to that of the micro mesh, the filtered field starts to capture the fluctuation effects of the micro field.

lim h→0 F h [ε(x)] = ε(x) (2.9) and lim h→∞ F h [ε(x)] = ε(x) (2.

Local problems on a RVE

In this section, I study the RVE problem including the computation of the tensor of localization and the identification of the non-local constitutive law on a macro mesh.

Governing equations

We consider a RVE whose volume is noted Ω and its contour ∂Ω. The localization problem on the RVE is therefore expressed as follows: For a given strain field ε(x), find the compatible field ε(x) which satisfies:

     ∇ • σ(x) = 0 ∀x ∈ Ω σ(x) = C(x) : ε(x) ∀x ∈ Ω F (ε(x)) = ε(x) ∀x ∈ Ω (2.11) (2.12) (2.13)
where F is the linear filter operator described in the previous section.

The micro field consists of the superposition of the filtered field and of a fluctuation, according to the relation:

ε(x) = ε(x) + ε(x) (2.14)
where ε(x) is the filtered part and ε(x) is the fluctuation part.

Applying the filter F to both sides of (2.14) results in:

ε(x) = F [ε(x)] + F [ ε(x)]
(2.15)

In addition, we have:

F [ε(x)] = ε(x) (2.16)
We can deduce:

F [ ε(x)] = 0 (2.17)
Suppose that the macroscopic strain is derived from the field of macroscopic displacements u(x) by:

ε ij (x) = 1 2 ∂u i (x) ∂x j + ∂u j (x) ∂x i (2.18)
The displacement field u(x) is interpolated from nodal displacements u p i on the macro mesh:

u i (x) = N macro ∑ p=1 M p (x)u p i (2.19)
where M p (x) are the shape functions associated with the nodes of the macro mesh.

Then:

ε ij (x) 1 2 N macro ∑ p=1 ∂M p (x) ∂x j u p i + ∂M p (x) ∂x i u p j
(2.20)

Local problems on a RVE

Equation (2.20) can be rewritten:

ε ij (x) N macro ∑ p=1 ε ijk (x)u p k (2.21)
where

ε ijk (x) = 1 2 N macro ∑ p=1 ∂M p (x) ∂x j δ ik + ∂M p (x) ∂x i δ jk (2.22)
where δ is the Kronecker symbol.

Boundary conditions

To enforce the last non-trivial condition in (2.13), we introduce an auxiliary strain field e(x) which should satisfy:

ε(x) = e(x) -F [e(x
∇ • (C(x) : [e(x) -F (e(x))]) = -∇ • (C(x) : ε(x)) ∀x ∈ Ω F [e(x) -F (e(x))] = 0 ∀x ∈ Ω (2.24) (2.25)
Condition (2.17) implies that:

ε(x) = 0 (2.26)
Condition (2.26) being satisfied for any value of e(x) , we choose:

e(x) = 0 (2.27)
This condition is verified by two types of conventional boundary conditions:

u e (x) = 0 ∀x ∈ ∂Ω (2.28) or u e (x) periodic ∀x ∈ ∂Ω (2.29)
where u e is the field of displacements associated with e(x) according to the relation:

e ij (x) = 1 2 ∂u e i (x) ∂x j + ∂u e j (x) ∂x i (2.30)
In this work, to be able to apply the method to non-periodic meshes of 3D woven composites, I chose the first type of boundary condition. The localization problem on the RVE becomes:

Find e(x) satisfying:

∇ • (C(x) : [e(x) -F (e(x))]) = -∇ • (C(x) : ε(x)) ∀x ∈ Ω u e (x) = 0 ∀x ∈ ∂Ω (2.31) (2.32)
Chapter 2. Filter-based computational homogenization

Iterative scheme

The presence of the filter F in the left term of equation (2.31) induces a difficulty: The associated stiffness matrix will be full. To facilitate the calculation, I propose an iterative scheme starting with an initial solution e 0 (x) = 0. The convergence of this scheme has been shown in Yvonnet and Bonnet (2014a). At each iteration n, we seek the field e n+1 (x) satisfying:

∇ • C(x) : e n+1 (x) = ∇ • (C(x) : F (e n (x))) -∇ • (C(x) : ε(x)) (2.33)
until convergence, for an error ∆ defined as:

∆ = (u e ) n+1 -(u e ) n (2.34)
At convergence, the microscopic strains is obtained by:

ε(x) = ε(x) + e n+1 (x) -F e n+1 (x) (2.35)
The weak form associated with the equation (2.33) reads:

Find (u e ) n+1 ∈ H 1 (Ω) satisfying: Ω ε (u e ) n+1 (x) : C(x) :ε (δu(x)) dΩ = Ω F ε (u e ) n (x) : C(x) : ε (δu(x)) dΩ - Ω ε(x) : C(x) : ε (δu(x)) dΩ ∀δu ∈ H 1 0 (Ω) (2.36)
We can rewrite the equation (2.36) in the form of a linear system as follows:

Ku = f + f n (2.37)
where K is the stiffness matrix:

K = Ω B T (x) C(x) B(x) dΩ (2.38)
B is the matrix of derivatives of finite element shape functions; and f and f n are the volume force vectors. These forces are respectively calculated from the macro field ε(x) and the filtered field F (e n (x)):

f = Ω B T (x) C(x) [ε(x)] dΩ (2.39) f n = Ω B T (x)C(x)[e(x)]dΩ (2.40)
where [ε(x)] is the vector form of the strain tensor ε(x), and [e(x)] is the vector form of the filtered strain field F (e n (x)).

2.4. Global problem on a coarse mesh

Calculation of the localization tensor

According to the principle of superposition, the solution of problem 2.12 can be expressed as follows:

ε ij (x) = N macro ∑ p=1 D p ijk (x)u p k (2.41)
where D p (x) is considered as the localization tensor ("equivalent" to A(x) in the classical numerical homogenization scheme). In 2D, the equation (2.41) can be rewritten as:

[ε(x)] = N macro ∑ p=1     D p 11 (x) D p 12 (x) D p 21 (x) D p 22 (x) D p 31 (x) D p 32 (x)     u p 1 u p 2 (2.42)
where [ε(x)] is the vector form of the second order strain tensor ε(x).

The first column of D p (x) is calculated from a prescribed macroscopic strain field ε ij (x), as per equation (2.21), with u p = [1 0] T , and the second column of D p (x) is calculated with u p = [0 1] T .

Non-local constitutive law

Substituting equation (2.41) in the constitutive equation of problem 2.12 we obtain:

[σ(x)] = N macro ∑ p=1 C(x) D p (x) u p (2.43)
where [σ(x)] is the vector form of the stress tensor σ(x), C(x) is the matrix form of the elasticity tensor C(x). Applying filter F to both sides of the equation (2.43) gives:

[σ(x)] = N macro ∑ p=1 F [C(x) D p (x)] u p = N macro ∑ p=1 G p (x) u p (2.44)
where :

G p (x) = F [C(x) D p (x)] (2.45)
We can note that G p (x) is "equivalent" to the effective elasticity tensor C of the classical numerical homogenization scheme.

The number of elementary problems has increased compared to the classical numerical method. A macro mesh containing N macro nodes requires the resolution of N macro problems times their dimension. But, all these elementary problems are independent; one is thus able to calculate them in parallel to reduce the computing time.

Global problem on a coarse mesh

We consider a macroscopic structure of domain Ω and boundary ∂Ω (Fig. 2.4). The boundary is composed of two complementary and disjoint parts, the part where Dirichlet We therefore have the problem to be solved:

                   ∇ • σ(x) = 0 ∀x ∈ Ω [σ(x)] = N macro ∑ p=1 G p (x)u p ∀x ∈ Ω u(x) = u d ∀x ∈ ∂Ω d σ(x) • n = f n ∀x ∈ ∂Ω n (2.46) (2.47) (2.48) (2.49)
Using classical procedure, one obtains the weak form associated with the above problem.

Find u ∈ H 1 (Ω) such that:

Ω [σ(u)] [ε(δu)] dΩ = Ω n F • δu dΓ + Ω d f • δudΓ ∀δu ∈ H 1 0 (Ω) (2.50) Or: N ∑ k=1 Ω k ∑ p∈Ω k G p (x) u p B(x) δu p dΩ = Ω n F • δu dΓ + Ω d f • δu dΓ ∀δu ∈ H 1 0 (Ω)
(2.51) Introducing classical FEM discretization in (2.51), a linear system follows:

K u = F (2.52)
where:

K = N ∑ k=1 ∑ p∈Ω k Ω k B T (x) G p (x) dΩ (2.53) F = ∂Ω d N T f dΓ + ∂Ω n N T F dΓ (2.54)
and where N is the matrix of finite element shape functions.

Reconstruction of local fields

Reconstruction of local fields

After the calculation of the macro structure, the microscopic strain and stress fields are relocalized by the following relations:

[ε(x)] = N macro ∑ p=1 D p (x)u p [σ(x)] = N macro ∑ p=1 C(x)D p (x)u p
(2.55)

(2.56)

In the next section, the results obtained by the filter-based method are compared with first-order homogenization using classical or enhanced relocalization as introduced in [START_REF] Kruch | Homogenized and relocalized mechanical fields[END_REF].

Classical relocalization is expressed as follows:

[ε(x)] = A(x)[ε] (2.57) and [σ(x)] = C(x)A(x)[ε] (2.58)
where A(x) is the localization tensor obtained by solving local problem on a RVE and [ε] is the homogeneous macro strain field, considered to be constant over the whole RVE.

The difference between classical and enhanced relocalization is that the macro strain field used to relocalize the micro fields in (2.57) or (2.58) is no longer homogeneous. It is replaced by a gradient field obtained by interpolating the strains of the macro structure on the RVE. The first step consists in repositioning each RVE on the analyzed domain (Fig. 2.5). The macroscopic strain will be interpolated on the RVE. We can see in Fig. 2.6 that instead of a homogeneous field, a gradient field is imposed on the elementary volume.

Macro mesh of the structure Micro mesh of the RVE The local strain and stress field are reconstructed using enhanced relocalization as follows:

[

ε(x)] = A(x)[ε i (x)] (2.59)
Chapter 2. Filter-based computational homogenization and

[σ(x)] = C(x)A(x)[ε i (x)] (2.60)
where [ε i (x)] is the macro strain field interpolated from the coarse mesh onto the RVE (Fig. 2.6c).

Numerical examples

I have implemented the filter-based numerical homogenization method presented in this chapter in a self-contained Matlab code, using linear triangular and tetrahedral elements in 2D and 3D, respectively. The approach is also developed for non-structured meshes, especially at the macro scale, to describe more complex structures. The code is implemented using parallel computing in order to reduce computational cost, in particular at the local scale where the number of local problems depends on the number of degrees of freedom of the macro mesh. The results of the filter-based method below are obtained with this Matlab code. The reference solution in the first example is calculated using the finite element method in Matlab whereas, in the second example, it is obtained using the Z-set computing software at Safran Tech. We define a global energy error as follows:

Err E = Ω ε re f (x) -ε num (x) : C(x) : ε re f (x) -ε num (x) dΩ Ω ε re f (x) : C(x) : ε re f (x)dΩ (2.61)
where ε re f denotes the strain field obtained by the reference solution, and ε num denotes the corresponding approximated strain field obtained from the homogenization schemes.

To alleviate numerical difficulties, I restrict the numerical examples to quasi-static and linear regime under the assumption of small perturbations.

Tension test of a 2D woven composite beam

We study the tension of the woven composite beam described in Fig. 2.7.a. The considered structure consists of six repeated RVE. The fine mesh of the reference structure contains 345,257 linear triangular elements, corresponding to 346,188 degrees of freedom. The fine mesh of the RVE contains 57,255 linear triangular elements, corresponding to 57,698 degrees of freedom. The coarse mesh of an RVE used for the macro calculation is composed of 732 linear triangular elements, corresponding to 812 degrees of freedom. The mechanical properties of each constituent of the composite are defined as follows:

-The matrix is isotropic elastic (see The strain-stress relation is defined as:

        ε 11 ε 22 ε 33 2ε 12 2ε 13 2ε 23         =            1/E 11 -ν 12 /E 22 -ν 13 /E 33 0 0 0 -ν 12 /E 11 1/E 22 -ν 23 /E 33 0 0 0 -ν 13 /E 11 -ν 23 /E 22 1/E 33 0 0 0 0 0 0 1/G 12 0 0 0 0 0 0 1/G 13 0 0 0 0 0 0 1/G 23                    σ 11 σ 22 σ 33 σ 12 σ 13 σ 23        
(2.62)

These properties are defined in frame oriented by the main direction of the woven threads. The global errors of different components of the stress and strain tensors using the 3 homogenization schemes are summarized in Table 2.3. It can be seen that using the enhanced relocalization scheme does not significantly reduce global errors. On the other hand, using the filter-based homogenization scheme reduced global errors significantly, especially in terms of ε 12 , where the error is decreased by more than 50%.

3 point bending of a woven composite coupon

We define the problem shown schematically in Fig. 2.11. The structure is a small coupon of dimensions L × b × h (i.e. 50 × 10 × 1.1 mm 3 ), subjected to a 3-point bending test. The mesh of a periodic subdomain at the meso1 scale is generated with the TexGen software [START_REF] Lin | Modelling and simulating textile structures using TexGen[END_REF], on a domain with dimensions 2 × 2 × 0.22 mm 3 , using linear tetrahedral elements (Fig. 2.12). The parameters used to generate this representative unit cell are described in Table 2.4. The entire structure is composed of 25 × 5 × 5 periodic subdomains. Each subdomain contains 463,462 linear tetrahedral elements corresponding to 251,769 degrees of freedom. The macro mesh used in this example is composed of 16,600 trilinear hexahedral elements corresponding to 60,600 degrees of freedom. The matrix and woven threads are considered linear elastic isotropic and linear elastic orthotropic, their mechanical properties are listed in Table 2.5 and 2.6, respectively. Young's modulus E (MPa) Poisson's ratio ν 4000 0.3 To avoid managing a contact algorithm and to preserve the linearity of the study, the loading is simplified in the form of pressure fields. These fields are applied to each of the three contact lines: x c 1 = 3 mm, x c 2 = 25 mm and x c 3 = 47 mm respectively, and governed by the quadratic form:

p i (x) = p 0 i 1 - x -x c i x p

2

(2.63) with x p their half-width (1 mm, see Fig. 2.13) and p 0 i maximum pressure amplitude (with p 0 2 = 20 MPa and p 0 1 = p 0 3 = p 0 2 /2). The support at each point being defined by [x c ix p , x c i + x p ], the pressure is always positive. In addition, the displacement according to z of points located in (x = 0, z = 0) and (x = L, z = 0) are blocked. Theoretically, it was necessary to block the movement according to z of only 3 points to have an unique solution. In practice however, a problem subjected to a balance of forces as precarious as this one leads to important stress concentrations due to numerical approximations. This is why this solution was preferred, despite the slight hyperstatism it implies. It can be seen in the Fig. 2.14 that despite a global calculation on a coarse mesh, the filter-based method was able to partially capture the mesostructure effects of the weaving pattern on the skin of the structure, even on macro fields, which cannot be obtained by first-order homogenization. In Fig. 2.15, it can be noted that the filter-based method is quite sensitive to the mesostructure and the meshes, with a very perturbed relocalized stress field.

Conclusion

In this chapter, I have investigated the use of the filter-based computational homogenization method introduced in Tognevi et al. ( 2016) for the linear elastic analysis of woven composite structures. For that purpose, I have developed an in-house Matlab code. Furthermore, I have extended the method to unstructured meshes to study complex, realistic microstructures, and have implemented it in 3D. I have compared the method with classical first-order homogenization and an enhanced relocalization procedure in various situations, including a 3D woven microstructure. I have shown that, in the case of non-separated scales (bending, concentrated loads...), significant error reduction can be obtained as compared to first-order homogenization. However, even though a significant reduction of macroscopic quantities is obtained, errors on relocalized mesoscopic fields remain non-negligible (around 30 %) and large computational costs are associated. The main reasons are the cost of local problems (requiring an iterative solution procedure due to the non-local operator) and the number of local problems to be solved to construct the macro non-local operator, especially in 3D. In addition, a sensitivity to the macro mesh has been noticed. For all these reasons, a new original method has been developed, which is presented in the next chapters of this thesis.

Chapter 3

Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

In this chapter, I introduce a new multiscale computational approach called Coarse Mesh Condensation Multiscale (CMCM) method. It aims at solving a heterogeneous structure through a solution basis obtained from subdomain problems. In this chapter, I restrict the method to linear problems and show that in this case it leads to an efficient oneiteration scheme where local solutions (in subdomains that compose the structure) can be solved in parallel through off-line calculations. Here, the solutions of fine scale subdomain problems are related to the degrees of freedom (or dofs) of a coarse mesh which has much fewer degrees of freedom than the finer problem (by several order of magnitudes). This chapter is adapted from our published papers Le et al. (2020a,b).

General description

The key idea of the CMCM method is to solve the structural problem taking into account the microstructural kinematics. This local information is obtained by off-line calculations on subdomains, which can be performed in parallel. First, the heterogeneous structure is decomposed into subdomains Ω α , α = 1, 2, ..., N sub where N sub denotes the number of subdomains, and such that Ω = α Ω α (see Fig. 3.1 (a)). Then, each subdomain is discretized with a mesh of finite elements associated with the fine scale. In addition, a coarse mesh is defined to discretize the structure, independently of the heterogeneity (see Fig. 3.1 (b)). Dirichlet boundary conditions are then defined over each subdomain, and associated with a small number of parameters, gathered in a vector g α (see section 3.2). We define these parameters such that they minimize the error in the least square senses between the strain field in the subdomain and the strain field in the elements of the coarse mesh covering the subdomain. In the linear case, a linear relationship can be established between displacements of the coarse mesh and the strain field in the subdomains. Note that this approximation depends on the location of point x in the element of the coarse mesh and is then richer than the simple linear combination of subdomain elementary solutions (see details in section 3.2.1). The solution in the subdomains can be condensed at the nodes of the coarse mesh and the problem can be solved on this coarse mesh. The fine scale solution is obtained by post-processing. The main steps are summarized below and described in more details in the following.

(i) Solve localization problems in each subdomain Ω α (parallel computations).

Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures (ii) Determine the relationship between the dofs of the coarse mesh u E and the reduced vector g α in each subdomain Ω α (parallel computations).

(iii) Find u E by solving the global problem on the coarse mesh.

(iv) Re-localize fine scale fields from the global solution.

u i (x) = ∑ k g α k d α ik (x) ∀x ∈ ∂Ω α (3.1)
where g α is the reduced vector. The size of g α is assumed to be small compared to the total number of dofs on the boundary ∂Ω α in order to perform a reduced condensation of the internal dofs.

In this work, I adopt the following particular form for relation (3.1):

u i (x) = ε ij x j + 1 2 G ijk x j x k ∀x ∈ ∂Ω α (3.2)
where ε ij is a macroscopic strain tensor and G ijk is a macroscopic, third-order, second gradient of displacement tensor. It can also be re-expressed as a function of an effective 3.2.

Step (i): Localization problems (homogeneous) first gradient of the strain tensor ∇ε ijk through (see e.g. [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF]; [START_REF] Auffray | A complete description of bidimensional anisotropic strain-gradient elasticity[END_REF]):

G ijk = ∇ε ijk + ∇ε ikj -∇ε jki (3.3)
In addition, we restrict the approximation (3.2) to the x j x k , (j = k) terms to reduce the number of parameters, i.e. products of x 2 i terms are set to zero. The corresponding obtained solutions will be referred to as Second-order CMCM in the examples (section 3.6). For comparison, another approximation of (3.2) keeping only the linear terms (i.e. setting all G term products to zero) will be considered and referred to as CMCM in the examples.

Special attention must be paid to the case of a homogeneous subdomain, i.e. containing only one material phase. In that situation, we must ensure that the boundary conditions (3.2) do not induce spurious fluctuations, i.e. that the solution in Ω α derives from (3.2) in the form:

ε(x) = ε + ∇ε • x ∀x ∈ Ω α (3.4)
Considering that one subdomain contains a single material phase with elastic properties C 0 , applying the Hooke's law and taking the divergence of (3.4) one obtains:

∇ • C 0 : ε(x) = ∇ • C 0 : [ε + ∇ε • x] = f (3.5)
which is verified if

f i = C 0 ijkl ∇ε klj (3.6)
Hence, body forces must be prescribed in addition to the boundary conditions to ensure null fluctuations in the case of a homogeneous subdomain. One possible choice for C 0 in the general (non-homogeneous case) is (see Yvonnet et al. (2020a)) C 0 ijkl = C, where C is the homogenized elastic modulus defined by:

C = 1 V Ω α C(x) : A(x)dΩ (3.7)
and where A(x) is the classical localization tensor. The localization problem is defined as follows:

Given g α = {ε, ∇ε}, find ε(x) in Ω α such that:

         ∇ • (C(x) : ε(x)) = f ∀x ∈ Ω α f i = C ijkl ∇ε klj ∀x ∈ Ω α u(x) = ε • x + 1 2 G : x ⊗ x ∀x ∈ ∂Ω α (3.8) (3.9) (3.10)
With this definition, if the subdomain is homogeneous, the solution is self-balanced for the above quadratic boundary conditions. Problem (3.8)-(3.10) is then classically solved by finite elements as follows.

The weak form of the problem to be solved on the subdomain Ω α is given by: Find u ∈ H 1 (Ω α ) satisfying the boundary conditions (3.10) and such that:

Ω α ε(u) : C(x) : ε(δu)dΩ = Ω α f(x) • δudΩ ∀δu ∈ H 1 0 (Ω α ) (3.11)
Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

where H 1 and H 1 0 are the usual Sobolev vector spaces. Introducing classical FEM discretization in (3.11), we obtain a linear system in the form:

K α u = F α (3.12) with K α = Ω α B e T (x) C(x) B e (x) dΩ and F α = Ω α N e T (x) f dΩ (3.13)
where N e and B e are matrices of shape functions and of shape functions derivatives on the subdomain's fine mesh, respectively; and C(x) is the matrix form associated with C(x). It is worth noting that, for one subdomain Ω α , the linear systems to be solved for each component of g α involve the same K α matrix, which thus may be assembled only once.

Subdomain crossing interfaces

In case the subdomain boundaries cross the interfaces between matrix and inclusions, or the inclusion's position is close to the subdomain's edges, large errors can be generated by the boundary conditions (3.10). In that case, an extension of the method consist in prescribing the boundary conditions on the external boundary ∂ Ω α of an extended subdomain Ω α (see Fig. 3.2), while the localization tensor is only evaluated in the initial subdomain Ω α . We define a parameter β = h ext /L 0 , where L 0 is the initial length of the subdomain Ω α (considered as square in this work for the sake of simplicity). Note that this process is similar to the oversampling technique in MsFEM [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF].

hext L 0 FIGURE 3.2: Illustration of an extended subdomain

It is worth noting that in the case where the inclusion's position is close to the subdomain's edges, periodic implementation of the boundary conditions (3.10) can also be effective in reducing errors. However, I want to find a more general technique which can be applied to different cases. Therefore, the extended-subdomain technique is preferred.

The inconvenient of this extended-subdomain technique is that it increases the local problem size. For example in 3D case, as maximum extended-subdomain problem is 27 times larger than the initial one. Alternatives to this technique are possible, for example constructing an error corrector at the boundary of the subdomain or adding more deformation modes to the subdomain. However, due to time constraints, these alternatives are kept for future works.

Step (i): Localization problems

Reduced vector g α in the 2D case

In 2D, after cutting terms in x 2 i , (3.10) leads to:

u 1 (x) = ε 11 x 1 + ε 12 x 2 + 1 2 G 112 + G 121 x 1 x 2 (3.14) u 2 (x) = ε 12 x 1 + ε 22 x 2 + 1 2 G 212 + G 221 x 1 x 2 (3.15)
After simplifications, we obtain:

u 1 (x) = ε 11 x 1 + ε 12 x 2 + ∇ε 112 x 1 x 2 (3.16) u 2 (x) = ε 12 x 1 + ε 22 x 2 + ∇ε 221 x 1 x 2 (3.17)
In the Second-order CMCM case, the boundary conditions of each Ω α subdomain are parameterized by 5 independent scalar values gathered in the reduced vector:

g α = [ε 11 , ε 22 , ε 12 , ∇ε 112 , ∇ε 221 ] T (3.18)
As the problem to be solved over Ω α is linear, applying the superposition principle allows the strain field to be expressed in Ω α as:

[ε(x)] = A α (x) g α ∀x ∈ Ω α (3.19)
where [ε(x)] is the column-vector notation of the second-order strain tensor ε = 1 2 ∇u + ∇ T u and A α (x) is the localization matrix (see section 3.2.1).

In each Gauss point of the fine mesh, the matrix A α (x) is therefore a 3 × 5 matrix whose columns [ε i (x)] are the strain vector solution of (3.8)-(3.10) for g α i = 1 and g α j = 0, i = j:

A α (x) = [ε 1 (x)], [ε 2 (x)], [ε 3 (x)], [ε 4 (x)], [ε 5 (x)]
(3.20)

In the first-order CMCM case, only the first three terms of g α , and consequently the first three columns of A α , are computed.

Reduced vector g α in the 3D case

In 3D, using (3.10) and (3.3), we obtain, after simplifications and cutting terms in x 2 i :

u 1 (x) = ε 11 x 1 + ε 12 x 2 + ε 13 x 3 + ∇ε 112 x 1 x 2 + ∇ε 113 x 1 x 3 + G 123 x 2 x 3 (3.21) u 2 (x) = ε 12 x 1 + ε 22 x 2 + ε 23 x 3 + ∇ε 221 x 1 x 2 + ∇ε 223 x 2 x 3 + G 213 x 1 x 3 (3.22) u 3 (x) = ε 13 x 1 + ε 23 x 2 + ε 33 x 3 + ∇ε 331 x 1 x 3 + ∇ε 332 x 2 x 3 + G 312 x 1 x 2 (3.23)
Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

In the Second-order CMCM case, g α contains 15 independent parameters (and in the first-order CMCM case, only the first 6 are needed): Thus 15 (resp. 6) linear problems need to be solved on each subdomain. Moreover all these problems can be solved in parallel. The matrix A α (x) is in that case a 6 × 15 (resp. 6 × 6) matrix.

g α = ε 11 , ε 22 ,
It is worth noting that in the case where the strain gradient parameters are applied, the subdomain needs to be centered in order to obtain the correct deformation mode.

Step (ii): Relation between coarse mesh-displacements and subdomain boundary conditions

To relate the solutions in the subdomains to the coarse mesh, a relationship between the vector of boundary conditions g α of subdomain Ω α and the dofs in each element Ω E of the coarse mesh, denoted by u E , must be established. For this purpose, we minimize the distance between the strain approximation provided by (4.49) and the strain in each element of the coarse mesh individually (see Fig. 3.1), which is given by:

[ε(x)] = B E (x) u E (3.24)
where B E is a matrix of shape functions derivatives of one element Ω E of the coarse mesh.

Then the problem to be solved is given by

g α = Arg min J (3.25) with J = Ω E ([ε(x)] -[ε(x)]) 2 dΩ = Ω E A α (x) g α -B E u E 2 dΩ (3.26)
where Ω E is the domain associated with element E in the coarse mesh, and (a) 2 = a • a, a being a vector. Minimizing J with respect to g α gives (writing

g α i ≡ g i , A α ij ≡ A ij , B e ij ≡ B ij and u E inc ≡ u i for alleviating notations): ∂J ∂g m = 0 m = 1, 2, ..., N g (3.27) or; Ω E 2A ij ∂g j ∂g m A ij g j -B ij u j dΩ = 0 (3.28)
Then:

Ω E A im (x)A ij (x)dΩg j = Ω E A im (x)B ij (x)u j dΩ m = 1, 2, ..., N g (3.29)
We obtain a linear system of equations in the form:

H αE g αE = L αE u E (3.30) 3.4.
Step (iii): Global problem and then:

g α = (H α ) -1 L αE u E (3.31)
with:

H αE = Ω E (A α (x)) T A α (x) dΩ (3.32)
and

L αE = Ω E (A α (x)) T B E (x) dΩ (3.33)
where the size of H α is the size of g α squared; and the size of L α is the size of g α times Nd E with Nd E is the number of degrees of freedom of a macro element E.

Finally, we obtain the approximation of the fine scale strain in each element Ω E of the coarse mesh within a subdomain Ω α , using (4.49) and (5.18), as:

[ε(x)] = T αE (x) u E (3.34) with T αE (x) = A α (x) H αE -1 L αE (3.35)
where the size of T αE (x) is 3 × Nd E in 2D and 6 × Nd E in 3D.

The calculations that are performed on the fine mesh, namely offline calculations are shown in Algorithm 1.

1: LOOP over all subdomains Ω α (parallel) 2:

-Solve (3.8)-(3.10) to obtain A α (x) and store 3:

-Compute H αE using (3.32) and store

4:

LOOP over all macro elements E ∈ Ω α 5:

-Compute L αE using (3.33) and store 6:

END LOOP 7: END LOOP

Algorithm 1: Offline calculations.

Step (iii): Global problem

We now consider the problem to be solved on the coarse mesh over the entire structure (see Fig. 3.1 (b)). Assuming zero body forces for conciseness sake only, the governing equations of the problem in the structure are given by:

           ∇ • σ(x) = 0 ∀x ∈ Ω σ(x) = C(x) : ε(x) ∀x ∈ Ω u(x) = u d ∀x ∈ ∂Ω d σ • n = f n ∀x ∈ ∂Ω n (3.36) (3.37) (3.38) (3.39)
where u d and f n are prescribed displacements and tractions on the corresponding Dirichlet and Neumann boundaries ∂Ω d and ∂Ω n . Using vector forms of strain tensors, the weak form of the system of equations (3.36)-(3.39) can be expressed as follows:

Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

Find u ∈ H 1 (Ω) such that:

Ω [ε(u)] T C(x) [ε(δu)] dΩ = ∂Ω n f n • δu dΓ ∀δu ∈ H 1 0 (Ω) (3.
δu E T ∑ α Ω α T αE (x) T C(x) T αE (x) dΩ u E = δu E T ∂Ω n N T (x) f n dΓ ∀δu E ∈ H 1 0 (Ω) (3.42)
Owing to the arbitrariness of δu E T , we then obtain a linear system of equations in the form:

K u = f (3.43)
where:

K = ∑ α Ω α T αE (x) T C(x) T αE (x) dΩ (3.44) f = ∂Ω n N T (x) f n dΓ (3.45)
and N(x) are the classical finite element shape functions of elements in the coarse mesh.

In practice, K is assembled from elementary matrices computed in each micro element e covered by macro element E defined in its associated domain Ω E of the coarse mesh (see Fig. 3.3). The Dirichlet and Neumann boundary conditions are classically prescribed on the coarse mesh. Once u is known, the strain field on the fine mesh in each subdomain can be reconstructed using (3.34) and the stress using:

[σ(x)] = C(x) T αE (x) u E (3.46)
The algorithm for solving the coarse mesh problem is summarized in Algorithm 2.

Handling incompatible micro and coarse meshes

1: LOOP over all coarse mesh elements E (parallel) 2:

-Find micro elements e ∈ Ω E 3:

-Find the subdomain α contaning e 4:

-Given H αe , L αE (x):

5:

LOOP over all micro elements e ∈ Ω E 6:

-Compute T αE using (3.32); (3.33) and (3.35) 7:

-

Compute K E = Ω E T αE (x) T C(x)T αE (x)dΩ 8: Assemble K E in K 9:
END LOOP 10: END LOOP Algorithm 2: Online calculations. Integration of L α in (3.33) requires special attention in the case when one subdomain is cut by the boundary of a coarse mesh element Ω E (see Fig. 3.4). In that case, a consistent integration requires subdividing the cut micro elements Ω e for integration purpose, as shown in Fig. 3.4 (b). Noting the set of micro elements S 0 and the set of subdivided micro elements S 1 , a consistent integration of L α reads:

Handling incompatible micro and coarse meshes

L α = ∑ e∈S 1 ∑ k [A α (x k )] T B E (x k ) w k V e (3.47)
where x k and w k denote the Gauss point positions and weights, respectively, and V e is the volume of the micro element in Ω e . In the present work, I used the Matlab© command "intersect" to construct the set of subdivided elements. This algorithm returns the geometric intersection of 2 polygons as inputs. An illustration of this process is shown in Fig. 3.5.

Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures However, this operation can be costly and burdensome in the 3D case. Alternatively, an approximated integration scheme is proposed as

L α ∑ e∈S 0 ∑ k [A α (x k )] T B E (x k ) V e w k (3.48)
where S 0 is the original set of micro elements. In this integration scheme, the local values of A α (x) are affected according to the position of Gauss points in the micro elements (see Fig. 3.4 (c)), but no subdivision of micro elements is performed for the sake of simplicity.

Then here, the interface between two macro elements Ω E is not explicitly described. The errors induced by this approximated integration scheme will be evaluated in the section 3.6.3.

Numerical examples

In this section, I analyze and demonstrate the potential of the method through both academic and simplified industrial examples of realistic scale. In the context of periodic subdomains, the local problem (3.8)-(3.10) has to be solved only once (see section 3.2). In the following, three solutions are considered for comparison:

• A reference solution obtained by a direct finite element calculation, where the structure is fully meshed, covering all the details of the heterogeneities. The reference solution of the largest 3D problem is obtained with the iterative AMPFETI method Spillane (2016); [START_REF] Bovet | Adaptive multipreconditioned FETI: scalability results and robustness assessment[END_REF].

• A solution using the present method without using strain gradient (noted CMCM).

• A solution using the present method with strain gradient (noted second-order CMCM).

A global energy error is defined as:

Err E = Ω ε re f (x) -ε CMCM (x) : C(x) : ε re f (x) -ε CMCM (x) dΩ Ω ε re f (x) : C(x) : ε re f (x)dΩ (3.49)
where ε re f denotes the strain field obtained by the reference solution, and ε CMCM denotes the corresponding strain field obtained from CMCM. Alternatively, a global L 2 -error is defined as:

Err L2 = Ω u re f (x) -u CMCM (x) . u re f (x) -u CMCM (x) dΩ Ω u re f (x) . u re f (x) dΩ (3.50)
where u re f denotes the displacement field obtained by the reference solution, and u CMCM denotes the corresponding displacement field obtained from CMCM. The reconstruction of the local displacement field can be found in Appendix A.

2D square composite structure with periodic circular inclusions

I consider the academic structure depicted in Fig. 3.6 (a), containing 4 periodic cylindrical fibers. The structure is subdivided into 4 subdomains (see Fig. 3.9 (b)). The size of the square structure is L = 180 mm, and the diameter of the fibers is D = 54 mm. The matrix is assumed to be isotropic elastic with coefficients E mat = 1 Mpa and ν mat = 0.25. The inclusions are also isotropic, with Poisson's ration ν inc = 0.25. I analyze several values of the inclusion's Young modulus: E inc = 5 MPa, E inc = 10 3 MPa (large ratio) and an extreme case E inc = 10 6 MPa to analyze the capabilities of the method in that case. I will especially monitor the case when the subdomains cross the interfaces, as is known to be a difficult one for iterative FETI domain decomposition methods (see [START_REF] Spillane | Solving generalized eigenvalue problems on the interfaces to build a robust two level FETI method[END_REF]).

The structure is subjected to the following Dirichlet boundary conditions: Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures I first solve the problem with the present method without the strain gradient enhancement (CMCM). In order to investigate the convergence of the proposed method, 2 cases will be considered in this example:

u(x) = u x = x 2 y/α u y = -x 3 /3α ∀x ∈ ∂Ω (3.
• For a fixed number of subdomains, I vary the number of coarse mesh elements.

• For a fixed number of coarse mesh elements, I vary the number of subdomains.

Convergence with respect to the number of coarse mesh elements

In this first case, I set the number of subdomains to 4 and increase the number of coarse mesh elements to analyze the convergence of the error. The case E inc = 10 6 MPa is first chosen (most defavorable case). In this work, I use 4-node bilinear elements and 8-node trilinear elements in 2D and 3D respectively for the coarse mesh. Five coarse meshes are used, whose numbers of elements and dofs are provided in Table 3.1. Fig. 3.7 shows a comparison of the micro strain field ε 11 (x) obtained with the CMCM method with the 5 meshes and the reference solution. We can note that the local solution, reconstructed from (5.18), is significantly improved when the coarse mesh is refined.

In Fig. 3.8, the global energy and L 2 errors are plotted as a function of the number of coarse mesh dofs, and for several values of β, the parameter which defines the size of the oversampling (see section 3.2.1). We can see from Fig. 3.8 that the convergence is observed with both respect to the number of coarse mesh dofs and to β.

Influence of the number of subdomains

In this case, the number of coarse mesh elements is fixed and equal to 502, corresponding to 1305 dofs. The number of subdomains is then increased as shown in Fig. 3.9. In the case of 16 subdomains, we can note that the interfaces cross the inclusions. In this test, I also investigate the influence of β on the relocalized solution for the case with 16 subdomains (Fig. 3.9 (c)). Several values of β (β = 0.16, 0.33, 0.47, 0.82 and 1) are used. Figs. 3.11 and 3.12 show the stress fields with respect to β in the cases E inc = 10 6 MPa and E inc = 5 MPa. The evolution of global energy error and L 2 error versus β is given in Fig. 3.10. We can see that the obtained global errors are very high when the interfaces cross the inclusions (i.e. β = 0.16, 0.33 and 0.47). However, when the extended subdomains cover the whole inclusions (i.e. β = 0.82 and 1), the global errors are significantly reduced. Error also decreases with respect to the phase properties contrast. More interestingly, for contrasts of the order of 5, the error is low even for β = 0.

I now use the converged solution in the case of 16 subdomains (i.e. β = 1) and compare to the cases of 1 and 4 subdomains. global error increases. This is due to the fact that some errors are localized at the interfaces between subdomains, thus more subdomains we have, more errors are likely to Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures be committed. However, it is worth noting that, in the third case where interfaces cross inclusions, the global error is 0.032 which is still reasonable considering a high contrast between inclusion and matrix properties (10 6 ). In conclusion, the numerical results show that refining the coarse mesh indeed reduces the global error while increasing the number of subdomains increases the global error, due to higher chances that the subdomains cross the interfaces. When inclusions cut the interfaces, we have seen that the use of larger subdomains (β > 0) (see section 3.2.1) drastically reduces errors. If this case can be avoided, e.g. by means of appropriate choices of subdomains not crossing the interfaces, the errors can also be reduced. These errors also reduce with the contrast of phase properties between the inclusions and the matrix. It is also worth noting that defining an error indicator in this framework and an associated mesh adaptation would be useful and might deserve investigations in future works. Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures three points with coordinates (x c i , y c i ) from left to right as follows: (0.5, 0), (10.5, 1) and (20.5, 0) and are provided as:

p i (x) = p 0 i 1 - x -x c i x p 2 (3.52)
where p 0 2 = 20 MPa, p 0 1 = p 0 3 = 10 MPa and x p = 1 mm. To remove rigid body motions, the nodal y-displacements at nodes (0, 0) and (L, 0), and the nodal x-displacement at node (0, H) are blocked. Three different coarse meshes are considered (Fig. 3.17). The number of elements and the corresponding number of dofs are listed in Table 3 Here, the parameter β (see section 3.2.1) was chosen to be zero. Different components of the local strain field ε(x) are compared with the reference solution, the CMCM and the Second-order CMCM in Figs. 3.18,3.19 and 3.20. To better quantify the accuracy of each solution, the strain and stress 11 components are plotted on the line (x = L/2) in Fig. 3.21. In these figures, only the most converged solutions, on coarse mesh 3, are displayed. In Fig. 3.22, the logarithmic error between the reference solution on a fine mesh and the reconstructed solution using the CMCM method is shown. We can note that in this case, the Second-order CMCM gives a good convergence to the reference solution when refining the coarse mesh, while the CMCM does not converge.

In this example, the method with strain gradient leads to a more accurate solution at almost the same computational costs than the first-order method, as we recall that no additional degrees of freedom are introduced in the macro problem for that purpose. In this example, I consider a square structure as depicted in Fig. 3.23a. It is composed of 30 circular, randomly distributed inclusions. The size of the structure is L = 1 mm. The diameters of the inclusions are also random, ranging from 0.002 to 0.02 mm. The structure is decomposed into 9 subdomains, as shown in Fig. 3.24a. To define these subdomains, the coordinates of the elements centers are tested as belonging to square domains. As the local mesh does not conform to regular square domains, the boundaries of the subdomains are not regular (see Fig. 3.24a). Since the microstructure is not periodic, the subdomain problems are fully computed during off-line calculations. The computational cost for these problems therefore increases, but as all the subdomain problems are independent, they can be performed in parallel, and the increase can in fact be null if a sufficient number of cores is available. In this example, the extended subdomain technique described in section 3.2.2 is applied since the subdomain interfaces cut the inclusions. The structure is subjected to the following boundary conditions:

u(x) = u x = x 2 y/α u y = -x 3 /3α ∀x ∈ ∂Ω (3.53)
with α = 10 -3 and ∂Ω is the boundary of the structure. The comparison of local strain and stress between the reference solution and CMCM solutions using different values of β is plotted in Figs. 3.27 and 3.28, respectively. Fig. 3.29 shows the evolution of global energy errors and online computational time respectively as a function of the number of subdomains, for E i /E m = 10 6 . We can see that the global error slightly increases when more subdomains are used. However, choosing a higher value of β, significantly reduces the errors. The increase of the number of subdomains decreases the computational time for the online (coarse mesh) calculations. It is because that the coarse problem does not change, therefore the online calculations are more parallelizable. From the above results, we can conclude that the present technique can be applied to non-periodic microstructures with non-conforming coarse meshes.

Large scale simulation involving 1.3 Billion of dofs

In this example, I demonstrate the capabilities of the method to handle fully detailed structures with large dimensions, and a larger scale ratio between structure dimensions and microstructural details, as found in some industrial applications. The structure has the same characteristics as the one presented in 2.6.2, but involves 49 × 10 × 10 subdomains (see Fig. 3.31). The material properties and the boundary conditions are also the same. Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

The complete model involves roughly 1.3 billion dofs. In such a situation, most available FEM solvers are not able to perform the simulation even with high performance computing. The coarse mesh used for the CMCM method, however, involves 96,800 trilinear hexahedral elements corresponding only to 321,489 dofs (Figure 3.32), indicating that there is still room for possible refinement of this mesh to improve accuracy. The parameter β (see section 3.2.1) was chosen as zero. To provide a more accurate solution in the vicinity of the load, the coarse mesh is refined in some regions near the pressure fields, as depicted in Fig. 3.32(a). For illustration, I relocalize the local strain and stress fields in two block columns of the structure (see Fig. 3.33). As I have already shown that, for such bending-dominated cases, the Second-order CMCM method improves the accuracy without much additional costs, I have here only used the Second-order CMCM method. The relocalized strain and stress fields on the chosen subdomains are shown in Figs. 3.34,3.35 and 3.36. It can be seen that the CMCM method is able to capture the effects of the weaving patterns under 3 point bending (ε 33 (x) and σ 33 (x)). The results also show that the relocalized fields are less sensitive than those obtained by the filter-based method presented in 2.6.2. From the obtained results, it can be seen that in all three cases, the relocalized strain and stress fields have well captured the effects of three loadings types on the microstructure. In addition, even with a small value of β, I was able to reduce significantly the Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures discontinuities at subdomain interfaces.

Conclusion

In this chapter, I have introduced a new method for parallel computations of large heterogeneous structures, namely Coarse Mesh Condensation Multiscale (CMCM). In this technique, the key idea is to construct an approximation of the full field solution in the structure at the scale of heterogeneities by combining solutions computed on a fine mesh on subdomains decomposing the structure and a solution computed on a coarse mesh. By analogy with the FE 2 method [START_REF] Feyel | A multilevel finite element method (FE 2 ) to describe the response of highly non-linear structures using generalized continua[END_REF], elements of the coarse mesh provide information for applying Dirichlet on the boundaries of subdomains covering the structure. As a result, the full-field solution is approximated over the whole structure at the cost of a finite number of calculations on subdomains which can be solved in parallel and a coarse mesh problem with a low number of dofs. In this chapter, I have only treated the linear case. In this situation, matrix relationships can be obtained to relate the dofs of the coarse mesh and the dofs of the fine mesh subdomains solutions, resulting in a oneiteration scheme only. An extension of the approach to strain gradient has been proposed to enhance the solution in case of global bending of the structures. The accuracy of the method has been tested through benchmark problems involving subdomains crossing the interfaces in the case of high contrast between the phase, which is known to constitute a difficult case for iterative domain decomposition methods. I have also presented an industrial-sized example of a composite beam involving 1.3 billion dofs which has been computed on a standard 32-core standard workstation to show the potential of the method.

The method has also been extended in this work to non-periodic heterogeneous structures which may arise e.g. from micro tomography images of heterogeneous structures. Non-periodic subdomains allow decomposing arbitrary, non-periodic heterogeneous structures. Non-conforming coarse mesh can be used. Then, it is not required that the coarse mesh conforms with the subdomain boundaries. A numerical analysis has been performed on benchmarks. It has been shown that when heterogeneities cross the boundaries of the subdomains, discontinuities can occur in the reconstructed local solutions. Extended subdomain is employed to drastically reduce the errors in that case. I have shown that using a non-periodic coarse mesh allows reducing the errors in the reconstructed solution by local refinement around singularities such as concentrated loads. Finally, I have shown that this method can be used to perform computations on micro tomography models of heterogeneous structures, involving hundred of millions of dofs, on a standard computer.

Chapter 4

Extension of the CMCM method to thermoelastic problems

In this chapter, the CMCM method introduced in the previous chapter is extended to thermoelastic problems. Structures used in aircraft engines can be submitted to high temperatures. Due to the local heterogeneities of thermal coefficients and local thermal gradients, the structure can deform as compared to the room temperature configuration. It is therefore very important to study the thermoelastic behavior of woven composite structures. As the thermoelastic problem is a weak coupling problem between thermics and mechanics, I first extend the CMCM to thermal diffusion problems. Then, I add to the elastic CMCM the ability to account for thermal eigenstrains which can be obtained by the solution of the thermal problem. The different developments and numerical FEM implementation details are presented, as well as benchmarks and application problems to woven composite structures. In this section, I introduce the algorithm of the CMCM method for the thermal diffusion problem. I consider a heterogeneous structure defined in a domain Ω ⊂ R D , D = 2, 3, with boundary ∂Ω = ∂Ω T ∪ ∂Ω q , ∂Ω T ∩ Ω q = ∅, where ∂Ω T and ∂Ω q denote Dirichlet and Neumann boundaries (see Fig. 4.1 (a)). The objective of this problem is Chapter 4. Extension of the CMCM method to thermoelastic problems to determine the temperature at each point of the domain Ω. The different steps of the method are described below.

Coarse Mesh Condensation Multiscale (CMCM) method for thermal problem

Parallel solving of local problems

As in the elastic CMCM formulation, I first decompose the structure into N subdomains Ω α , α = 1, 2, ..., N such that Ω = N α=1 Ω α (see Fig. 4.1 (a)). The boundaries of Ω α are noted ∂Ω α . Then, each subdomain is meshed at the scale of the heterogeneities. This mesh is named in the following the fine mesh. The thermal problem to be solved on each subdomain Ω α reads:

∇ • q(x) = 0 ∀x ∈ Ω α q(x) = -k(x) • ∇T(x) ∀x ∈ Ω α (4.1) (4.2)
Temperature boundary conditions are prescribed over ∂Ω α according to:

T(x) = P(x) h α ∀x ∈ ∂Ω α (4.3)
where h α is a vector of scalar parameters and P(x) is a matrix of functions depending on

x. The boundary conditions are expressed in the form:

T(x) = x 1 x 2 P(x) ∇T α 1 ∇T α 2 h α (4.4)
and in 3D:

T(x) = x 1 x 2 x 3 P(x)    ∇T α 1 ∇T α 2 ∇T α 3    h α (4.5)
Then on each subdomain Ω α , the following local problems are solved: Given h α , find T(x) and ∇T(x) in Ω α such that:

∇ • [-k(x) • ∇T(x)] = 0 ∀x ∈ Ω α T(x) = P(x)h α ∀x ∈ ∂Ω α (4.6) (4.7)
As the problem is linear, we can express the local temperatures and temperature gradients on the fine mesh as:

T(x) = R α (x) h α ∀x ∈ Ω α (4.8) and ∇T(x) = G α (x) h α ∀x ∈ Ω α (4.9)
where R α (x) and G α (x) are the localization matrices relating the local temperature and temperature gradient fields, respectively, to the values in h α . In 2D, these matrices can be expressed as:

R α (x) = T 1 (x), T 2 (x) (4.10) 4.1. Coarse Mesh Condensation Multiscale (CMCM) method for thermal problem and G α (x) = ∇T 1 (x), ∇T 2 (x) (4.11)
and in 3D:

R α (x) = T 1 (x), T 2 (x), T 3 (x) (4.12) and G α (x) = ∇T 1 (x), ∇T 2 (x), ∇T 3 (x) (4.13)
where T i (x) and ∇T i (x) are the temperature and temperature gradient vector solutions, respectively, of the local problem whose boundary conditions are given by h

α i = 1, h α j = 0, j = i.
It can be seen that the number of elementery problems on each subdomain Ω α is equivalent to the length of h α , which is 2 in 2D and 3 in 3D. In total we have 2 × N and 3 × N problems to be solved on 2D and 3D, respectively. It is worth recalling that these local problems are all independent, thus they can be solved in parallel to reduce the computational costs.

Coarse mesh condensation

In the next step, the structure is discretized over the whole domain Ω by a coarse mesh (see Fig. 4.1 (c)). It has been shown in Chapter 3 that the coarse mesh does not necessarily conforms with the boundaries of the subdomains ∂Ω α , i.e. elements of the fine mesh can intersect those of the coarse mesh.

Relation between local temperature gradient and coarse mesh temperature

Given nodal temparatures T E in the coarse mesh, the temperature gradient on the coarse mesh is defined by: ∇T

(x) = B E th (x) T E (4.14)
where B E th (x) is the finite element shape function derivatives for thermal problem of element E defined on the coarse mesh. I define the distance in the sense of the L 2 -norm between the temperature gradient defined over the fine mesh ∇T(x) and the temperature gradient defined over each element of the coarse mesh ∇T(x) as:

L E = Ω E ∇T(x) -∇T(x) 2 dΩ (4.15)
where (v) 2 = v • v. Using (4.9) and (4.14), (4.15) can be rewritten as:

L E = Ω E G α (x) h α -B E th T E 2 dΩ (4.16)
Then the minimization process can be performed independently in each element of the coarse mesh:

h α = argminL E = argmin Ω E G α (x) h α -B E th T E 2 dΩ (4.17)
Chapter 4. Extension of the CMCM method to thermoelastic problems

In each element E of the coarse mesh we obtain the linear system:

Ω E (G α (x)) T G α (x) dΩ H α th h α = Ω E (G α (x)) T B E th (x)dΩ L αE th T E (4.18)
Then we can express h α as a function of the coarse mesh nodal displacements as:

h α = (H α th ) -1 L αE th T E (4.19)
Finally, using (4.9) and (4.19) and a relationship between local temperature gradient and coarse mesh nodal temperatures can be obtained as:

∇T(x) = M α th (x) T E ∀x ∈ Ω α (4.20) with M α th (x) = G α (x) (H α th ) -1 L αE th (4.21)

Relation between local temperature and coarse mesh temperature

We define the distance between the temperatures defined over the fine mesh and each element of the coarse mesh in least-square sense as:

L E T = Ω E T(x) -T E 2 dΩ (4.22)
Using (4.8), (4.22) can be rewritten as:

L E T = Ω E R α (x) h α -T E 2 dΩ (4.23)
Doing the same minimization process as in 4.1.2.1, I obtain a relationship between fine mesh and coarse mesh temperatures according to:

T(x) = Y α (x) T E (4.24)
where

Y α (x) = R α (x) Ω E (R α (x)) T R α (x) dΩ -1 Ω E (R α (x)) T dΩ (4.25)

Solving the problem on the coarse mesh

At this step, we can see that the local temperature gradient (on the fine mesh) can be fully expressed as a function of the coarse mesh nodal temperatures. The objective is then to express the problem to be solved on the coarse mesh. The full-field problem on the fine mesh is given by: 4.2. Coarse Mesh Condensation Multiscale (CMCM) method for thermoelastic problem

           ∇ • q = 0 ∀x ∈ Ω q = -k(x) • ∇T(x) ∀x ∈ Ω T(x) = T d ∀x ∈ ∂Ω T q • n = q n ∀x ∈ ∂Ω q (4.26) (4.27) (4.28) (4.29)
The corresponding weak form is given by: Find T ∈ H 1 (Ω) such that:

Ω k(x) • ∇T(x) • ∇δT(x)dΩ = ∂Ω F q n • δTdΓ ∀δT ∈ H 1 0 (Ω) (4.30)
Introducing (4.20) in (4.30), we obtain:

δT E T ∑ α Ω α (M α th (x)) T k(x) M α th (x) dΩ T E = δT E T ∂Ω F N T th (x) q n dΓ ∀δT E ∈ H 1 0 (Ω) (4.
31) Finally, owing to the arbitrariness of δT E , we obtain the linear system involving only the coarse mesh dofs as:

K th T = f th (4.32)
where:

K th = ∑ α Ω α (M α th (x)) T k(x) M α th (x) dΩ (4.33) f th = ∂Ω t N T th (x) q n dΓ (4.34)
where N th (x) is the thermal finite element shape function matrix of each element of the coarse mesh.

Once the global problem is solved, the local temperature gradient and temperature fields can be reconstructed using (4.20) and (4.24), respectively.

Coarse Mesh Condensation Multiscale (CMCM) method for thermoelastic problem

I now consider the thermoelastic problem in the same domain Ω as in the previous section. It is worth noting that the thermoelastic problem is a one-way coupling problem, which means that the thermal problem has an effect on the elastic problem but not the opposite. In the case where the temperature field over the domain Ω is unknown, the CMCM computation for thermal problem needs to be done prior to the thermoelastic problem.

Parallel solving of local problems

Similar to the CMCM for thermal problem, I consider a heterogeneous domain denoted by Ω whose boundary is denoted ∂Ω with ∂Ω = ∂Ω T ∪ ∂Ω U ∪ ∂Ω q ∪ ∂Ω f , where ∂Ω U and ∂Ω f are partial boundaries where mechanical Dirichlet and Neumann boundary conditions are prescribed, respectively; and ∂Ω T and ∂Ω q are partial boundaries where ). The boundaries of Ω α are denoted by ∂Ω α . In the followings, I consider a thermoelastic problem to be solved on a fine mesh of each subdomain Ω α using finite element discretization.

On each subdomain Ω α , I can define the local problem: where α(x) is the coefficient of thermal expansion, ∆T(x) = T(x) -T 0 is the temperature change and I is the identity matrix.

         ∇ • C(x) : (ε(x) -ε th (x)) = 0 ∀x ∈ Ω α ε(x) = ε e (x) + ε th (x) ∀x ∈ Ω α ε th (x) = α(x) ∆T I ∀x ∈ Ω α
The weak form of the above problem can be written as:

Find u ∈ H 1 (Ω) such that: Ω C(x) : ε(u) : ε(δu)dΩ - Ω C(x) : α(x) ∆T I : ε(δu) dΩ = 0 ∀δu ∈ H 1 0 (Ω) (4.38) or Ω C(x) : ε(u) : ε(δu)dΩ = Ω C(x) : α(x) ∆T(x) I : ε(δu) dΩ ∀δu ∈ H 1 0 (Ω) (4.39)

Local problems based on elastic strain

The boundary conditions of the local elastic problem can be written as:

u(x) = D(x)g α ∀x ∈ ∂Ω α (4.40)
where g α is a vector of scalar parameters and D(x) is a matrix of functions depending on x. For example, in this work, I use the following boundary conditions in 2D:

u(x) = x 1 0 1 2 x 2 0 x 2 1 2 x 1 D(x)   ε α 11 ε α 22 2ε α 12   g α (4.41)
and in 3D:

u(x) =   x 1 0 0 0 1 2 x 3 1 2 x 2 0 x 2 0 1 2 x 3 0 1 2 x 1 0 0 x 3 1 2 x 2 1 2 x 1 0   D(x)         ε α 11 ε α 22 ε α 33 2ε α 23 2ε α 13 2ε α 12         g α (4.42)
which are identical to classical KUBC boundary conditions in homogenization (see e.g. [START_REF] Yvonnet | Computational Homogenization of Heterogeneous Materials with Finite Elements[END_REF]).

It can be seen from the above boundary conditions that there are 3 local problems based on elastic strain in 2D and 6 in 3D. It is worth noting that when the local problems based on elastic strain are solved, the thermal strain is null, therefore the body forces induced from this type of strain are zeros. Consequently, the local problems based on elastic strain can be rewritten as:

Given g α , find T(x) and ∇T(x) in Ω α such that:

∇ • (C(x) : ε(x)) = 0 ∀x ∈ Ω α u(x) = D(x)g α ∀x ∈ ∂Ω α

Local problem with thermal eigenstrains

In addition to the local problems based on elastic strain, one more problem needs to be solved on each subdomain to prescribe the thermal eigenstrains, for both 2D and 3D cases. In order to determine the boundary conditions, the local temperature field is required. I define the two following cases:

• The temperature is known and is constant over the whole domain Ω α .

In this case, the expression of the body force term is straightforward:

f = Ω α C(x) : α(x)∆T I : ε(δu) dΩ (4.45)
Using finite element discretization, the body forces can then be rewritten as:

f = Ω α B T (x) C(x) α(x) ∆T I dΩ (4.46)
• The temperature is not constant over the domain Ω α .

In this case, a thermal problem calculation detailed in Section 4.1 needs to be solved to determine the temperature at each point of the fine mesh. Then the body forces Chapter 4. Extension of the CMCM method to thermoelastic problems can be expressed as:

f = Ω α C(x) : α(x) ∆T(x) I : ε(δu) dΩ (4.47)
Using finite element discretization, the body forces can be rewritten as follows:

f = Ω α B T (x) C(x) α(x) ∆T(x) I dΩ (4.48)

Relation between local fields and solutions of local problems

In this chapter, I assume linear problems. Using the superposition principle, I can express the local strain (on the fine mesh) of each elementary problem as:

[ε(x)] = A α e (x) g α + ε α th ∀x ∈ Ω α (4.49)
where A α e (x) is a localization matrix relating the local strain field in the fine mesh to the values in g α and A α th (x) is the solution of the local problem based on thermal strain. Note that higher order terms can also be included in g α . However, due to time constraints, I have not considered this case in the present study.

Equation (4.49) can be rewritten as:

[ε(x)] = [A α e (x), [ε α th ]] A α total g α 1 g α total (4.50)
The number of local problems in each subdomain is equal to the length of g α plus 1, which is here 4 in 2D and 7 in 3D. Therefore, I have 4 × N and 7 × N problems to be solved in 2D and 3D, respectively. However, each of these problems is computed on a small part of the structure, and can thus be solved efficiently. Further, all problems are independent and can be solved in parallel to reduce computational time. It can be noted that the computational time of local problems scales linearly with the number of available processors. Then, we can deduce that the computational time for solving all subdomains problems is 4 × N/ N p in 2D and 7 × N/ N p in 3D with N p is the number of processors. The solutions of the local problems are stored for the next step. In 2D, the columns of A α e (x) are formed with the solutions of each elementary problems based on elastic strain as:

A

α e (x) = [ε 1 (x)], [ε 2 (x)], [ε 3 (x)] (4.51)
where [ε i (x)] is the strain vector solution of the local problem whose boundary conditions are defined associated with g α i = 1, g α j = 0, j = i. In 3D, A α e (x) has then 6 columns. As in the elastic problem, when a subdomain boundary crosses an inclusion, an extended subdomain technique can be employed to reduced the errors at the boundary (see section 3.2.1).

Coarse mesh condensation

Similar to the thermal problem, the structure is now discretized over the whole domain Ω by a coarse mesh (see Fig. 4.2 (c)). Given nodal displacements u E in the coarse mesh, 4.2. Coarse Mesh Condensation Multiscale (CMCM) method for thermoelastic problem the strain field on an element E of the coarse mesh is defined by:

[ε(x)] = B E (x) u E (4.52)
where B E (x) is the finite element shape function derivative matrix of element E defined on the coarse mesh. For an arbitrary point x in Ω, x belongs both to (i) a coarse mesh element E and (ii) a fine mesh element e lying in a subdomain Ω α . I define the distance in the sense of the L 2 -norm between the strain field defined over the fine mesh [ε(x)] and the strain defined over the coarse mesh as:

J = Ω ([ε(x)] -[ε(x)]) 2 dΩ (4.53)
where (v) 2 = v • v. Using (4.50) and (4.52):

J = Ω A α total (x) g α total -B E (x) u E 2 dΩ (4.54)
The above equation can be rewritten as:

J = ∑ E Ω E A α total (x) g α total -B E (x) u E 2 dΩ (4.55)
Then the minimization process can be performed independently in each subdomain:

g α total = Argmin J E = Argmin Ω E A α total (x) g α total -B E (x) u E 2 dΩ (4.56)
We obtain a linear system in each element E of the coarse mesh of the form:

Ω E (A α total (x)) T A α total (x) dΩ H α g α total = Ω E (A α total (x)) T B E (x) dΩ L αE u e (4.57)
I can express g α as a function of the coarse mesh nodal displacements as:

g α = (H α ) -1 L αE u E (4.58)
Finally, using (4.49) and (4.58), a linear relationship between local strain in the fine mesh and the coarse mesh dofs can be established as:

[ε(x)] = M α (x) u E ∀x ∈ Ω α (4.59) with M α (x) = A α total (x) (H α ) -1 L αE (4.60)

Solving the problem on the coarse mesh

At this step, we can note that the local strain (on the fine mesh) can be fully expressed as a function of the coarse mesh nodal displacements. The objective is then to express the Chapter 4. Extension of the CMCM method to thermoelastic problems problem to be solved on the coarse mesh. The full-field problem on the fine mesh is given by: It can be seen that the problem on the coarse mesh contains only elastic boundary conditions. The reason is that the thermal boundary conditions have been taken into account when solving the thermal problem in order to determine the thermal eigenstrains ε th (x). Using vector forms for second-order tensors and matrix form for fourth-order tensor, the corresponding weak form reads:

               ∇ • C(x) : (ε(x) -ε th (x)) = 0 ∀x ∈ Ω ε th (x) = α(x) ∆T I ∀x ∈ Ω u(x) = u d ∀x ∈ ∂Ω u σ • n = f n ∀x ∈ ∂Ω f
Find u ∈ H 1 (Ω) such that:

Ω [ε(δu)] T C(x) [ε(u)]dΩ = Ω [ε(δu)] T C(x) α(x) ∆T(x) I dΩ + ∂Ω f f n • δu dΓ ∀δu ∈ H 1 0 (Ω) (4.65)
Introducing (4.59) in (4.65), we obtain:

δu E T ∑ α Ω α (M α (x)) T C(x) M α (x) dΩ u E = δu E T ∑ α Ω α (M α (x)) T C(x) α(x) ∆T(x) I dΩ+ δu E T ∂Ω F N T (x) f n dΓ ∀δu ∈ H 1 0 (Ω) (4.66)
Finally, owing to the arbitrariness of δu E , we obtain the linear system involving only the coarse mesh dofs as:

K u = f (4.67)
where:

K = ∑ α Ω α (M α (x)) T C(x) M α (x) dΩ (4.68) f = ∑ α Ω α (M α (x)) T C(x) α(x) ∆T(x) I dΩ + ∂Ω t N T (x)f n dΓ (4.69)
and N(x) are the classical finite element shape function matrix of elements in the coarse mesh. Once the global problem is solved, the local total strain field in each subdomain can be reconstructed using (4.59) and the local elastic stress field is reconstructed using:

[σ(x)] = C(x) M α (x) u E -ε α th ∀x ∈ Ω α (4.70) 4.3. Numerical examples

Numerical examples

The objective of the following numerical examples is to illustrate the efficiency of the present CMCM method in solving thermoelastic problems of heterogeneous materials.

In each example, two solutions are compared:

• A reference solution, denoted by "ref" in the different figures, is obtained by a Direct Numerical Simulation (DNS)-finite element calculation on a fully meshed structure at the scale of the heterogeneities.

• A solution obtained by the CMCM method.

For a better visualization of the results, I define a global energy error Err E and a global L 2 error Err L2 such as:

-For the thermal problem

Err E = Ω (∇T re f (x) -∇T CMCM (x)) : k(x) : (∇T re f -∇T CMCM (x))dΩ Ω ∇T re f (x) : k(x) : ∇T re f (x)dΩ (4.71) Err L2 = Ω (T re f (x) -T CMCM (x)) • (T re f -T CMCM (x))dΩ Ω T re f (x) • T re f (x)dΩ (4.72)
where ∇T re f denotes the temperature gradient field obtained by the reference solution; ∇T CMCM denotes the temperature gradient field obtained from CMCM; T re f is the temperature field obtained from reference solution and T CMCM is the reconstructed temperature field obtained from CMCM.

-For the thermoelastic problem, the errors are defined by Eqs. (3.49) and (3.50).

Tension test of 2D structure

In this example I consider a thermoelastic problem whose structure and boundary condition are shown in Fig. 4.3a. The size of the square structure is L = 3 mm, the diameter of each circular inclusion is 0.54 mm. The fine mesh used to describe the structure including all inclusions contains 21,986 linear triangular elements, corresponding to 22,388 degrees of freedom (Fig. 4.3b). The coarse mesh used to describe the structure at the macro scale is composed of 162 quadratic triangular elements, corresponding to 722 degrees of freedom (Fig. 4.3c). The mechanical and thermal properties of each component of the structure are listed in Table 4.1. The structure is decomposed into 9 periodic subdomains. Since the temperature is unknown and is not constant over the whole structure, a CMCM calculation for thermal problem presented in Section 4.1 needs to be done in order to determine the temperature at each Gauss point of the fine mesh.

Chapter 4. Extension of the CMCM method to thermoelastic problems Matrix Inclusion Young's modulus (MPa) 1 10 6 Poisson coefficient 0.3 0.3 Thermal conductivity (W.m -1 .K -1 ) 1 50 Coeff of thermal expansion (K -1 ) 5 × 10 -5 10 -5 The CMCM solutions for thermal problem with β = 0 and 1 are presented in Fig. 4.4. The error curve of the L 2 error in function of β for the thermal problem is shown in Fig. 4.6b. From these figures, we can note that both reference and CMCM re-localized solutions are very similar, for both β = 0 and 1. Indeed, the values of L 2 error for β = and 1 are quite low, 1.7 × 10 -4 and 9.9 × 10 -5 , respectively. Therefore if one only wants to determine the local temperature field for the use of the thermoelastic problem, β = can be applied in order to reduce computational time.

The comparison between temperature gradient fields of reference solution and relocalized CMCM solutions are presented in Fig. 4.5. The curve of the global errors as a function of β is plotted in Fig. 4.6. It can be seen from these figures that in this case, β = is not sufficient to obtain a good accuracy, and β = 1 is needed to obtain a significant error reduction, especially at the interfaces between subdomains. It can be seen from these figures that increasing the value of β does not decrease significantly the global error. One possible explanation is the presence of high temperature gradients, which are not well taken into account into the solution approximation. One future improvement of this scheme could be to include temperature gradients into the approximation of local basis solutions. Note that I also did not include here the strain gradients as in the previous chapter, due to time constraints. I then believe that there is room for improvement of the CMCM solutions in the present thermoelastic framework.

Conclusion

In this chapter, the Coarse Mesh Condensation Multiscale (CMCM) method, introduced in Chapter 3 has been extended to thermal diffusion and thermoelastic problems of heterogeneous materials. As in the elastic case, the micro temperature solutions are computed in subdomains by solving a series of local calculations which can be performed in parallel. The introduction of thermoelastic strains in the mechanical problems can then be defined in a second step, by introducing the thermal eigenstrains in the mechanical problem through additional calculations on the mechanical subdomains problems. Therefore, the whole procedure can be maintained as fully parallel. Numerical benchmarks in 2D have been investigated to quantify the errors introduced by the CMCM approximation in thermoelastic cases. The results show good accuracy for local temperature, strain and temperature gradient fields, even though an extended subdomain technique is required to remove local errors at the subdomain interfaces. I have however noticed a poor convergence of the global error. This might be explained by the fact that I did not include strain and temperature gradients terms in the definition of the boundary conditions of the local problems. This point constitutes a possible improvement of this extension and is reported in the perspectives of this work.

Chapter 5

Extension of the CMCM method to nonlinear problems

In this chapter, I investigate the extension the CMCM method to nonlinear problems. One strong difficulty inherent to nonlinear problems is that the superposition principle, which is used as a foundation in the CMCM to decompose the local solution into subdomain basis solutions, cannot be applied. To circumvent this difficulty, the method is applied to the linearized problem which arises in a classical Newton solution procedure. At each iteration of the Newton algorithm, the same procedure as in the thermoelastic case can be applied, and the parallel solving procedure still holds, but the local solutions on the subdomains must be re-computed at each iteration of the Newton algorithm. It can be observed that, in a nonlinear context, the CMCM method partly resembles to the FE 2 method introduced by [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF]. The analogies and differences between these two methods are also discussed in this chapter. I provide tentative numerical examples on a elastic-viscoplastic composite structure to test the capability of the method to handle a complex nonlinear behavior and the resulting accuracy.

Reference problem on the structure fine mesh

First, I express the local balance equations of the nonlinear problem as:

∇ • σ(u(x)) + f = 0 (5.1) with σ(x, t) = ξ(ε(x, t)) (5.2)
where ξ is a nonlinear operator describing the nonlinear behavior of the considered material.

The weak from of the above problem is given by:

Find u ∈ H 1 (Ω) such that: Ω σ(u(x)) : ε(δu(x))dΩ = Ω f • δudΩ + ∂Ω F • δudΓ ∀δu ∈ H 1 0 (Ω) (5.3)
which can be rewritten as:

Ω σ(u(x)) : ε(δu(x))dΩ - Ω f • δudΩ - ∂Ω F • δudΓ = R(u(x)) ∀δu ∈ H 1 0 (Ω) (5.4) with R(u(x)) = 0 (5.5)
In order to solve the above problem in a nonlinear context, I employ the iterative Newton method. Using Taylor's expansion at first order, I express the residual R as:

R(u k + ∆u) = R(u k ) + D ∆u R(u k ) + o(∆u.∆u) (5.6)
where u k is a known solution at the k th iteration, ∆u is the difference between the solution at k th iteration and the previous one and D ∆u R(u k ) is the directional derivative of R(u k ) with respect to ∆u, defined by

D v f (u) = d dα { f (u + αv)} α=0
(5.7)

Combining (5.7) and (5.6) and neglecting high order terms o(∆u.∆u) we obtain:

D ∆u R(u k ) = -R(u k ) (5.8)
Substituting (5.4) in (5.8) we have:

D ∆u R(u k ) = D ∆u   Ω σ(u k (x)) : ε(δu(x))dΩ   = Ω D ∆u σ(u k (x)) : ε(δu(x))dΩ = Ω ∂σ ∂ε u k : ε(∆u) : ε(δu)dΩ = Ω C tan u k : ε(∆u) : ε(δu)dΩ
(5.9)

The weak form of the linearized problem is then expressed as:

Find u ∈ H 1 (Ω) such that: Ω C tan u k : ε(∆u) : ε(δu)dΩ = - Ω σ(u(x)) : ε(δu(x))dΩ + Ω f • δudΩ + ∂Ω F • δudΓ ∀δu ∈ H 1 0 (Ω) (5.10)
This linearized problem will be solved by the elastic CMCM formulation described in Chapter 3. Solving the above problem gives a correction to be applied on the known displacement field. The algorithm stops when an error criterion is reached. Here, I adopt the norm of the residual on the right hand side of (5.10) to be lower than a prescribed tolerance. A parallel solving procedure of (5.10), similar than in the thermoelastic formulation of Chapter 4, is described in the following.

Local problem on subdomains

Local problem on subdomains

As in the elastic CMCM formulation, a series of local problems needs to be solved to obtain the behavior information of the structure at the fine scale at each step of the Newton iteration. For this purpose, the structure is again decomposed in several subdomains (see section 3.1). It is worth noting that the linearized problems are to be solved on each subdomain, and since they are fully independent, parallel computing can be applied to reduce computational time. The only difference with the linear problem is that this operation must be repeated at each iteration of the Newton method. The boundary conditions remains in the same form than for the linear problem, and are expressed by (3.2).

The local problem on a subdomain α can then be expressed as follows.

Given g α = {ε, ∇ε}, find ε(x) in Ω α such that:

         ∇ • σ α (u(x)) + f = 0 ∀x ∈ Ω α σ α (x, t) = ξ(ε α (x, t)) ∀x ∈ Ω α u(x) = εx + 1 2 G : x ⊗ x ∀x ∈ ∂Ω α
(5.11)

(5.12)

(5.13)

Applying the iterative Newton method presented in the previous section, I obtain the linearized problem to be solved at iteration k as:

Find u k ∈ H 1 (Ω α ) such that: Ω C tan u k-1 : ε(∆u k ) : ε(δu k )dΩ = - Ω σ(u k-1 (x)) : ε(δu k (x))dΩ + Ω f • δu k dΩ ∀δu k ∈ H 1 0 (Ω α ) (5.14) with u k (x) = ε • x + 1 2 G : x ⊗ x ∀x ∈ ∂Ω α (5.15)
It can be noted that the resolution of the above problem requires the knowledge of C tan , which is calculated based on the micro strain field reconstructed from the previous iteration.

Global problem on a coarse mesh

I employ the same methodology as in the linear thermoelastic problem but applied to the linearized problem which arises in the Newton solving procedure. It is worth noting that the procedure is very similar to the thermoelastic problem of Chapter 4, as terms in the residual play the role of thermal eigenstrains.

The micro strain field of the micro problem can be expressed as:

ε(∆u) = ∆ε = A α (x)H α (5.16)
where A α (x) is the localization matrix obtained by solving localization problems on each subdomain. In order to perform a condensation procedure, I minimize the distance between the micro and macro strain fields in each macro element of the coarse mesh in the 5.5. Analogies between FE 2 and CMCM 1: Initialize u = 0 2: Initialize A α (x) = 0 3: Initialize tol = 10 -6 4: WHILE Err > tol 5:

LOOP over all subdomains Ω α (parallel) 6:

-Reconstruct micro strain field using (3.34) 7:

-Compute C tan (x) and store 8:

-Solve localization problems to obtain A α (x) and store 9:

-Compute H α using (3.32) and store 10:

END LOOP

11:

LOOP over all macro elements E 12:

-Find micro elements e ∈ Ω E 13:

-Find the subdomain α contaning E 14:

-Compute L αE using (3.33) 15:

-Compute T α (x) using (3.35) and store 16:

-Given T α (x):

17:

LOOP over all micro elements e ∈ Ω E 18:

-Compute C tan 19:

-Compute K

E tan = K E tan + Ω E (T α (x)) T C tan (x) T α (x) dΩ 20: -Compute R E = R E -Ω E (T α (x)) T σ(u k ) dΩ + Ω E N T (x) f dΩ + ∂Ω F N T (x) F n dΓ 21: END LOOP 22: -Assemble K E tan in K tan 23: -Assemble R E in R 24:
END LOOP 25:

-Solve K tan ∆u = R to obtain ∆u 26:

-Update u = u + ∆u 27:

-Evaluate Err = norm(R) 28: END WHILE Algorithm 3: CMCM algorithm in a nonlinear context

Analogies between FE 2 and CMCM

In this subsection, I show analogies between the present method and the FE 2 method, originally proposed in [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF] to solve nonlinear multiscale problems with separated scales. In FE 2 (see Fig. 5.1 (b)), a relationship between macro strain ε and the macro stress σ is obtained numerically at each Gauss point of the coarse mesh by: (i) prescribing Dirichlet boundary conditions depending on ε on the boundary of an RVE associated with the element; (ii) solving the local problem on the RVE; (iii) averaging the stress in the RVE to obtain σ. In CMCM (see Fig. 5.1 (a)), a group of elements is covering a subdomain. Solving a minimization problem, Dirichlet boundary conditions are prescribed on the boundary of the subdomain. The solution is then solved on the subdomain to obtain the stress field.

The main differences are summarized as follows:

CMCM FIGURE 5.1: (a) Schematic description of FE 2 method Feyel (1999) and of CMCM

• In FE 2 , the local calculations are computed on an RVE associated with the microstructure. In CMCM, the local calculations are performed in subdomains whose union is the fully detailed heterogeneous structure

• In FE 2 , the parameters defining the boundary conditions on the local problem (RVE) are associated with a single integration point in the mesh at the macroscale. In CMCM, the parameters associated with boundary conditions to be applied to the local problem (subdomain) depend on the strain field of all elements of the coarse mesh covering the subdomain • In FE 2 , the local problem is used to define the macro strain -macro stress relationship. In CMCM, the local problem is used to relate the strain in the coarse mesh to the full stress field in all elements in the subdomain (fine scale) covered by the coarse mesh elements and the notion of macro stress is no more necessary Both methods share the following features:

• All local computations can be performed in parallel

• The local solution (fine scale) can be reconstructed given the macro (coarse) mesh solution

• The macro (coarse) mesh problem involves much fewer degrees of freedom as compared to solving the full problem However it is necessary to make it clear that CMCM is not a special case of FE 2 . Also, note that, in CMCM, the full microstructure has to be known everywhere in the structure, in contrast with FE 2 . However, we recall that this method applies to structures where scales cannot be separated, and where the size of heterogeneities becomes not negligible as compared to the size of the structure, as for example in woven composites. In industrial applications, we are going toward more realistic structures, either by X-ray tomography or process-induced simulation. In this section, I evaluate the accuracy of the CMCM method in a nonlinear context. for this purpose, I consider a 2D beam made of elastic-viscoplastic constituents as described in Fig. 5.2. A more detailed description of the material parameters used in this example can be found in [START_REF] Bower | Applied mechanics of solids[END_REF].

The strain tensor is classically decomposed into elastic and plastic parts as:

ε(x) = ε e (x) + ε p (x) (5.23)
where ε e (x) and ε p (x) are elastic and plastic strain, respectively. The elastic strain is related to stress using the Hooke's law of linear elasticity:

σ(x) = C(x) : ε e (x) (5.24)
In this example, the elastic domain is determined using von Mises criterion defined as:

f = σ eq = 3 2 σ d (x) : σ d (x) ≤ 0 (5.25)
where σ d (x) = σ(x) -1 3 Tr(σ(x))I is the deviatoric stress tensor where I is the second order indentity tensor.

The creep state within each phase of the structure is modeled using von Mises flow potential with power-law rate sensitivity:

ψ(σ eq , σ 0 ) = ε0 σ eq σ 0 m (5.26)
The plastic strain rate is given by:

εp = ε0 σ eq σ 0 m 3 2 σ d σ eq (5.27)
where σ 0 is the flow stress and n is material constant. The isotropic hardening rule, often used to model transient creep, which specifies the evolution of state variables with plastic straining is defined as follows [START_REF] Clifton | High strain rate behavior of metals[END_REF]:

σ 0 = Y 1 + ε p ε 0 1/n
(5.28)

Chapter 5. Extension of the CMCM method to nonlinear problems where ε p e = ε0 σ eq σ 0 m dt is the cumulative plastic strain; ε 0 and ε0 are reference strain and strain rate, respectively; Y, m and n are material constants.

The material properties of the matrix and the inclusions are summarized in The curve describing the relation between applied force and maximum displacement of the reference and the CMCM solution over 10 loading steps is plotted in Fig. 5.3. The strain-stress evolution over 10 loading steps is also illustrated in Figs. 5.4 and 5.5 for 11and 22-components of the stress. We can note that the force-displacement relationship is well reproduced by the CMCM in this nonlinear case. Regarding stress response, the 11-stress response gives reasonable agreement while there are noticeable errors for the 22-component. I examine in more details the local fields in Figs. 5.6, 5.7 and Fig. 5.8, where local strain and cumulated plastic strain are depicted as a function of β. We can note that errors are non-negligible, especially near the boundary conditions. This can be explained by the fact that in the present chapter, I did not add local solutions related to strain gradients as in the elastic case, due to time constraints. I then believe that this point could be improved in a near future by adding these preliminary calculations in the off-line procedure.

Conclusion

Conclusion

In this chapter, another contribution is proposed, by extending the CMCM to nonlinear problems. The main idea is to apply the CMCM to the linearized problem which arises in a classical Newton method solution procedure. The different steps of the method are maintained, but solutions in the subdomains must be computed at each Newton iteration.

A first test on an elastic-viscoplastic composite structure has proved the applicability of the method and provided encouraging results. I have reported local errors which do not vanish with an increase of the β parameter. As in the thermoelastic extension, I believe that these errors are due to the fact that I did not include strain gradient terms in solving local subdomain problems. Due to time constraints, more in-depth investigations of this point is reported to the perspectives.

Conclusion and Perspectives Conclusion

In this Ph.D., I have developed multiscale approaches using finite elements to solve composite structures without scale separation and involving structure gradients. These features are especially important in industrial applications such as woven composites. More specifically, I have proposed and developed an original and new numerical method, called the Coarse Mesh Condensation Multiscale (CMCM) method, which allows taking into account the microstructural kinematics in heterogeneous structures, while maintaining advantages of multiscale methods such as parallel calculations and use of a coarse mesh to solve the full problem with a classical finite element approximation. The different contributions of this Ph.D. are summarized below.

In the beginning of the Ph.D. period, I have first investigated an existing method, called filter-based computational homogenization [START_REF] Tognevi | A multi-scale modeling method for heterogeneous structures without scale separation using filter-based homogenization scheme[END_REF]. This method has the advantage not to depend on the scale separation assumption and allows reconstruction of local fields even when strain gradients are involved. I have also extended the method to 3D cases using unstructured meshes at both micro and macro scales. As compared to previous studies, I have provided a more in-depth analysis of errors in problems with microstructure gradients and strong global strain gradients in elastic problems. I have shown that the method allows significantly reducing errors, particularly on macroscopic quantities of interest, but also on local fields as compared to classical first-order homogenization method with enhanced reconstruction techniques. However, one drawback lies on the computational costs, as the local fields must be computed using an iterative scheme even in the linear case. The second drawback is that even though reduced, the errors on relocalized quantities remain non-negligible and I have found a sensitivity related to the macroscopic mesh, especially in an example of woven structure.

For this reason, I have developed a new technique in the second part of this study, which constitutes the main contribution of this Ph.D.. The Coarse Mesh Condensation (CMCM) method has been introduced Le et al. (2020a,b) to efficiently solve large heterogeneous structures including all microstructural details, while avoiding the issues of scale separation, characteristic of homogenization-based multiscale methods. The main ideas of the method are as follows: (i) the structure is fully meshed a the level of microstructural details; (ii) the mesh is decomposed into subdomains which cover the entire structure; (iii) local problems are solved on the subdomains with appropriate boundary conditions, to construct a basis of local solutions; these problems are independent and can be solved in parallel; (iv) the local solutions are related to the degrees of freedom of a coarse mesh covering the structure by a condensation operator derived from a non-iterative least-square minimization principle; (v) the problem is solved on the coarse mesh, which features a much smaller number of degrees of freedom; (vi) local solutions can be reconstructed efficiently to obtain the microstructural field details in all the structure. An in-depth error analysis has been conducted, showing that the use of two ingredients is important for the accuracy of the method: (i) the use of "extended subdomains" Chapter 5. Extension of the CMCM method to nonlinear problems in the local problem solving step to avoid discontinuity errors across non-periodic subdomains and (ii) the use of strain gradient-type boundary conditions in the same step to capture global strain gradients within the structure. I have demonstrated the capability of the technique by solving a variety of increasingly complex examples.

Within the developments of this method, I have proposed several other contributions. First, I have extended the method to non-periodic microstructures, which can arise e.g. in the analysis of micro tomography-images models. For this purpose, I have introduced a technique to handle coarse meshes which do not conform to the local subdomains with an appropriate integration technique, relying on the subdivision of fine mesh elements. Another advantage of this improvement is that the coarse mesh can be locally refined to capture localized fields. The method has been demonstrated by solving a problem arising from an X-ray image of a concrete microstructure specimen.

Second, I have applied the method to thermoelastic problems, for the analysis of thermally-induced strains in aircraft engines structures. The CMCM method has been extended, by first providing both formulations and numerical implementations for thermal diffusion problems, and then to elastic problems with thermal eigenstrains.

Finally, I have investigated the extension of the method to nonlinear problems. One difficulty was that the superposition principle, which is a basis in the CMCM method to express the subdomain bases solutions, could no longer be applied. To circumvent this difficulty, I have proposed a methodology where the main steps of the method are applied to the linearized problem in a classical iterative Newton solving procedure. The complexity increases, as the parallel calculations on the subdomains must be re-computed at each iteration. Even though I have only conducted preliminary studies, encouraging results have been provided. For both thermoelasticity and nonlinear problems, only the first-order (barring strain-or temperature-gradient terms) version of the method has been tested.

In this Ph.D., all codes have been developed and implemented in a in-house Python code, with the aim of connecting it to commercial pre-and post-processors and local solvers. Technical developments have been performed, such as preliminary assessments of performance in parallel computing, which are at the core of this method. The calculations have been implemented using different parallel algorithms such as Divide-andconquer or Master/slave. The Divide-and-conquer paradigm has the possibility to fully exploit the computing resources but requires large computer memory. The Master/slave paradigm was found to maintain approximately the same performances while using much less memory and appears to be the better choice.

Perspectives

This Ph.D. opens many perspectives described as follows.

The first possible progresses include improvement of the error reduction with the CMCM method. I have shown in the linear case that significant reduction of local errors could be obtained by an extended subdomain technique. However, this method implies larger subdomains, which may be very costly in 3D. One enhancement of the method would include a correction treatment at the interface to avoid the extended subdomain technique.

Other improvements include the error reduction in thermal and nonlinear problems. For these problems, I did not include, for time reasons, higher-order terms in the subdomain boundary conditions. Regarding nonlinear structures, it is worth noting that there 5.7. Conclusion are no special restrictions regarding the local behavior of phases. Therefore, other nonlinearities could be considered, such as finite displacements and strains, viscoplasticity, damage, among others, which may arise in service-use of composite structures.

Another exciting possibility is to apply the present method to crack propagation. The CMCM method allows the analysis of micro fields in cracked structure without introducing the crack in the fine mesh (see e.g. Fig. 5.9). Combining the ideas developed in solving nonlinear problems and the possibility of refining the coarse mesh, investigations of micro and macro cracks initiation in complex heterogeneous structures with CMCM could open the route to critical engineering applications. Restricting to linear analysis, which is the most advanced area for this method so far, the CMCM method could be applied to more complex industrial composite structures, for example a complete part like a fan blade of the aircraft engine, where strong microstructure gradients are observed. These problems are mainly solved in a thermoelastic context, due to the intense thermal cycles of in-service aircraft engines. Some mesoscale defects could also be introduced in such structures to take into account processing anomalies or damage tolerance.

The CMCM method has been implemented in a Python code allowing to solve large composite structures in parallel. However, it has not been fully optimized. One option to improve the code is to optimize the local solvers. The idea is to use the linear and nonlinear solvers developed in the Z-set suite. In fact, since these solvers are already optimized and implemented in C++, the speedup of the method can be enormous. Another option is to optimize the global solver, especially in nonlinear problems, where the calculations of certain operators are repeated in the current version.

B.2. Parallel computing paradigms

1: data = comm.bcast(data, root=0) Broadcast the initial data to all processors 2: Decompose fine mesh in subdomains and divide corresponding data 3: Get number of processors using "comm.Get_size()" 4: Do the work and save the result 5: LOOP over all processors apart from 0

Send results from all other processors to 6:

processor 0 for post-processing

B.2.2 Master/slave paradigm

The Master/slave paradigm is also a very efficient approach in parallel computing. Considering a set of jobs that needs to be executed by a set of processors. The latter are divided into one master and the rest as slaves. In the Master/slave paradigm, all the so-called light works such as pre and post-processing of data, are set to be done by the master; and all the heavy works such as problem solving are set to be done by the slaves. In other words, the master's job is to pre-process data and then distribute one job at a time to each slave; when the slave finishes its work, it sends back the result to the master for post-processing. This process continues until all required jobs are done. The difference between the Master/slave and the Divide-and-conquer paradigms is that, in the former, there is one processor dedicated to data processing (i.e. the master). Therefore, the initial data will be fully stored only in the master, not in all processors as with the Divide-and-conquer paradigm. However, the Master/slave paradigm is more difficult to implement and communications between processors are more complex.

Sample implementations of the Master/slave paradigm for the master and the slaves are presented in Algorithm 5 and 6, respectively. 1: Decompose fine mesh in subdomains and divide corresponding data 2: Get number of processors using "comm.Get_size()" 3: LOOP over all N slaves Seed the slaves, send one unit of work to each slave 4:

Distribute N first problems to N slave processors 5: END LOOP 6: WHILE true DO

  ) . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Fine and coarse mesh of a RVE . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Illustration of a field in 2D: (a) microscopic; (b) filtered on a coarse mesh . 2.3 Illustration of filtered fields in 1D with different size of the macro mesh . . 2.4 Illustration of: (a) a macroscopic structure; (b) a RVE . . . . . . . . . . . . . 2.5 Positioning of a RVE on the macro structure . . . . . . . . . . . . . . . . . . 2.6 (a) Macro mesh of the structure; (b) Macro strain field ε 11 ; (c) macro strain field ε 11 interpolated on a RVE . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 (a) Structure of the 2D composite beam ; (b) RVE . . . . . . . . . . . . . . . 2.8 (a) Macro mesh of the structure ; (b) macro mesh of the RVE . . . . . . . . 2.9 Comparison of the displacement field along y of: (a) reference solution on a fine mesh; (b) first-order homogenization on a coarse mesh and (c) filterbased homogenization on a coarse mesh . . . . . . . . . . . . . . . . . . . . 2.10 Micro strain field ε 12 : (a) reference solution; (b) filter-based homogenization; (c) first-order homogenization; (d) first-order homogenization with enhanced relocalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Considered woven composite coupon and boundary conditions; the red cylinders are where the pressure fields are prescribed . . . . . . . . . . . . 2.12 Visualisation of the RVE and its mesh . . . . . . . . . . . . . . . . . . . . . 2.13 Distribution of pressure applied to the center of the specimen . . . . . . . 2.14 Strain field on the skin of the reference solution (a), of the macroscopic solution obtained by filter-based method (b) . . . . . . . . . . . . . . . . . ix 2.15 Comparison of the local stress fields σ 22 between the reference solution (a), first-order homogenization (b) and filter-based homogenization (c) on a small section at the center of the structure . . . . . . . . . . . . . . . . . . . 3.1 Schematic description of the CMCM method; (a) decomposition of the structure into subdomains; (b) discretization of the structure with a coarse mesh; (c) elementary (off-line) problems to be solved over the subdomains 3.2 Illustration of an extended subdomain . . . . . . . . . . . . . . . . . . . . . 3.3 Subdomain Ω α , coarse mesh element E associated with a domain Ω E and fine mesh element e associated with a domain Ω e . . . . . . . . . . . . . . . 3.4 Subdomain Ω α cut by a coarse mesh element and different integration schemes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Remeshing operation in the case where micro elements are cut by a macro element . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 (a) 2D square composite structure: geometry; (b) subdomain used for offline calculations; (c) subdomain finite element mesh . . . . . . . . . . . . . 3.7 ε 11 (x) of reference solution and obtained with the CMCM method for different coarse meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Evolution of global energy (a) and L 2 (b) errors with respect to the number of coarse mesh dofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Structure divided into (a) 1 subdomain, (b) 4 subdomains and (c) 16 subdomains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Evolution of energy error (left) and L 2 error (right) in function of different values of β in the case of 16 subdomains with contrast E inc /E mat = 10 6 , 10 3 and 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 σ 11 (x) using different values of β in the case of 16 subdomains with E inc /E mat = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 σ 11 (x) using different values of β in the case of 16 subdomains with E inc /E mat = 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Relocalized stress σ 11 (x) (MPa) E inc /E mat = 10 6 : (a) reference solution, (b) proposed methos with 1 subdomain, (c) proposed method with 4 subdomains and (d) proposed method with 16 subdomains . . . . . . . . . . . . 3.14 Relocalized stress σ 11 (x)(MPa) with E inc /E mat = 10 3 : (a) reference solution, (b) proposed methos with 1 subdomain, (c) proposed method with 4 subdomains and (d) proposed method with 16 subdomains . . . . . . . . . 3.15 Evolution of the global energy error in function of the number of subdomains 3.16 (a) 2D three-point bending composite beam: geometry and boundary conditions (b) subdomain finite element mesh . . . . . . . . . . . . . . . . . . . 3.17 Coarse meshes for global calculation of the beam: (a) coarse mesh 1, (b) coarse mesh 2, and (c) coarse mesh 3 . . . . . . . . . . . . . . . . . . . . . . 3.18 Relocalized strain solution ε 11 (x) . . . . . . . . . . . . . . . . . . . . . . . . 3.19 Relocalized strain solution ε 22 (x) . . . . . . . . . . . . . . . . . . . . . . . . 3.20 Relocalized strain solution ε 12 (x) . . . . . . . . . . . . . . . . . . . . . . . . 3.21 Comparison of reference, CMCM, and Second-order CMCM solutions of local (1, 1) components long the line (x = L/2) . . . . . . . . . . . . . . . . 3.22 Minimum displacement (left) and maximum stress (right) using different coarse meshes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.23 Random structure geometry (a) and its finite element mesh (b) . . . . . . . x 3.24 Reference structure decomposed into 9 subdomains (left) and macro mesh used for global calculations (right) . . . . . . . . . . . . . . . . . . . . . . . 3.25 Macro strain field in the xdirection using: (a) approximated integration scheme, (b) consistent integration scheme and (c) difference between these two fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.26 Comparison of global error between approximated and consistent integration scheme in the case of E i /E m = 10, 10 3 and 10 6 . . . . . . . . . . . . . . 3.27 Strain fields ε 11 (x): (a) reference solution, (b) relocalized CMCM solution for β = 0 and (c) relocalized CMCM solution for β = 0.3 and E i /E m = 10 6 3.28 Stress fields σ 11 (x)(MPa): (a) reference solution, (b) relocalized CMCM solution for β = 0 and (c) relocalized CMCM solution for β = 0.3 and E i /E m = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.29 Evolution of error and computational time in function of number of subdomains for E i /E m = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.30 Evolution of global energy error with respect to the number of macro and micro dofs for E i /E m = 10 6 . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.31 Structure decomposed into 49 × 10 × 10 subdomains and involving 1.3 × 10 9 dofs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.32 Coarse mesh: (a) global view and (b) zoom-in at the center of the structure 3.33 Subdomains chosen for relocalization of strain and stress fields . . . . . . 3.34 Relocalized solution: ε 11 (x) (left) and ε 22 (x) (right) . . . . . . . . . . . . . . 3.35 Relocalized solution: ε 33 (x) (left) and σ 11 (x)(MPa) (right) . . . . . . . . . . 3.36 Relocalized solution: σ 22 (x)(MPa) (left) and σ 33 (x)(MPa) (right) . . . . . . 3.37 Structure geometry and a typical cross section . . . . . . . . . . . . . . . . 3.38 Three loading types of the structure and their corresponding deformed shapes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.39 Decomposition of the microstructure (a) and the block used for local field comparison (b) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.40 Relocalized strain solution in a cross section of the relocalization block for (a) β = 0 and (b) β = 0.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.41 Relocalized stress and strain fields for 3 considered cases . . . . . . . . . . 4.1 Decomposed structure and different meshes used in CMCM for the thermal problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Decomposed structure and different meshes used in CMCM for the thermoelastic problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Reference structure and finite element meshes used for calculations . . . . 4.4 Temperature fields of reference solutions and relocalized solution using CMCM for thermal problem with β = 0 and 1 . . . . . . . . . . . . . . . . 4.5 Temperature gradient fields of reference solutions and relocalized solution using CMCM for thermal problem with β = 0 and 1 . . . . . . . . . . . . . 4.6 Global energy and L 2 errors in function of beta for the thermal problem . 4.7 Local energy error of the thermal problem for different values of β . . . . . 4.8 Total 11-strain fields of reference solutions and relocalized solution using CMCM for thermoelastic problem with β = 0 and 1 . . . . . . . . . . . . . 4.9 11-stress fields of reference solutions and relocalized solution using CMCM for thermoelastic problem with β = 0 and 1 . . . . . . . . . . . . . . . . . . 4.10 Global energy and L 2 errors in function of beta for the thermoelastic problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 5.1 (a) Schematic description of FE 2 method Feyel (1999) and of CMCM . . . . 5.2 Viscopplastic 2D heterogeneous beam: geometry and boundary conditions 5.3 Displacement-force curve of the reference and CMCM solutions over 10 loadings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 X-Strain-stress curve of the reference and CMCM solutions over 10 loadings 5.5 Y-Strain-stress curve of the reference and CMCM solutions over 10 loadings 5.6 Evolution of the X-strain along the line y = H/2 at the eighth loading time 5.7 Micro strain fields ε 11 (x) of (a) reference solution; (b) CMCM solution with β = 0 and (c) CMCM solution with β = 1 . . . . . . . . . . . . . . . . . . . 5.8 Cumulative plastic strain fields ε p e of (a) reference solution; (b) CMCM solution with β = 0 and (c) CMCM solution with β = 1 . . . . . . . . . . . 5.9 22-strain fields of a mode-I crack opening structure: (a) reference solution and (b) CMCM relocalized solution . . . . . . . . . . . . . . . . . . . . . . . B.1 Illustration of the collective communications Snir et al. (1998) . . . . . . . B.2 Divide-and-conquer parallel computing flowchart . . . . . . . . . . . . . . B.3 Master/slave parallel computing flowchart . . . . . . . . . . . . . . . . . . B.4 Meshed structure used for the comparison . . . . . . . . . . . . . . . . . . B.5 Wall time of two paradigms as a function of (a) the number of processors with 16 subdomains; (b) the number of subdomains using 12 processors . B.6 Total memory usage of two paradigms in terms of the number of subdomains112 xii List of Tables

4. 1

 1 Mechanical and thermal properties of the components. . . . . . . . . . . . 5.1 Material properties of 2 phases of the structure. . . . . . . . . . . . . . . . . xiii Notations x : Position vector at the fine scale u(x) : Displacement vector at the fine scale T(x) : Temperature at the fine scale ε : Linearized strain tensor σ : Cauchy stress tensor G : Tensor of second gradient of displacements ∇ε : Strain gradient tensor ∇T : Temperature gradient at the fine scale C(x) : Elastic tensor k(x) : Thermal conductivity tensor C 0 (x) : Auxiliary elastic tensor ε : Macroscopic linearized strain tensor G : Macroscopic second gradient of displacements tensor ∇ε : Macroscopic strain gradient tensor g α : Vector of parameters defining Dirichlet boundary conditions on the boundary of subdomain Ω α h α : Vector of parameters defining Dirichlet boundary conditions for the thermal problem on the boundary of Ω α Ω α : Subdomain α A α (x) : Localization tensor in Ω α [ε(x)] : Vector form associated with ε [σ(x)] : Vector form associated with σ C(x) : Matrix form associated with C(x) Ω e : Domain associated with an element e in the fine mesh Ω E : Domain associated with an element E in the coarse mesh u E : Vector of nodal displacements in one element E of the coarse mesh T E : Vector of nodal temperatures in one element E of the coarse mesh H

  FIGURE 2: Images of fan blades made in 3D woven composites of LEAP 1B (version for Boeing 737) (source: Safran)

FIGURE 3 :

 3 FIGURE 3: Illustration of different scales when studying 3D woven composites (source: Safran)

FIGURE 1 . 1 :

 11 FIGURE 1.1: Example of a regular grid used for the FFT method Zeman et al. (2017)

FIGURE 1 . 2 :

 12 FIGURE 1.2: Comparison between Cosserat, Cauchy continuum and the reference solution of the y-displacement along a horizontal line of a bending structure (Forest and Sab (1998))

FIGURE 1

 1 FIGURE 1.3: Second-order computational homogenization scheme (Geers et al. (2010))

FIGURE 1 . 4 :

 14 FIGURE 1.4: Variation of the component C 2222 of the sixth-order elasticity tensor in function of the volume fraction τ (Tran et al. (2012))

FIGURE 1 . 5 :

 15 FIGURE 1.5: Illustration of one-dimensional basis functions[START_REF] Efendiev | Multiscale finite element methods: theory and applications[END_REF]

FIGURE 1 . 6 :

 16 FIGURE 1.6: Schematic description of the HMM method (Abdulle et al. (2012))

  ; Xu et al. (2010); Anciaux et al. (2012); Talebi et al. (2013); Tu et al. (2014). Recently, a non-overalpping coupling technique has been proposed to deal with nonseparated scale heterogeneous structures Wangermez et al. (2020), with incorporating second-order effects in global-local analysis. The main objective of the approach is to construct an approximated solutions only at certain zones of interest (see eg Fig. 1.7).

FIGURE 1 . 7 :

 17 FIGURE 1.7: Illustration of different scales used in the approach Wangermez et al. (2020)

•

  Smoothing: reducing high frequency errors • Residual computation: computing residual after smoothing step • Restriction: downsampling the residual error to a coarse grid • Interpolation: interpolating a correction into a finer grid • Correction: adding prolonged coarse grid solution onto the finer grid The efficiency of the AMG algorithm in solving large systems of equations has been demonstrated in the literature Zecchin et al. (2012); Clees (2005); Geenen et al. (2009);

FIGURE 2

 2 FIGURE 2.2: Illustration of a field in 2D: (a) microscopic; (b) filtered on a coarse mesh

  FIGURE 2.3: Illustration of filtered fields in 1D with different size of the macro mesh

  FIGURE 2.4: Illustration of: (a) a macroscopic structure; (b) a RVE

FIGURE 2

 2 FIGURE 2.5: Positioning of a RVE on the macro structure

FIGURE 2

 2 FIGURE 2.6: (a) Macro mesh of the structure; (b) Macro strain field ε 11 ; (c) macro strain field ε 11 interpolated on a RVE

FIGURE 2

 2 FIGURE 2.7: (a) Structure of the 2D composite beam ; (b) RVE

FIGURE 2

 2 FIGURE 2.8: (a) Macro mesh of the structure ; (b) macro mesh of the RVE

First

FIGURE 2

 2 FIGURE 2.11: Considered woven composite coupon and boundary conditions; the red cylinders are where the pressure fields are prescribed

FIGURE 2 .

 2 FIGURE 2.12: Visualisation of the RVE and its mesh

  FIGURE 2.13: Distribution of pressure applied to the center of the specimen

  FIGURE 3.1: Schematic description of the CMCM method; (a) decomposition of the structure into subdomains; (b) discretization of the structure with a coarse mesh; (c) elementary (off-line) problems to be solved over the subdomains

  FIGURE 3.3: Subdomain Ω α , coarse mesh element E associated with adomain Ω E and fine mesh element e associated with a domain Ω e

  FIGURE 3.4: Subdomain Ω α cut by a coarse mesh element and different integration schemes

  FIGURE 3.5: Remeshing operation in the case where micro elements are cut by a macro element

  FIGURE 3.6: (a) 2D square composite structure: geometry; (b) subdomain used for off-line calculations; (c) subdomain finite element mesh

  FIGURE 3.7: ε 11 (x) of reference solution and obtained with the CMCM method for different coarse meshes

3

  FIGURE 3.9: Structure divided into (a) 1 subdomain, (b) 4 subdomains and (c) 16 subdomains

(

  FIGURE 3.17: Coarse meshes for global calculation of the beam: (a) coarse mesh 1, (b) coarse mesh 2, and (c) coarse mesh 3

  FIGURE 3.22: Minimum displacement (left) and maximum stress (right) using different coarse meshes

  FIGURE 3.23: Random structure geometry (a) and its finite element mesh (b)

FIGURE 3 .

 3 FIGURE 3.27: Strain fields ε 11 (x): (a) reference solution, (b) relocalized CMCM solution for β = 0 and (c) relocalized CMCM solution for β = 0.3 and E i /E m = 10 6

FIGURE 3 .

 3 FIGURE 3.28: Stress fields σ 11 (x)(MPa): (a) reference solution, (b) relocalized CMCM solution for β = 0 and (c) relocalized CMCM solution for β = 0.3 and E i /E m = 10 6

FIGURE 3 .

 3 FIGURE 3.31: Structure decomposed into 49 × 10 × 10 subdomains and involving 1.3 × 10 9 dofs

  FIGURE 3.32: Coarse mesh: (a) global view and (b) zoom-in at the center of the structure

FIGURE 3 .

 3 FIGURE 3.33: Subdomains chosen for relocalization of strain and stress fields

  FIGURE 3.38: Three loading types of the structure and their corresponding deformed shapes

  FIGURE 4.1: Decomposed structure and different meshes used in CMCM for the thermal problem

Chapter 4 .

 4 FIGURE 4.2: Decomposed structure and different meshes used in CMCM for the thermoelastic problem

  FIGURE 4.3: Reference structure and finite element meshes used for calculations

  FIGURE 4.7: Local energy error of the thermal problem for different values of β

  FIGURE 5.2: Viscopplastic 2D heterogeneous beam: geometry and boundary conditions

  FIGURE 5.3: Displacement-force curve of the reference and CMCM solutions over 10 loadings

  FIGURE 5.7: Micro strain fields ε 11 (x) of (a) reference solution; (b) CMCM solution with β = 0 and (c) CMCM solution with β = 1

  FIGURE 5.9: 22-strain fields of a mode-I crack opening structure: (a) reference solution and (b) CMCM relocalized solution

  Algorithm for the Divide-and-conquer paradigm

  FIGURE B.3: Master/slave parallel computing flowchart

  .recv(source=MPI.ANY_SOURCE, tag=tag, status=status) (receive result of any slave) 11: comm.send(new_data, dest=status.Get_source(), tag=tag) Send new job for the 12: non-working processor 13: END WHILE 14: LOOP over all N slaves Tell all the slaves to exit by sending an empty message with a 15: specific tag "0" 16: comm.send(0, dest=rank, tag=0) 17: END LOOP Algorithm 5: Algorithm for the master in Master/slave paradigm 1: status = MPI.Status() 2: WHILE true DO 3: data = comm.recv(source=0, tag=MPI.ANY_TAG, status=status) (result, dest=0, tag=0) Send result back to the master 10: END WHILE Algorithm 6: Algorithm for the slave in Master/slave paradigmB.3 Performance comparison of two parallel computing paradigmsIn this section, to compare the performance of two parallel computing paradigms presented in the previous section, I consider a meshed structure depicted in Fig.B.4. The mesh contains 367,496 linear tetrahedral elements, corresponding to 195,057 degrees of freedom. For the sake of simplicity, I only compare the performance of two paradigms on solving local problems in parallel. First of all, to investigate the performance of the two paradigms as a function of the number of processors, I divide the mesh into 16 subdomains with approximately the same size. Next, in order to investigate the performance as a function of the number of subdomains, I keep the number of processors at 12 and vary the number of subdomains. The curves describing the computational time in function of the number of processors and subdomains are plotted in Fig. B.5. The memory usage spent using the two parallel computing paradigm are illustrated in Fig. B.6.

2

  Images of fan blades made in 3D woven composites of LEAP 1B (version for Boeing 737) (source: Safran) . . . . . . . . . . . . . . . . . . . . . . . . . 3 Illustration of different scales when studying 3D woven composites (source: Safran) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1.1 Example of a regular grid used for the FFT method

1

  Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Mechanical properties of the matrix. . . . . . . . . . . . . . . . . . . . . . . 2.2 Mechanical properties of the woven threads. . . . . . . . . . . . . . . . . . 2.3 Global error on selected quantities of interest of 3 homogenization schemes. 2.4 Parameters used in TexGen for the generation of the texture form. . . . . .

2.5 Mechanical properties of the matrix. . . . . . . . . . . . . . . . . . . . . . . 2.6 Mechanical properties of the woven threads . . . . . . . . . . . . . . . . . . 3.1 Number of elements and dofs for each coarse mesh. . . . . . . . . . . . . . 3.2 Number of elements and of dofs for each mesh used in the calculation (T3: linear triangular element; Q4: bilinear quadrangular element) . . . . . . . 3.3 Material properties of concrete . . . . . . . . . . . . . . . . . . . . . . . . .

TABLE 1 :

 1 Notations

  Table 2.1).

	Young's modulus E (MPa) Poisson ratio ν
	4000	0.3

TABLE 2 .

 2 1: Mechanical properties of the matrix.-The woven threads are orthotropic elastic (see Table2.2).

	E 11 (MPa) E 22 (MPa) E 33 (MPa) ν 12 ν 23	ν 31	G 12 (MPa) G 23 (MPa) G 31 (MPa)
	194400	8200	8200	0.3 0.3 0.0126	7000	3100	7000

TABLE 2.2: Mechanical properties of the woven threads.

TABLE 2

 2

.5: Mechanical properties of the matrix.

  ε 33 , ε 12 , ε 13 , ε 23 , ∇ε 112 , ∇ε 113 , ∇ε 221 , ∇ε 223 , ∇ε 331 , ∇ε 332 , G 123 , G 213 , G 312

	T

TABLE 3 .

 3

1: Number of elements and dofs for each coarse mesh.

  .2).

		subdomain Coarse mesh 1 Coarse mesh 2 Coarse mesh 3
	Element type	T3	Q4	Q4	Q4
	Number of elements	5000	21	84	336
	Number of dofs	5202	88	258	850

TABLE 3 .

 3

2: Number of elements and of dofs for each mesh used in the calculation (T3: linear triangular element; Q4: bilinear quadrangular element)

TABLE 4 .

 4 1: Mechanical and thermal properties of the components.

Table 5

 5

	.1.

Meso is the scale of the woven threads in woven composites, whereas the micro scale is a smaller scale where carbon fibers are considered.

(a) Global energy error (b) L 2 error

(a) Reference solution (b) β = 0.16 (c) β = 0.32

(a) Reference solution (b) 1 subdomain (c) 4 subdomains (d) 16 subdomains

(a) Reference solution (b) CMCM solution with β = 0 (c) CMCM solution with β = 0.3

(a) Evolution of global error as a function of the number of subdomains (b) Computational time during online calculations as function of number of subdomains

(a) Structure decomposed into 16x16x16 subdomains. (b) A block used for relocalization and post-treatment.

Acknowledgements

Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures The difference in terms of macroscopic strain solution between consistent and approximated integration schemes is shown in Fig. 3.25 (see section 3.5). Fig. 3.26 shows the evolution of the global error in the energy norm as a function of β for three cases of property contrast, using both consistent and approximated integration schemes. When the mechanical contrast ratio E i /E m is high, a high value of the error is obtained. However, increasing the value of β decreases the global error. It is worth noting that in the case of 9 subdomains, the global error curves converge when β = 0.3 for E i /E m = 10 6 , and β = 0.2 for E i /E m = 10 and 10 3 . In addition, it can be seen that using a consistent integration scheme does reduce the global error in all three cases. However the gain is not significant in view of the added complexity brought by the consistent integration. For this reason, I apply the approximated integration schemes in the next examples, and more specifically in the 3D examples. The total simulation time was 20 minutes for each subdomain off-line calculation (performed in parallel) and the coarse mesh calculation took roughly 5 days, on 32 cores.

Chapter 3. Coarse Mesh Condensation Multiscale method for linear heterogeneous structures

Then, the present method has a very high potential to study fully detailed composite structures without any assumptions on scale separation and periodicity.

3D tomography of concrete

In this example, the CMCM is applied to a structure generated from micro tomography, including a large number of dofs to show the potential of the method. The structure is a sample made of concrete (see Fig. 3.37). The sample contains 3 phases, sand grains, plaster matrix and voids. The properties of each phase are provided in Table 3.3. To apply the proposed method, the structure is decomposed into 4096 subdomains as depicted in Fig. 3.39. The corresponding mesh for the whole structure contains 148,604,196 elements, corresponding to 448,351,500 dofs. The coarse mesh contains 390,462 degrees of freedom, which is approximately 1000 times less than the reference microstructure. To limit computational costs, I have only investigated β = 0 and β = 0.2. In this example where all subdomains are different, the computational times related to off-line calculations are quite expensive. However, once off-line calculations have been conducted, the coarse mesh problem can be solved for several loading cases at low computational costs.

In this example, I have tested 3 types of loading: tension, shear and concentrated load, as depicted in Fig. 3.38. Chapter 5. Extension of the CMCM method to nonlinear problems least-square sense. Such distance is expressed as:

where ∆u e is the macro strain field of macro element e. Minimizing J gives us a relation between the parameter vector g α and the macro displacement of the coarse mesh:

where H α and L αE can be computed from (3.32) and (3.33), respectively.

Then the relation between the micro strain field and the macro displacement field can be deduced as:

where T α (x) is computed from (3.35). Substituing (5.19) in the linearized problem (5.10) I obtain a linear system of equations to be solved on the coarse mesh:

where

We can note that in order to solve the problem on the coarse mesh at each iteration, it requires the knowledge of C tan , which means that the reconstruction of the micro displacement field is nescessary. The construction of such field can be found in Appendix A.

Once the global problem is solved, the stress in the whole domain is known and the internal variables of the nonlinear problem can be updated.

Algorithm of the CMCM method in a nonlinear context

The overall algorithm for the CMCM solving procedure in the nonlinear context is summarized in Algorithm 3.

Appendix A

Reconstruction of the local displacement field

The CMCM method allows all local fields to be reconstructed, including displacement field in all subdomains, as described in the following. I denote by x a local coordinate system related to one subdomain. For example, in the case where each subdomain is associated with a periodic unit cell, the local coordinate system might include an origin at the center of the unit cell. The displacement field is the global coordinate system x is expressed as:

where ũ(x) denotes the remaining fluctuations of the local displacement, which is non zero in subdomains containing heterogeneities; û(x) is the displacement related to the rigid body motions (translation t and rotation ω × r P with r P is the vector of any point on the axis of rotation to P). This term is determined following the fact that the micro and macro domains have the same rigid body motions, such that:

and

However in the presented numerical examples, û(x) is very small. Therefore, it is neglected to avoid numerical complications.

The localization problem being solved in the local coordinate system attached to the subdomain, we have:

where U α (x ) = u 1 (x ), u 2 (x ), ..., u N g (x ) . From (A.1) and (A.4) we have:

Substituting (A.5) in (A.1) we obtain:

Appendix A. Reconstruction of the local displacement field

And so:

Implementation of parallel computing in Python

In this appendix, I summarize the algorithm to implement the CMCM method in Python using parallel computing with MPI. As introduced in the previous chapters, the CMCM method is well suited for parallel computing in at least two of its main steps. In the next sections, a brief introduction to MPI and its application to our problem will be presented.

B.1 Introduction to parallel computing using MPI in Python

Message Passing Interface (i.e. MPI) is a standardized and portable message-passing tool used to connect between processors in a machine or between different machines in a cluster. MPI can be applied to many programming languages for parallel computing such as C, C++, Fortran or Python. In this Ph.D., I have chosen Python for programming and the "mpi4py" open-source package for implementing MPI in Python.

B.1.1 Collective communications

Collective communications transmit data among all processors. The collective communications of MPI can be listed as follows:

• Broadcast which broadcast data from one processor to all the processors

• Gather which collects all data from all processors to one processor

• Scatter which scatter data from one processor to all processors

• Scatter/Gather data from all processors to all processors (or alltoall).

These collective communications are illustrated in

B.1.2 Inter-communication

Before communicating between processors, a group of communicators, usually called comm needs to be created. In Python, after importing the MPI module from the mpi4py package using the command:

the group of communicators "comm" can be defined using:

Inter-communication, as its name represents, is the way to communicate between processors. In this code, there are two main inter-communications that are used to connect between processors which are send and receive. These inter-communications are often used in parallel computing paradigms whose processors have different roles. One simple example is the Master/slave paradigm, where the master's role is to pre-process and distribute data while the slave's role is to solve the problems. The send communication is used when the master wants to send a problem setup to a slave, or a slave wants to send a problem result back to the master: comm . send ( data , dest = processor_rank , tag = data_tag ) where data is the data needed to be sent from one to another processor, dest stands for destination, which is the rank of the receiving processor and tag is the tag of the data in case where there are several data packages being sent to one processor. The rank of each processor can be obtained locally using the command:

The receive communications is used when a slave wants to receive a problem setup from the master or the master wants to receive a problem result from a slave: result = comm . recv ( source = processor_source , tag = tag_result , status = processor_status )

where source is the processor where the result comes from and status is the current status of the receiving processor.

B.2 Parallel computing paradigms

There are many parallel computing paradigms that have been proposed in the literature. For example, from [START_REF] Pritchard | Mathematical models of distributed computation[END_REF]; [START_REF] Hansen | Model programs for computational science: A programming methodology for multicomputers[END_REF]; [START_REF] Wilson | Parallel programming for scientists and engineers[END_REF], some parallel computing paradigms can be filed in the following categories:

• Processor farms based on replication of independent jobs

• Geometric decomposition based on the structure of the input data

• Algorithmic parallelism based on the use of data flow

• Divide-and-conquer based on multi-branch recursion

• Master/slave

In this Ph.D., I have chosen two paradigms which are "Divide-and-conquer" and "Master/slave". The differences in definition and performance between these two algorithms will be investigated in the next subsections.

B.2.1 Divide-and-conquer paradigm

The Divide-and-conquer paradigm is a parallel computing algorithm based on multibranch recursion. The paradigm works by recursively breaking down a problem into two or more similar sub-problems until they are able to be solved directly (see e.g. The distribution of the local problems onto processors can be done randomly. One of the main difficulties of this paradigm is that, when the sub-problems do not have the same size, there might be extra waiting time between resolutions of sub-problems. More specifically, if processor 1 receives two sub-problems with small size, while processor 2 receive two problems with bigger size, then evidently processor 1 will need to wait for processor 2 to finish its job before sending back the results. Another difficulty of this paradigm is the distribution of data onto processors. Since there is no master processor to operate the decomposition of data, the initial data needs to be distributed to all processors. Therefore, when dealing with problems containing large input data, this paradigm often exhibits poor performance. On the other hand, the advantage of this parallel computing algorithm is that it is very easy to implement, and the communications between processors are simple. A sample implementation of the Divide-and-conquer paradigm is presented in Algorithm 4. It can be seen that in terms of wall time, the two paradigms exhibit approximately the same performance when I vary the number of processors or the number of subdomains. On the other hand, the memory usage of two paradigms is very different (see Fig. B.6). The Master/slave paradigm uses a lot less memory than Divide-and-conquer. Indeed, since the Master/slave paradigm uses a dedicated processor for pre-processing data, said data is partially distributed to all processors; unlike the Divide-and-conquer paradigm, where the pre-processing of data needs to be done on each processor.

Appendix B. Implementation of parallel computing in Python

B.4 Conclusion

In this appendix, I have introduced the implementation of the parallel computing of the CMCM method in Python, using two parallel computing paradigm namely Divide-andconquer and Master/slave. From a comparison using a 3D example, it can be concluded that in terms of computational time, the two chosen paradigms exhibited approximately the same performance. However, in terms of memory usage, the Divide-and-conquer paradigm used a lot more memory than the Master/slave paradigm, because the initial data needs to be broadcasted to all processors for pre-processing. For example, if we have a large input data, typically in industrial applications with hundreds of gigabytes, the computer memory might not be sufficient. In conclusion, in this Python code, the Master/slave has been chosen to be the default parallel computing paradigm. The Divide-and-conquer paradigm is still an option for users, and can be used for example for small-size problems, such as 2D academic problems.