Renée El

Melhem Bât

Pascal Blaise

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Philippe Delachartre
email: philippe.delachartre@insa-lyon.fr

M Philippe
email: philippe.normand@univ-lyon1.fr

Mme Sylvie

M Hamamache Kheddouci

M Stéphane Benayoun
email: stephane.benayoun@ec-lyon.fr

Stéphanie Cauvin

M Jocelyn Bonjour

M Christian Montes

Keywords: parametrized architectured material, homogenization, uncertainty quantification, multiscale topology optimization, robust opimitization

Topology optimization is a systemic design that requires simulation and optimization of a system for a single or multiple physics coupling processes. However, it is limited in the engineering sense regarding the absence of uncertainties and limitations on applied monophase material.

The foundation of this dissertation is to combine homogenization and stochastic processing into topology optimization to formulate a robust multiscale topology optimization approach. Accordingly, this Ph.D. dissertation concerns (1) the multiscale and multiphysics performance of heterogeneous materials/structures embedded with microstructured material, taking into account the uncertainties, (2) for further optimizing the heterogeneous structure at different scales to satisfy target performance.

These microstructures may arise from the processing of biological materials, or from dedicated engineered materials, e.g., aerogels, foams, composites, acoustics metamaterials, etc. We parametrize architectured material; study the performances of the microstructure at the macroscopic scale by homogenization method. Then, the homogenization model can be considered as a stochastic model with uncertainties exhibited in the unit cell. It can be built from a polynomial chaos development. In addition, these parametrized micro geometry features can be mapped into homogenized properties space, which can be utilized as design variables to control the macrostructure performance.

Afterward, we combined the topology optimization, homogenization, and uncertainty quantification to (1) design macro topology and micro material distribution to maximize structure stiffness (2) reduce the structure sensitivity to presented uncertainties (e.g., loading and material properties). This proposed general framework has the advantage and compatibility ability in solving optimization problems considering the (1) multiple parametrized architectured cells, (2) complex loading problem, (3) hybrid uncertainty, etc., with an affordable computation manner.

To my family, my son Yuanlang.

L'optimisation de la topologie est une conception systémique qui nécessite la simulation et l'optimisation d'un système pour un ou plusieurs processus de couplage physique. Cependant, il manque de sens technique en ce qui concerne l'absence d'incertitudes et de limitations sur le matériau monophasé appliqué. Le fondement de cette thèse est de combiner l'homogénéisation et le traitement stochastique dans l'optimisation de la topologie pour formuler une approche robuste d'optimisation de la topologie à plusieurs échelles. En conséquence, ce doctorat. La thèse porte sur (1) les performances multi-échelles et multiphysiques de matériaux/structures hétérogènes noyés dans des matériaux à microstructures, en tenant compte des incertitudes, (2) pour optimiser davantage la structure hétérogène à différentes échelles pour satisfaire les performances cibles. Ces microstructures peuvent provenir du traitement de matériaux biologiques ou de matériaux d'ingénierie dédiés, par exemple des aérogels, des mousses, des composites, des métamatériaux acoustiques, etc. Nous paramétrons le matériau d'architecture ; étudier les performances de la microstructure à l'échelle macroscopique par la méthode d'homogénéisation. Ensuite, le modèle d'homogénéisation peut être considéré comme un modèle stochastique avec les incertitudes présentées dans la cellule unitaire. Il peut être construit à partir d'un développement de chaos polynomial. De plus, ces caractéristiques de microgéométrie paramétrées peuvent être mappées dans un espace de propriétés homogénéisé, qui peut être utilisé comme variables de conception pour contrôler les performances de la macrostructure. Par la suite, nous avons combiné l'optimisation de la topologie, l'homogénéisation et la qualification des incertitudes pour (1) concevoir la macrotopologie et la distribution des micromatériaux pour une rigidité maximale de la structure (2) réduire la sensibilité de la structure aux incertitudes présentées (par exemple, le chargement et les propriétés des matériaux). Ce cadre général proposé a la capacité d'avance et de compatibilité pour résoudre les problèmes d'optimisation en considérant [START_REF] Abueidda | Mechanical properties of 3D printed polymeric cellular materials with triply periodic minimal surface architectures[END_REF] les cellules d'architectures paramétrées multiples, (2) le problème de chargement complexe, (3) les hybrides non certifiés, etc., avec un mode de calcul abordable.

Mots-clés : matériau d'architecture paramétré, homogénéisation, qualification d'incertitude, optimisation de topologie multi-échelle, optimisation robuste Chapter 1

Motivation and research objective

Many natural materials such as wood, bones, shells, etc., synthesize spatially varying material properties and structural hierarchy to achieve appealing functions. Inspired by this, numerical hierarchical topology optimization provides a tool for optimizing the macro topology and the microarchitecture of the material (shown in figure 1.1), which has more advanced performance than a mono scale topology optimization. It is due to increasing the design freedom in microarchitectured material. In addition, additive manufacturing (AM) provides a possible tool to fabricate such a delicate multiscale structure.

While the power of the multiscale topology optimization is demonstrated, deterministic conditions are assumed, meaning that the various uncertainties in real engineering problems are disregarded. This will lead to that designs being usually less valuable in real applications.

To this end, we take into account the uncertainties into multiscale topology optimization to formulate the robust optimized macro topology, and complex microarchitecture, leading to more robust design with the presence of uncertainties.

Topology optimization

Topology optimization (TOP), first developed by Bendsose and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF], is a systemic design that requires simulation and optimization of a system for single or multiple physics coupling processes. Within active engineering domains, the most common cases are: mechanical, thermal, fluid flow [START_REF] Dilgen | Topology optimization of turbulent flows[END_REF], thermo-mechanical coupling [START_REF] Desmorat | Structural rigidity optimization with an initial design dependent stress field. Application to thermo-elastic stress loads[END_REF], electrical-mechanical coupling [START_REF] Canfield | Topology optimization of compliant mechanical amplifiers for piezoelectric actuators[END_REF] (piezoelectric materials) heat-flow coupling [START_REF] Dede | Multiphysics topology optimization of heat transfer and fluid flow systems[END_REF], etc. The popularity of the TOP is due to its ability to provide a novel, 'good' optimum design without any geometry constraints, compared to shape optimization or parametric optimization. It provides the primary design that Additive Manufacturing can produce. The different approaches to TOP relied on the description of the structure topology. They are the Density Description TOP and the Boundary Description TOP (Level Set [117], Method Morphable Component/Void [START_REF] Guo | Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[END_REF]). Herein, we mainly focused on Density-based topology optimization; the TOP community complained that this method is a nightmare for large-scale topology optimization since the number of design parameters is equivalent to the number of elements in the spatial discretization. However, it is has been dramatically improved by, on the one hand, parallel computing. On the other hand, the development of multiscale topology optimization providing plenty of microarchitecture details will enhance the performance of macro topology. Multiscale topology optimization is an alternative way to solve large and full-scale topology optimization, which is one of the topics we are interested in.

Density based topology optimization

The basic idea of TOP is to seek the best material layout by deciding the density ρ of each material point or the integration points on the design domain Ω to be either 0 (void) or 1 (solid) while satisfying the constraints. In this case, this problem can be written as

min ρ∈{0,1} c(ρ) = U T F s.t. KU = F V (ρ) < V * (1.1)
where the compliance c is classically used as a cost function to test TOP algorithms, and will be used herein. F and U are the loading condition and structural displacement response. K is the structural stiffness matrix. V * is the volume fraction of solid material ratio to the global design space size Ω, it is classically used as constraints in TOP algorithms. Directly solving such a problem is rather challenging, because firstly, the binary representation of macro topology is not well-posed and lacks optimum. Secondly, the objective function is non-differential with respect to the design variable (the density ρ). A possible solution is the gradient-free approaches (e.g., genetic algorithm, Particle Swarms [START_REF] Rao | Engineering optimization: theory and practice[END_REF]) that only require the evaluation of the objective function and the constraints to search towards a possible global optimum heuristically. The drawbacks of such optimization obviously are the computation cost. Often, it requires high order of magnitude objective function realization. It seems not a brilliant choice for large sale topology optimization, with millions of design variables.

Material interpolation method

A gradient-based method requires the gradient descent or the sensitivity information of objective function to the design variables, enabling it to fast search towards a 'good' but not guaranteed global optimum. First, for obtaining a differentiable objective function, the ρ is relaxed between 0 and 1. However, the existence of intermediate pseudo-density ρ lacks physical interpretation; a critical aspect of the density-based topology optimization is to steer the design to converge towards a binary well-connected macro topology representation. To this end, a most popular and well-developed method called the Simple Interpolation Material Penalized (SIMP) method [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF], where the intermediate density can be mathematically interpreted as a material, which the corresponding Young's module E ρ can be expressed as:

E ρ = ρ p E 0 (1.2)
where, E 0 is the solid material Young's modulus, p is a penalty factor. Literature [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF] shows that, when p = 1, the classical minimization of the compliance problem is convex, it lose the convexity when p > 1. It indicates that the global optimum is not easy to search for the appearance of many local optimum. However, p = 3 is a common choice that can guarantee the design converge binary (an 0/1 density represented structure) [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF]. A common practice to enhance the possibility to find the global optimum is the multiple starts from different initial guesses. Another way is to continuously increase the p from 1 to a target value. Other similar material interpolation methods are also proposed:

• Bendsøe and Kikuchi [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] proposed a homogenization method to model the intermediate density elements using periodically distributed microstructures, but the optimum design is far from easy manufacturing.

• Another power law-like model is nicknamed as Rational Approximation of Material Properties (RAMP) [102]. The idea behind RAMP is to eliminate the zero gradients in the low-density elements to speed up the optimization process. Comparison of the performance of different RAMP, SIMP, and Homogenization interpolation methods in topology optimization have been investigated in [START_REF] Dzierżanowski | On the comparison of material interpolation schemes and optimal composite properties in plane shape optimization[END_REF].

• With the development of peak function material interpolation scheme [128], multiple materials can easily be incorporated into the topology optimization without increasing the number of design variables.

• Another multiphase material (assuming m candidate materials) interpolation method is called the Ordered multi-material SIMP interpolation [135]. However, it will introduce m(m -1)/2 design variables in a single element, leading the computational cost to a heavy burden.

• For extension of SIMP to nonlinear material interpolation, reference [START_REF] Maute | Adaptive topology optimization of elastoplastic structures[END_REF] incorporated von Mises plasticity model into SIMP, where three plasticity model parameters were penalized by different factors.

Optimization algorithms

The Relaxation and Penalization density approach (R/P) approach will attempt to drive the design variable (i.e. the density ρ) almost completely binary. The main advantage of this approach is that, since it obtained a differentiable objective function and sensitivity information, some gradient-based non-linear programming (NLP) [START_REF] Boggs | Sequential quadratic programming[END_REF][START_REF] Forsgren | Primal-dual interior methods for nonconvex nonlinear programming[END_REF][START_REF] Rojas-Labanda | Benchmarking optimization solvers for structural topology optimization[END_REF] can be used to interactively update the design variables until converge criteria are satisfied. (e.g. threshold the change of two successive variables or objective function values, maximum iteration number, and KKT condition). In TOP community, the most commonly used optimization schemes are the Optimality Criteria (OC) method [112], the Method of Moving Asymptotes (MMA) and the Globally converged MMA (GCMMA) [107]. However, the primary drawback of OC method is that it cannot handle multiple constraints [START_REF] Haftka | Elements of structural optimization[END_REF]; nevertheless, MMA can tackle this problem. However, the gradient based optimization algorithm has a common issue with the local optimum. This may be overcome by using the gradient-free global search algorithms, i.e. genetic algorithm (GA), simulated annealing (SA) ant colony, etc.

Filter method

Gradient-based density topology optimization with SIMP material interpolation often suffers from two common numerical instability problems: checkerboard pattern (figure 1.2(b)) and mesh dependency (figure 1.2(c)). The checkerboard pattern corresponds to a scattered structure with too many small holes. Several techniques have been proposed to overcome the checkerboard problem by high-order elements in FEM simulation and filtering restriction [START_REF] Jog | Stability of finite element models for distributed-parameter optimization and topology design[END_REF].

Filtering has excellent advantages in restricting the R/P convergence towards a continuously connected design structure and overcoming the mesh-dependency issue. In standard density based topology optimization, the ρ j is the design variable associated with each element j. When using filter method, ρj is introduced as a filtered density, ρ = F (ρ). Then the filtering ρj , is then used to interpolate material properties in equation 1.2. The chain rule, therefore, solves all derivatives to the density, reads as following:

∂ * ∂ ρ = ∂ * ∂ρ ∂ρ ∂ ρ (1.3)
A linear density filter [START_REF] Bourdin | Filters in topology optimization[END_REF][START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF] is common used:

ρ i = j w ij ρ j (1.4)
where ρ is the design variable sets, subscript i denote the central element, subscript j denotes the j element within a certain distance to element i, w ij are the weighting factors. Conic weights w ij is defined as:

w ij = R-d(i,j) k∈N i (R-d(i,k)) j ∈ N i 0 j / ∈ N i (1.5)
where R is the prescribed filter radius. Some alternative filters are existing, for example

• Sensitivity filter is the initially proposed method [START_REF] Sigmund | On the design of compliant mechanisms using topology optimization[END_REF], widely used in early literature. It is not a density filter in a sense indicated above. However, the objective funtion c is filtered using: ∂ c ∂ρ i = j ∂c ∂ρ j w ij ρ j ρ i (1.6)

• Heaviside filter [START_REF] Guest | Achieving minimum length scale in topology optimization using nodal design variables and projection functions[END_REF] is proposed for obtaining 0/1 binary solutions, based on continuous approximation of the Heaviside function.

1 -ρi = 1 -e -β(1-j w ij ρj ) + e -β 1 -

j w ij ρj (1.7)
where the β > 0 is the control parameter to approximate the continuous Heaviside step, ρ is the filtered density using equation 1. [START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF]. Then results should be mapped to 0/1 physical density. Whereas β = 0, it is the exact linear density filter.

• Helmholtz-type PDE [START_REF] Lazarov | Filters in topology optimization based on Helmholtztype differential equations[END_REF] achieves a similar effect of mesh independence to sensitivity and density filtering. The PDE filter method saves the need of storing the neighboring information, therefore, has a significant improvement for parallel computation.

Extensions of the Density TOP in Mechanics application

To summarize, the density-based TOP solved by gradient-based optimization algorithm consists of three main parts, i.e. structural analysis, sensitivity analysis, design variables updated, as shown in figure 1. [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF]. Each primary step has several related sub-steps to improve TOP performance and overcome the so-called checkerboard pattern and mesh-dependency issues. This framework will be extended in the following section to solve a multiscale multiphase multiphysics robust topology optimization. Density-based topology optimization is easy to implement and compatible with different FEM softwares. Moreover, minimizing and/or constraint compliance problems have been well investigated. The figure 1. [START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF] shows the flowchart for development and application in mechanical design. In minimizing compliance problems, a typical issue corresponding to stress concentration may appear. However, minimizing the globally maximum quantity is a challenge, relying on obtaining a differentiable objective function. Refs [START_REF] De Leon | Stress-constrained topology optimization for compliant mechanism design[END_REF] adopted the pnorm, p-mean, and the Kreisselmeier-Steinhauser (K-S) functions as aggregated measures of the maximum stress. Consequently, the differential mini-max function (minimize the global maximum quantity) can be mathematically formulated. Different quasi-static failure, fatigue criteria [START_REF] Huang | Evolutionary topology optimization of continuum structures with an additional displacement constraint[END_REF] are embedded into the sensitivity analysis, preventing damage in structure within loading conditions. References [START_REF] Fedulov | Residual strength estimation of a laminated composite with barely visible impact damage based on topology optimization[END_REF]113] consider the material degradation, where stress exceeds a given stress, in the density-based topology optimization. Consequently, fail-safe design [START_REF] Huang | Evolutionary topology optimization of continuum structures with an additional displacement constraint[END_REF][START_REF] Jansen | Topology optimization of fail-safe structures using a simplified local damage model[END_REF][START_REF] Kranz | An empirical study on stress-based fail-safe topology optimization and multiple load path design[END_REF]115] which incorporates the failure uncertainties e.g. unknown various failure shapes, sizes, and locations into the topology optimization has been developed to obtain still a good design while damage happens.

Multiscale topology optimization

Most of the TOP is limited to designing single-phase material or multi-phase structure. However, the macrostructure embedded with complex architectured material has been attractive recently because of a significant advance in engineering application. Such structural materials are commonly exhibited in nature. Such as wood, bones shells, which exhibit spatially varying material properties on macro-scale due to the graded architecture in micro-scale. Typically, the multiscale design structure has the desired performance compared to a single material structure. For example, the reinforced composite materials, sandwich panels, and foam structure are the engineered multi-scale systems used to enhance structural, mechanical performance. Therefore, multi-scale optimization is getting attractive recently and significantly contributes to real engineering applications.

Many works have made outstanding contributions to this topic. We categorize the techniques into three different categories: full-scale topology optimization, sequential multiscale topology optimization, and concurrent multiscale topology optimization. They may exist different sub-techniques for different methods.

Checkerboard indeed occurs when optimization is performed on an intermediate scale design mesh problem (in figure 1.2 (b)), however, for full-scale approach optimization on a fine fixed mesh and in the absence of regularization (filtering density), the design with local architecture details appears naturally. In addition, the appearance of fine microstructure can be stimulated by local constraints. It is used to control the maximum length of material, avoiding one phase being accumulated to a large region. The local volume constraints have been adopted in full-scale topology optimization to design wherein micro-scale, there are continuous fiber structure [START_REF] Li | Full-scale topology optimization for fiberreinforced structures with continuous fiber paths[END_REF], bone-like porous structure [120], and shell infilled uniform, and non-uniform porous structure [121]. Note that it creates a large number of constraints (number of elements) to the optimization problem.

To avoid intensive computation on a full-scale simulation, scale separation is adopted herein; therefore, a multiscale topology optimization approach is proposed. An approach to implement the multiscale topology is to consider the macro topology and micro topology optimization sequentially [START_REF] Imediegwu | Multiscale structural, thermal and thermo-structural optimization towards three-dimensional printable structures[END_REF][START_REF] Schury | Efficient two-scale optimization of manufacturable graded structures[END_REF]. In that way, the macrostructure is designed using Free Material Optimization. Once the macro problem is solved, then the inverse homogenization [START_REF] Sigmund | Materials with prescribed constitutive parameters: an inverse homogenization problem[END_REF] is used to decide the macro architecture.

An alternative approach is concurrent multi-scale topology optimization. Since many works have different ways on optimizing microarchitectured materials, we categorize them into the online optimized microstructure, the offline database, and the parametric cell.

• Online optimized microstructure.

Since the two-scale separation is adopted, the optimization formulation involves one global macro problem at structure scale and many local problems at micro-scale. In each optimization iteration, the global problem is used to optimize the macro topology of density or the level set representations, the local problem designs the architecture to meet the desired homogenized property [START_REF] Sivapuram | Simultaneous material and structural optimization by multiscale topology optimization[END_REF]123,124]. The online optimized microstructure material has the most design freedom because it allows any number of unique cell structures applied to the design problem.

• Parametric cell.

Predefined cell structures are controlled by one or several geometric parameters, which have been applied to design multi-scale structures efficiently [START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF][START_REF] Maconachie | SLM lattice structures: Properties, performance, applications and challenges[END_REF]. Such parametric cells, for example, as shown in figure 1.4, are the Triple Periodic Minimum Surface(TPMS), the different types of the strut-based lattice structure and the optimized architecture to achieve the maximum bulk module with respect to different volume fractions. Moreover, one can refer to literature [START_REF] Metz | Towards multi-dynamic mechano-biological optimization of 3D-printed scaffolds to foster bone regeneration[END_REF] for different 3D-printed scaffold materials used in mechanobiological.

Due to the complexity of microarchitectures, FEM based homogenization is used to evaluate offline effective material properties of the parametric cell; afterward, one can construct semi-analytically surrogate models to FEM-based homogenization model by Response Surface Methodologies [START_REF] Imediegwu | Multiscale structural, thermal and thermo-structural optimization towards three-dimensional printable structures[END_REF]. It is convenient to construct a differential function that maps the geometry parameters to homogenized properties, leading the optimization efficient, but with significant optimized space limitation to specific parametric cells. In that way, the micro-scale optimization is no more extended topology optimization but a kind of parametric optimization.

Hybrid types of parametric cells can be embedded in the same design system; however, an important material transition approach [START_REF] Sanders | Optimal and continuous multilattice embedding[END_REF] is used to smooth the adjacent but different types of microarchitectures.

• Offline material database: a data-driven approach. We think the offline material database is somehow an extension technology development to the online optimizing microstructure approach and parametric cell approach mentioned above, rather than a unique novel method for solving the multiscale topology optimization problem. For example, [START_REF] Djourachkovitch | Multiscale topology optimization of 3D structures: A micro-architectured materials database assisted strategy[END_REF] indeed apply the online optimizing approach, however, in which the online optimized elastic tensor is mapped to a corresponding cellular material in the database, which is preliminarily obtained by a sufficient number of microscopic topology optimization.

[116, 133] are the parametric cell approach with multiple microscopic design parameters (5 and 3 parameters, respectively), to control either isotropic or anisotropic material properties by a so-called data-driven surrogate homogenization model training by machine learning method. In the context of data-driven multiscale topology optimization, the parametric microcell is not explicitly described, but the mixing of both quantitative (geometry variables of a specific unit cell) and qualitative variables (class of the unit cell). In addition, we think the data-driven approach is a new tool to design architectured material by mixing; finally, we can achieve an extended continuous material space and well-connected geometry, function graded micro architectured material.

We summarize the multi-scale topology optimization methods in terms of scale separation and different kinds of microarchitecture optimization. Full-scale optimization is the most straightforward way, and the local constraints consider the micro features (fiber or porosity, etc.) and the micro characteristic length leading the number of constraints to the problem equivalent to the number of elements. Multi-scale optimization is more efficient. One optimizes the microarchitecture which provides more freedom to the multiscale scale optimization, for one can design unique cell for every macro material point; however, it is less efficient, and poses more restriction issues. Using the parametric cell will reduce the design space, but it is the most efficient way, and this has been widely applied for bio-mechanic design, such as human-bones [START_REF] Metz | Towards multi-dynamic mechano-biological optimization of 3D-printed scaffolds to foster bone regeneration[END_REF].

However, the topic is not limited to this. Several critical review papers have categorized the multiscale topology optimization into hybrid or non-hybrid micro architectures [START_REF] Gao | Concurrent topology optimization of multiscale composite structures in Matlab[END_REF], and the different restrictions (material distribution, shape, connectivity, or orientation of the unit cell etc.) to microarchitecture [122].

Topology optimization considering uncertainties

In the deterministic topology optimization, one seeks the design variables X to minimize the structural performance f (X) (now denoted as f for a general purpose, not restricted to compliance) subject to the environment loading and the constraints. However, the uncertainties source ξ in material, geometry or excitation can be accounted for in topology optimization; the design objective function f (X, ξ) has a probability description function. Therefore, non-deterministic topology optimization is proposed to reduce the system performance's sensitivities to the uncertainties, while not scarifying a lot to the deterministic TOP (DTO) objective function. The approach for taking into account the uncertainties is categorized as either (1) worst-case orientated or (2) stochastic, [START_REF] Ben-Tal | Robust optimization-methodology and applications[END_REF]109].

The worst case is the common practice technology to minimize the objective function subjected to the worse condition among the uncertainty ξ. It has been wildly applied in the mechanical loading uncertainties topology optimization problem. The key of the method is how to define the worst-case condition, which optimization is subjected to. [START_REF] Papadrakakis | TOPOLOGY OPTIMIZATION UNDER LOAD-UNCERTAINTY-FINDING THE WORST LOAD[END_REF] considered a two-player Nash game for finding the worse loading case to formulate robust optimization.

[108] applied the idea of the aggregation to contract the linear matrix composed of the local uncertain load vector and the local displacement vector, then the worst load case is easily established as the eigenvalue corresponding to the maximum eigenvalue of the matrix. Ref [START_REF] Allaire | A linearized approach to worst-case design in parametric and geometric shape optimization[END_REF] proposed the unified framework on worst-case in parametric and shape optimization under the uncertainties of geometry, elastic modulus, applied body force. A fail-safe design is proposed [115] where the damaged compliance for the worst failure case is set as the optimization objective. However, it is not easy to define the worst-case within hybrid uncertainties, and multi-physics loading.

Stochastic method is the alternative way, aims to describe a probability function of f (X, ξ), or obtain the stochastic moments of f (X, ξ). The former one is the so-called the reliability-based topology optimization, where the goal is to minimize the probability f (X, ξ) of the structure response being greater than target value f * , subjected to uncertain loading and constraints.

The latter one is nicknamed robust optimization (RTO), where the typical formulation is to minimize the weighted sum of the mean and standard deviation of f (X, ξ), which aims to simultaneously minimize the mean and standard deviation of the objective function concerning random parameters ξ. The RTO formulation is in equation 1

.8, min X f R (X) = µ(f (X, ξ)) + λσ(f (X, ξ)) (1.8)
where λ is the weight factor and µ and σ are the mean and standard deviation. Typical non-gradient (Genetic Algorithm is the most used method [START_REF] Mallipeddi | Robust controller design for systems with probabilistic uncertain parameters using multi-objective genetic programming[END_REF]) and gradient based solution are both feasible to solve the problem. The gradient-based RTO solution is fascinating because the sensitivity information efficiently guides the optimization problem toward a 'good' minimum. Although straightforward in concept, RTO poses several numerical challenges. First is the mesh-dependency and checkerboard pattern to topology optimization, successfully solved using the filter method in section 1.3.3. The second is the many query problem to qualify the uncertainty ξ propagation to the structure performance f (X, ξ) and the sensitivity information of the ∂µ(X, ξ)/∂X and ∂σ(X, ξ)/∂X. The existing papers consider the loading amplitude, loading orientation, manufacturing error, geometry shape, material properties, etc., either sole or hybrid uncertainties into robust topology optimization. The uncertainties are either random variables or random fields, which can be discretized by using different series expansion methods, i.e. Karhunen-Loeve (KL), expansion optimal linear estimation (EOLE), orthogonal series expansion method (OSE), and Polynomial Chaos Expansion (PCE) with truncation. One may refer to literature [104] on details and the comparison of different series expansion methods.

Typically, two UQ methods were the non-intrusive and intrusive ones in the literature. The wildly applied non-intrusive UQ in RTO is the stochastic collocation points method, Monte Carlo simulation(MC), quadrature points design. Besides, ref [START_REF] Garcia-Lopez | An improved robust topology optimization approach using multiobjective evolutionary algorithms[END_REF] proposed a gradient-free genetic algorithm-based optimization method with Taguchi design 'samples'. However, non-intrusive UQ turns computation unaffordable, while considering the hybrid uncertainties problem; it will cause the dimension curse issue. Some papers relieve or solve this problem through three aspects, (1) reduce the cost of the FEM simulation [START_REF] Keshavarzzadeh | Parametric topology optimization with multiresolution finite element models[END_REF], (2) use sparse grid collocation points [START_REF] Chen | Level set based robust shape and topology optimization under random field uncertainties[END_REF][START_REF] Keshavarzzadeh | Robust multiphase topology optimization accounting for manufacturing uncertainty via stochastic collocation[END_REF][START_REF] Lazarov | Topology optimization considering material and geometric uncertainties using stochastic collocation methods[END_REF] (3)uncertainty dimension reduction [132]. In addition, the intrusive approach, where analytically expressed mean and variance, is explored by literature [START_REF] Kriegesmann | Robust design optimization with design-dependent random input variables[END_REF].

Research Objectives

The goal of this dissertation is to develop a density-based multiscale topology optimization, that can consistently incorporate various multi-physics loading, microarchitectured materials, uncertainties, for solving the robust multiscale topology optimization. The specific objectives are:

• To develop an asymptotic expansion homogenization method for scale separation (homogenization at the macro scale, heterogeneous at micro-scale). It allows evaluating full structure performance by a computationally acceptable multiscale homogenization method.

• To develop an uncertainty quantification method to qualify the uncertainty to propagate to different scales. The computation cost requirement is the main issue of this topic. Combination with the homogenization model (or multiscale homogenization method) and the Polynomial Chaos Expansion method, allow the stochastic analysis of the structural performance.

Dissertation Outline

This dissertation consists of 4 Chapters. The following summary provides a brief overview of each Chapter:

• Chapter 1: Introduction. This chapter provide the background and literature survey and research objectives. It mainly introduces SIMP density-based topology's key ingredients, i.e., material interpolation, filtering technology. We will extend this TOP framework for solving our problem in the following chapters.

• Chapter 2: Homogenization. We mathematically derive the classical homogenized thermal capacity, thermal conductivity, and non-classical thermal memory function (in context of large contrast of thermal conductivity of two phases) and numerically obtain these effective material parameters with given microarchitecture. We validate the FEM obtained homogenized variables by comparing the full-scale and macro (or homogenized, with and without thermal memory effect) transient thermal simulation. The mechanical and thermo-mechanical developments are briefly recalled in the appendix.

• Chapter 3: Uncertainty propagation across scales. We take the transient thermal problem as examples, uncertainties exhibited in the micro-scale, which are the properties of the material micro gemology, to perform the stochastic analysis of heterogeneous structure by a stochastic collocation points based PCE method. They include: (1) Uncertainty propagates to the homogenized material properties: thermal capacity, thermal conductivity, thermal memory function.

(2) One saves the computation cost by applying dimensional analysis and the probability transform method. The computation cost vs. precision is checked by the brute force Monte Carlo method. (3) Proposed a PCE 2 approach to perform uncertainty propagation to the full structure with and without thermal memory effects by embedding the PCE surrogate homogenized model.

• Chapter 4 Topology optimization: We contribute to extending density based topology optimization by embedding the homogenization model. We explore two different optimization algorithms: gradient-free genetic algorithm and gradient-based optimization: (1) In gradient-free GA, we proposed a two-resolution transition scheme to reduce the design variables but maintain the FEM evaluation accuracy of the structural design performance. Moreover, we apply linear search to guarantee the volume fraction design constraints. (2) In the gradient-based optimization method, we adapted a semi-analytically PCE surrogate homogenized model to a complex architecture material, letting one directly obtain derivation of the object to micro/macro design variables.

(3) Performance of the two adopted optimization algorithms has been compared in optimized fiber-reinforced composites material structure. (4) The effect of mass and volume fraction constraints has been investigated. [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness-Application to truss structures[END_REF] In the context of the thermo-mechanical topology optimization problem, adjoint sensitivity has been derived. Thus, the gradient method has been used, diversity of local optimum designs have been explored.

• Chapter 5: Robust optimization We applied stochastic collocation points based PCE method into multi-scale topology optimization to reduce the sensitivity of design targets (the deterministic one) to the presence of uncertainty. Thus, (1) The objective of robust optimization is formulated, and the sensitivity derivation can be evaluated; naturally, the gradient method is used to solve this problem. (2) It is used to solve the mechanical design problem with the loading uncertainty presented. Then, the accuracy of PCE approximated mean and standard deviation has been checked, and the improvement of the design 'robustness' has been observed. (3) It is further applied in the complex design systems, i.e., multiphysics loading, multi-scale structure, multiple independent uncertainties (material properties, thermo-mechanical loading uncertainty), where we proposed adaptive approach (a) progressively refined mesh and (b) succeed in increasing penalty factor to improve the convergence to a feasible manufacturing design. (4) Robust optimization with hybrid uncertainties problem has been analyzed.

Nomenclature

ρ = design variable density c = structure compliance V (ρ) = volume fraction as function of density ρ E ρ = elastic modulus as function of density ρ p = penalty factor to density ρ = filtered density w ij = conic weights f R = robust objective function f (X, ξ) = structure performance as function of design variable X and uncertainty ξ µ = mean value σ = standard deviation λ = weight factor to standard deviation

Introduction

Composite materials, which consist of multi-phases with different material properties (thermal, mechanical, etc.) made with complex micro structure, often exhibit attractive properties when compared to mono-phase material. They have a wild range applications in lightweight, rigid, deployable structures. The response of composites or heterogeneous material emerges from the physics of underlying micro structure, in particular, it strongly depends on the architectures, spatial distribution and material properties of component phases. Moreover, micro morphology and material properties evolve when microstructure suffers from a macro loading.

Studying the relationship between the micro phenomenon and macro response not only provides a tool to predict the effective material properties of heterogeneous structures, but also allows one to design the material micro structure to fulfill the target performance.

However, modeling the macro heterogeneous media problem is computation costly. Homogenization is a well developed approach to extract the homogenized material properties from the heterogeneous media; it is based on the scale separation assumption, as depicted in Figure 2.1, the macro scale media 2.1 (a) consists of periodic micro structure (shown in 2.1 (b)), often the size the period is small compared to macro media, denoted with goes to zero. A single period as in 2.1(c) is the so-called unit cell, which properties are equivalent to macro scale one.

Basically, the homogenization approach can be divided into the following categories: analytical, numerical and experimental models. For simple microarchitectural material, analytical homogenization is efficient and accurate. Numerical homogenization is even universal, implemented by solving Periodic Boundary Conditions on a unit cell, and has several advantages, especially the ease of implementation and consideration of complex microarchitecture.

In this chapter, we start with solving a macro transient thermal problem where the macro properties are derived from microarchitectured material. An Asymptotic Expansion based Homogenization (AEH) is explored; further, the effective thermal conductivity and thermal capacity are obtained by first-order homogenization. In addition, a non-classical macro thermal memory physics, originating from the large contrast of thermal conductivity of two phases in the unit cell, is homogenized by higher-order homogenization process, which is used to correct the macro thermal problem to evaluate a so-called thermal delay effect. The accuracy of the homogenized macro thermal memory function has been validated by comparing it to direct heterogeneous numerical simulation.

Deterministic reference and homogenized problems

We consider herein a structure made with a micro-architectured material, i.e. that exhibits a heterogeneous micro-structure. Moreover, we consider a periodic micro-structure, deduced from a so-called unit cell. Modeling the overall problem at the length scale of the microstructure is computationally not affordable. Indeed, the scale length ratio, denoted , is small. In such cases, a multiscale modeling is mandatory. The problem is then split into a micro-scale problem defined on the unit cell only, and a macro-scale homogenized problem defined on the full structure, but that is not intended to capture the micro-scale effects, and The periodic homogenization [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] is such a tool and is used herein. It also embeds a modeling aspect since it derives explicitly the macroscopic model. We consider two different cases for a linear transient thermal problem, with a micro-structure built with two material phases. First, a case where the phases are indeed heterogeneous but with thermal capacities and thermal conductivities of same order of magnitude. This is a classical situation, where the homogenized thermal model involves the same classical form, with a homogenized thermal capacity and a homogenized thermal conductivity. Second, a case with a large contrast of conductivities (with a ratio of the order of 2 ). In this last case, the macroscopic model exhibits a memory function [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF][START_REF] Auriault | Upscaling Heterogeneous Media by Asymptotic Expansions[END_REF] and involves a macroscopic function that captures this effect. This is also a macroscopic characteristic emerging from the micro-structure.

In this section, we do not expose the methodology to find such models, but briefly recall the results obtained in these two cases. The reader is therefore suggested to refer to the literature for more developments of the homogenization approach, e.g. [START_REF] Bensoussan | Asymptotic Analysis for Periodic Structures[END_REF][START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF], for thermal conductivity, capacity, elastic modulus and thermal expansion.

The periodic homogenization relies on three main stages:

• the spatial fields are described as functions of two spatial coordinates: the slow variable x, related to the variations at large scale (typically, the scale of the considered structure), and the fast variable, y, dealing with variations at the size of the unit cell (this last variable is zoomed out with a factor 1/ ). All quantities are periodic in every direction with respect to the fast variable, they are said y-periodic;

• the primary variables are expanded with respect to parameter ; for instance the thermal field is obtained:

θ(x, y, t) = θ 0 (x, y, t) + θ 1 (x, y, t) + 2 θ 2 (x, y, t) + . . . (2.1)
• all equations of the problem are consequently expanded in powers of and are split in each power separately (hence assuming that they are verified for any value of );

• first order homogenization is considered herein.

Homogenized thermal properties 2.3.1 Classical transient thermal problem

We consider the model problem of transient conductivity first in a heterogeneous medium consisting of two components periodically distributed in the domain Ω. The transient thermal problem reads:

q = -k∇θ in Ω, (2.2)
-∇ • q = ρc θ in Ω × [0, T ], (2.3) θ = θ d on ∂Ω × [0, T ] (2.4) θ(t = 0) = 0 in Ω (2.5)
where q is the heat flux vector, k is thermal conductivity, θ is the temperature scalar field, ρ is the mass density, c is the thermal capacity. θ d is the thermal boundary on ∂Ω. ∇• is the gradient operator with respect to the spatial domain. Note that here, properties of components are in same order of magnitude, while in another case, a large contrast in thermal conductivity will be analysed in the next section. To fix the idea of homogenization, we define two coordinates, namely the slow one x, describing the phenomenon at macro scale, and the fast coordinate y, evaluating the microscopic phenomenon i.e. the response in periodic unit cell. The scale ratio is expressed as :

y = x (2.6)
Equivalently, y is the microscopic variable and x is indeed the macroscopic variable. The domain Ω is highly heterogeneous, therefore, modeling the full scale of c(x, y) and k(x, y) is computationally expensive. Usually, with intrinsic periodicity of length scale , it is assumed that Ω is homogeneous domain, finding the effective or averaged c(x), k(x) with considering its high heterogeneity in scale y is the process of homogenization. As already stated, we started homogenization by expanding the temperature filed θ(x, y, t):

θ(x, y, t) = ∞ i=0 i θ i (x, y, t) (2.7)
where each term θ i (x, y, t) is function of both macro and micro variables x, y and time t.

Following the derivation rule of:

∇• = ∇ x • + 1 ∇ y • (2.8)
where ∇ x and ∇ y are the derivative operators with respect to variables x and y, the transient thermal problem becomes a series in term of :

1 2 ∇ y • (k∇ y θ 0 ) + 1 [∇ x • (k∇ y θ 0 ) + ∇ y • (k∇ x θ 0 ) + ∇ y • (k∇ y θ 1 )] + 0 [∇ x • (k∇ x θ 0 ) + ∇ x • (k∇ y θ 1 ) + ∇ y • (k∇ x θ 1 ) + ∇ y • (k∇ y θ 2 )] + • • • = ρc θ0 + ρc θ1 + • • • .
(2.9)

Identifying each coefficient w.r.t. different orders of leads to different equations. -2 terms leading to:

0 = ∇ y (k∇ y θ 0 ) (2.10)
It implies that θ 0 is the macroscopic temperature, independent of y coordinate.

The -1 terms is:

0 = ∇ x [k∇ y θ 0 ] + ∇ y [k∇ x θ 0 + k∇ y θ 1 ] (2.11)
This leads to a steady-state thermal problem on cell Y , allowing one to solve the y-periodic unknown θ 1 . It is linear w.r.t. the ∇ x θ 0 . Normally, it was solved by FEM with periodic boundary conditions on unit cell. As detailed in following section. We obtain:

∇ y θ 1 = -L y (∇ x θ 0 ) in Y (2.12)
where L y is a linear operator. the 0 0-order problem leads to the macroscopic transient thermal problem on Ω:

ρc θ0 = ∇ x [k∇ x θ 0 + k∇ y θ 1 ] + ∇ y [k∇ x θ 1 + k∇ y θ 2 ] inΩ θ 0 = θ d on ∂Ω θ 0 (t = 0) = 0 (2.13)
to get rid of the unknown θ 2 in Y domain, since ∇ y θ 1 = -L y (∇ x θ 0 ), one applies average operator<> Y on Y domain, to transform the equation into:

ρc Y θ0 = ∇ x [ k(1 -L y ) Y (∇ x θ 0 )] in Ω × [0, T ] (2.14)
With the homogenization procedure, one may solve a macroscopic transient thermal problem θ 0 with effective propeties ρc M = ρc Y and k M = k(1 -L y ) Y . Full solution temperature can be recovered approximately as θ ≈ θ 0 + θ 1 , namely with a first -order multiscale homogenization method. The numerical homogenization cases are implemented and validated in next section.

Non-classical thermal memory function

In the context of large contrast in thermal conductivity of components, a non-classical thermal memory effect arises. We expect to directly solve the thermal memory delay problem by homogenization procedure. Recalling the same definition of x, y, ∇, Y , Ω and , one may refer to the previous section, and subscript F and S are indicated as different material phases. It is assumed that the repeated unit cell Y domain is composed of two phases S and F , occupying the domains Y S and Y F ; Y Γ = Y F Y S is the interphase. Large contrast in thermal conductivity arises when k F ≈ 2 k S [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF], their capacities are of same order i.e. ρ F c F ≈ ρ S c S . Thermal memory effect is expected to be interpreted with the characteristic time, since macro characteristic time is

τ M = (c M /k M )(L/π) 2 (2.15)
where the c M and k M are the homogenized thermal capacity and conductivity, and L is the characteristic length of macro scale. The characteristic times of S and F phases are τ S = (ρ S c S /k S )(l S /π) 2 and

τ F = ρ F c F l 2 f k F π 2 ≈ ρ F c F L 2 k F / 2 π 2 ≈ τ S / 2 ≈ τ M (2.16)
Here introduces an intermediate thermal conductivity k F = k F / , and characteristic length corresponding to different scales are: l S ≈ l F ≈ l and L S = l S / ≈ L F = l f / ≈ L, where L and l are associated to macro and micro length, respectively. In such case, the transient thermal response of F phase may arise from the macro scale. A phenomena of highly thermal heterogeneity on the unit cell appears, one can capture this phenomena by relocalization, or a multiscale method. However, we expect to directly capture the so-called thermal memory effect by a homogenization procedure. For this purpose, the thermal equations then become:

ρ F c F θF = ∇(k F ∇θ F ) inY F × [0, T ],
(2.17)

ρ S c S θS = ∇(k S ∇θ S ) inY S × [0, T ], (2.18) θ F = θ S onY Γ , (2.19) k S ∇θ S n S = -k F 2 ∇θ F n F onY Γ . (2.20)
From the beginning, Equation 2.18 to solve the θ S at order -2 implies that, θ S0 is of macroscopic temperature field, does not rely on fast coordinate y.

At order of -1 , it gives:

∇ y (k S ∇ x θ S0 ) = -∇ y (k S ∇ y θ S1 ) in Y S (2.21) k S ∇ • θ S n S = 0 on Y Γ (2.22)
It is a steady state thermal problem on S, one can recall the first order solution θ S1 (or characteristic of S-microstructure) using linear operator in Equation 2.12 by:

∇ y θ S1 = -L Sy (∇ x θ S0 ) on Y S (2.23)
The numerical determination of L Sy requires: (1) as many resolutions of the previous linear micro problem on Y S as there are independant components in ∇ x θ S0 , i.e. 2 for 2D problem 3 for 3 D problem; (2) the storage of the same number of temperature gradient fields on Y S [START_REF] Dureisseix | Numerization of a memory effect for an homogenized composite material with a large contrast in the phase thermal conductivities[END_REF]. Indeed, It allows to recover the heat flux q S by:

q S = -k S (1 -L Sy )∇ x θ S0 (2.24)
Then the boundary condition of Equation 2.20 at -1 and 0 , respectively, are expressed as:

k S ∇ y θ S0 n S = 0 on Y Γ (2.25) k S (∇ y θ S1 + ∇ x θ S0 )n S = 0 on Y Γ (2.26)
which indicates that θ 0 on interphase Y Γ is independent to either x or y. Then, the order of 0 for θ F provides the following problem:

ρ F c F θF = ∇ y k F ∇ y θ F 0 in Y F × [0, T ] θ F 0 = θ S0 on Y Γ × [0, T ] (2.27)
It implies that θ F 0 only y dependent, otherwise no admitted solution to θ F 0 . Due to the θ S0 (x, t) as a boundary condition on interphase, one introduces θ F 0 (x, y, t) = θ S0 (x, t) + w(x, y, t) into Equation 2.27, generates a new problem

ρ F c F θS0 + ρ F c F ẇ = ∇ y (k F ∇ y w) in Y Γ × [0, T ] (2.28)
where w is y periodic and null on Y Γ . The order 0 of Equation 2.18, and order 1 of Equation 2.20(boundary condition), yield to the problem:

ρ S c S θS0 = ∇ x [k S ∇ x θ S0 + k S ∇ y θ S1 ] + ∇ y [k S ∇ x θ S1 + k S ∇ y θ S2 ] inΩ k S (∇ x θ S1 + ∇ y θ S2 ) = -k F ∇ y θ F 0 onY Γ (2.29)
One defines the partial averaging operators as:

• S = 1 vol(Y ) Y S •dY and • F = 1 vol(Y ) Y F •dY , so that • y = 1 vol(Y ) Y S •dY = • S + • F Combining
ρ S c S S θS0 = ∇ x ( k S (1 -L Sy ) S ∇ x θ S0 ) -ρ F c F θF0 F (2.30)
where, the macro parameters of are ρ S c S S = C M S and k S (1 -L Sy ) S = k M S , respectively. Once discretized finite elmements, the problem involves macro thermal problem with an additional term of θF0 on a full micro cell, i.e. a FE 2 procedure [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials[END_REF] can be used. For the efficiency and storage saving purposes, memory effect model [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF] is proposed to replace the term of ρ F c F θF0 F , turning Equation 2.30 into:

c M S θS0 = ∇ x (k M S ∇ x θ S0 ) + τ =t τ =-∞ θS0 (τ )β(t -τ )dτ (2.31)
where β(t) is a macro memory function at the macro scale [START_REF] Auriault | Effective macroscopic description for heat conduction in periodic composites[END_REF], vanishing in time t and independent of microscale boundary value problem. Considering the practical intial conditions, one defines θS0 (t ≤ 0) = 0, θF0 (t ≤ 0) = 0, θ S0 (t ≤ 0) = 0 and θ F 0 (t ≤ 0) = 0.

To achieve the micro dependent θ F 0 , a transient problem Equation 2.27 with a time dependent thermal boundary θ F 0 (t) = θ SM (t) on inter phase between two phases requires to be solved. Therein, Duhamel's theorem is adopted, instead of solving 2.27, one can first solve a auxiliary problem u(y, t) with constant boundary,

ρ F c F u = ∇ y k F ∇ y u inY F × (0, T ) u(t = 0) = θ c1 inY F u(t) = θ c2 onY Γ × (0, T ) (2.32)
where u(y, t) is the temperature disbribution inside the F phase, the θ c1 is any non-zero thermal field, and θ c2 = 0 is null on interphase, then relates the solution of θ F 0 (t) to auxiliary one u by Duhamel intergral:

w = τ =t τ =-∞ u(y, t -τ ) θS0 (τ ) d τ (2.33)
This memory function, obtained as:

β M (t) =< C F u(y, t) > F (2.34)
The accuracy of β achieve from asymptotic homogenization method is depended on the resolution of spatial and temporal discretizations.

Physically, the macro thermal memory effect is due to the large contrast of thermal conductivities on different phases material in unit cell. Mathematically, it originates from k F /k S = 2 but ρ F ≈ ρ S , can be interpreted by the characteristic time, that is τ F ≈ τ M ≈ τ S / 2 . Asymptotic expansion theory is applied, expanding the temperatures θ F and θ S w.r.t. orders of . the memory function is derivative from the solution of globally thermal equations(equation 2.17 -2.20). Summarily, one can concluded:

• θ S0 is the macro temperature.

• θ S1 is the microscopic temperature, ∇ y θ S1 can be recalling from ∇ x θ S0 and liner operator L Sy .

• θ F 0 is y-periodic. can be expressed as θ F 0 = θ S0 + w.

• By averaging, the macro thermal problem is derivatived (equation 2.30), where the difference compare to the classic thermal problem relies on: θ S0 is the quantity of interest; homogenized thermal capacity and conductivity are averaged on S phase; an additional component ρ F c F θF0 is introduced to correct such problem.

• To the end of storage saving, ρ F c F θF0 is expressed by a memory effect model

τ =t τ =-∞ θS0 (τ )β(t -τ )dτ
. β is indeed the macro thermal memory function.

• Submitting θ F 0 = θ S0 + w into equation 2.25, using Duhamel principle (equation 2.32), w =

τ =t τ =-∞ u(y, t -τ ) θS0 (τ ) d τ . Consequently, θ F 0 is solved. • By averaging, finally β M =< C F u(y, t) > F .

Numerical implementation

The thermal problem without memory effect is indeed the classical transient thermal problem in heterogeneous periodic media as one discussed in equation 2.14. Instead of modeling the full scale problem, we model a macro-scopic problem with homogenized thermal capacity and conductivity on an underlying unit cell. The transient thermal problem with memory effect is discussed in section 2.3.2, it is due to the large contrast in thermal conductivity of phases, and will be solved by a macro approach, where an additional macro thermal memory β M (t) are applied. Consequently, three homogenized C M k M β M (t) should be calculated offline, each has its own features leading to different strategies. Details on numerical implementations of homogenized c, k and β are illustrated. The proposed macro thermal memory solution is validated by a reference full scale modeling.

Macro thermal capacity C M

Analytical homogenization of thermal capacity is made by averaging:

c M = ρc y (2.35)
The computation cost is negligible, which resided in measuring the volume fraction of different phase on unit cell. On a structure consisting of S and F phases,

c M = n F (ρ F c F ) + (1 -n F )(ρ S c S )
, where n F is the volume fraction of F phase.

Macro thermal conductivity k M

The mathematically derived thermal conductivity based on AEH requires a localization, which is a steady state solution on the unit cell with periodic boundary conditions (PBCs). It indicates that k M is not only phases properties dependent, but also the micro architecture dependent thermal property. Meanwhile, with applying FEM method, the main procedures of homogenization technology are: 1. Determine the unit cell structure. The macro structure is assumed to be sufficiently homogeneous, i.e., in mathematics view, ≈ 0, at the same time, the length scale of unit cell should be smaller than the characteristic length of the macroscopic sample. The assumption suggests that, the macro scale structure consists of spatially repeated unit cells, thus periodicity at unit cell level is proposed.

FEM solution

Once determines the region of unit cell, one spatially discretizes the unit cell depicted in figure 2.3, with inclusion F and matrix S. Given the periodic unit cell Ω Y , the boundaries are decoposited into two opposite groups ∂Ω Y + and ∂Ω Y -, satisfying the constraints of

∂ω Y = ∂ω Y + ∪ ∂ω Y -and 0 = ∂ω Y + ∩ ∂ω Y -Perodicity Boundary Conditions (PBCs) are given as, θ k+ = θ k- (2.36)
k+ and k+ are subscripts indicating the periodic boundaries on the unit cell Ω Y . In practice, PBCs are introduced to the pair multi-points on opposite surfaces. Thus, one of the prerequisites to categorise the pair nodes on opposite surface into 3 sets, i.e. vertex, inner-edge and inner face. We only categorise the nodes on one face as an example, as shown in figure 2.2.

Figure 2.2:

The categorized sets for applied periodic boundary condition

• set-1 (vertex):

• set-2 (inner edges):

• set-3 (inner surfaces):

The pair-nodes in set-2 belong to edges but excluding the pair nodes on set-1. The pair-nodes on set-3 are the points belonging to the surfaces excluding the ones on set-2 and set-3. Note that there are no conjoint nodes of 3 sets to avoid over constrained finite element analysis system. By solving the PBCs on unit cell, then take averaging operator, K M is calculated:

kZ 1 = -kZ 0 with PBCs (2.37) Z 1 = LZ 0 (2.38)
where the Z 0 is the given unity gradient with respect to macro scale, and Z 1 is the calculated gradient temperature with respect to unit cell level. Once solving the 2.37, one can construct the so called linear operator L. To determine the different components of and3) it has been numerically simulated three times with applied different unity gradient thermal fields of Z i 0 , i = 1, 2and3, construct linear operators corresponding to different components L ij . Then take the averaging operator:

k ij M (k i1 M , k i2 M and k i3 M ) with i = 1, 2
k ij M = Ω Y k(1 -L ij ) d Ω Y (2.39)
A numerical case is proposed, a porous cell structure with two phases( in Figure 2.3), in which, the inclusion F is surrounded by matrix S phase, their dimensionless thermal capacity and conductivity are assumed to be (ρc) s = (ρc) f = 2.0 × 10 6 , k s = 20, k f = 10, respectively. It is assumed a perfect inter phase of inclusion and matrix that heat flux can totally transfer, effective thermal conductivity of the composites homogenized in k M : Other types of PBCs on unit cell are feasible, the only requirement for them are fulling the periodicity. Note that, this homogenization assumption relies on the concept of scale separation (heterogenous in micro scale, homogenous in macro scale) is not appropriate in highly heterogeneous macro scale region. • One applies the same materials as used in homogenized thermal conductivity problem.

• It is implemented on micro level F phase. To simplify, only 1/8 of the cell can be simulated, periodicity conditions that are automatically satisfied due to the planar symmetries.

• Expand the size of physical unit cell by ratio of . The characteristic of micro unit cell and macro structure are of the same order.

• Material properties c F and k F are applied; note that k F = k F / 2 , is not a physical thermal conductivity.

• Define initial and boundary conditions (in equation 2.32), θ c2 is null on the interphases between S and F materials. θ c1 = 1, a unity initial temperature is adopted.

• Time steps nt and duration t ≈ 3τ F , also with uniform time step t/nt. To satisfy instabilities maximum discrete principle, a thumb rule [114] to determine the t/nt is written:

t/nt = n × τ F /nt ≥ (ρ F c F e 2 )/(6k F ) (2.41)
where e is the element size.

• In each iteration, one obtained β M (t) by averaging operator in equation 2.34. With applied nt = 45, one achieves a good convergent β M compared to coarse time discretization. Since the accuracy of this proposed numerical solution to solve β M has been validated compared to a analytically solution [START_REF] Dureisseix | Numerization of a memory effect for an homogenized composite material with a large contrast in the phase thermal conductivities[END_REF]. In next section, we will validate it from a different perspective: comparing the solution of equation 2.31 to the solution obtained by full-scale modeling.

Macro thermal memory problem, correlation θ SM

In the context of large contrast in thermal conductivities, the homogenized model exhibits a memory function β M (t), to estimate the temperature evolution in fast S phase θ SM (x, t).

To assess the validity of this homogenized model, in this section we compare its prediction to the results provided by the full scale model as a reference, and to the results of the macroscopic model without memory function, for an affordable deterministic problem, with the boundary conditions of Figure 2.6,

θ d (t) = θ a + θ a sin 2π t T c - π 2 (2.42)
with θ a and T c being given amplitude and characteristic loading time. For this deterministic test case, the geometric parameter (volume fraction of F phase) is n = 0.4, and material properties are C S = C F = 1, k F = k S = 2 and k F = 10 -4 × k F , so = 100 herein. The macro thermal characteristic C M , k M and β M (t) are computed offline as described in the previous section. For this macroscopic model, the thermal characteristic time is

𝜃 𝑥 1 =0 = 𝜃 𝑑 (𝑡) 𝜃 𝑥 1 =𝐿 = 0 𝜃 𝑥 1 =0 = 𝜃 𝑑 (𝑡) 𝜃 𝑥 1 =𝐿 = 0 𝜃 𝑥 3 =0 -𝜃 𝑥 3 =𝑙 = 0 𝜃 𝑥 2 =0 -𝜃 𝑥 2 =𝑙 = 0 𝜃 𝑥 2 =0 -𝜃 𝑥 2 =𝑙 = 0 𝑥 1 𝑥 2 𝑙
τ M = C M L 2 /k M .
Figure 2.7 (a)-(c) reports the macroscopic temperature evolution θ SM /θ a with respect to time scale t/T c , for the macroscopic problems with and without memory function, at central location x 1 = L/2 as well as the reference temperature, obtained by averaging on the central unit cell fast phase S of the full reference temperature obtained at microscale. Quasi-identical results are obtained for the reference and the macroscopic problem with memory function, while some discrepancies occur both in amplitude and time shift for the macroscopic problem without memory function, hence emphasizing the interest of using the dedicated memory function as a macroscopic property.

We play a little trick here by reversing the inclusion and matrix material properties in the unit cell. In that situation, it requires to derive the homogenized β M on the matrix phase with a different geometry (labeled as Geo1 in figure 2.7 (d)). In the figure 2.7 (d), the green curve is exactly the case shown in figure 2.7 (b), the blue curve is the new case. It shows that the geometry has a significant influence on the homogenized β M and the macro thermal problem.

Conclusions

This chapter introduces applying the Asymptotic Expansion Homogenization to solve the effective thermal properties of the underlying heterogenous unit cell, which is repeatedly distributed on the macro scale. A non-classical thermal characteristic, named thermal memory function arising from larger contrast conductivities of two phases is demonstrated in this chapter. We have mathematically derived the solutions of macro thermal capacity, conductivity, and memory function by AEH, proposed distinct numerical solutions to different macro properties by FEM.

We implement three macro transient thermal cases, being the direct numerical case on heterogeneous structure, the macro case with effective capacity and conductivity, the macro case with thermal memory function, to validate the accuracy of our proposed thermal memory function to describe the thermal lag phenomenon.

As benefits from the spatial homogenization technology, we are able to analyze the complex heterogeneous structure performance using an affordable computational approach. Moreover, it is the fundamental development for the further uncertainty qualification and optimization problem.

Introduction

Materials with engineered micro-structures are still an emerging topic, since they can be optimized for given functionalities, and therefore may exhibit higher performances when compared to standard bulk materials, [START_REF] Ashby | Designing hybrid materials[END_REF][START_REF] Fleck | Micro-architectured materials: past, present and future[END_REF]. Nevertheless, their manufacturing may also be costly, and their recycle-ability an issue. They can therefore be targeted not for mass production, but tuned for on-demand applications, and eventually for a dedicated structural part, or a dedicated critical location in a mechanical part.

As a counterpart, there are some uncertainties involved, due to a not-so-well controlled production flow, that should be taken into account for the overall design problem. For optimization stage, this leads to robust optimization issue, [START_REF] Beyer | Robust optimization-a comprehensive survey[END_REF]. For all design phases, this leads to the uncertainty quantification issue, [119,129].

In this chapter, we will focus on uncertainty propagation of virtual material prototyping, issuing the question of scale propagation of uncertainties due to the micro-scale material behavior, and geometry of the micro-structure. Two different cases built upon the previous chapter are concerned: (i) a classical macro case of a limited contrast in thermal conductivities of the material phases, and (ii) the large contrast case, leading to different characteristic transient thermal times, and a nonstandard macroscopic behavior that is also studied for uncertainty propagation throughout scales.

Since we promote the non-intrusive approach, no stochastic finite element method is used but a deterministic parametric FE analysis, together with a polynomial chaos expansion (PCE) approach with collocation both at different spatial scale. Generally, PCE showed better computational efficiency, when compared with standard and modified Monte Carlo techniques.

Review of uncertainty propagation and quantification

A deterministic model owns no capability of dealing with the engineering problem with uncertainty variables. However, a stochastic model with analytical expression is a kind of challenge in practical application. For this purpose, Uncertainty Quantification (UQ) is proposed to obtain a relative approximation to original model. The general concept of uncertainty propagation is to consider the code between random input and output variables, read as:

Y = M(X) (3.1)
In this context, Equation 3.1 described random Y due to propagation of random inputs X through a numerical model M. With adequate realizations of M(X), the probability density function (pdf) of Y can be inferred, some quantities of interest (e.g. mean, variance) are obtained.

The principle objectives of Uncertainty Quantification and Propagation are listed as follow:

• Validation. Simulation model must be validated in context of real engineering application, where exist multiple experimental and numerical uncertain ranges.

• Variance analysis. The variance of the response around its mean value characterizes the robustness of the numerical prediction and the system's controllability. They (mean and variance) provide essential information relevant to the design or the optimization of the system, taking into consideration uncertainties.

• Risk analysis. Based on the probability laws of the input data, it is often desired to get the probabilities of the system response within a threshold to guarantee the health of the system. Such probability low described outputs that can be used to risk or reliability assessment analysis .

• Optimization considering uncertainties. By providing mean and variance or probability law of the outputs data, it is possible to manage the system to be less sensitive to the uncertainty by optimization. A so-called Robust or Reliability optimization is proposed.

To these purposes, some classical UQ methods are reviewed, including Monte Carlo, spectral method.

Monte Carlo: a brute force non-intrusive method

The Monte Carlo (MC) [100] is a sampling method of realizing the model in equation 3.1 with quite a large number of random sample points. The general procedure of MC is summarized and can be easily implemented by:

• Design of experiment through sampling method: The MC is kind of non-intrusive method that takes the M as a black box with inputs of n samples, i.e,

X = {x i , i = 1, 2, ..., n}
. Random or quasi-random sampling method (e.g. sobol sequence, latin hypercube and Niederreiter sequence etc.) are adopted.

• Compute the system response y i = M(x i ), (i = 1, 2, ..., n).

• Evaluate the mean µ y and standard deviation σ y by:

µ y = 1 n M i=1 y i (3.2) σ y = 1 n 2 n i=1 y 2 i -µ 2 y (3.3)
• Fitting the Probability Density Function (PDF) f Y (y). Since f Y (y) is an unknown distribution, non-parametric fitting method such as kernel density can be used to estimate the PDF. Subsequently, Cumulative Distribution Function (CDF) of y is obtained by integrating f Y (y).

The MC is a non-intrusive approach (no need to rewrite the original or deterministic model), for building such a case is effortless. Unfortunately, the rate of convergence of MC is 1/ √ n, which is relatively slow, however independent of the dimension of the uncertainty vector. The computation cost mainly relies on the number of samples n; thus an appropriate sampling method, exploring through the whole sampling space with fewer points, may help relatively improve the efficiency of MC.

Spectral Methods

As outlined above, MC is the non-intrusive stochastic collocation method, the uncertainty must be sampled with an adequately fine resolution to evaluate the variation of outputs. In contrast, the spectral method is based on a completely different approach, namely reconstruction of the equation 3.1 expanded such as:

Y = +∞ j=0 y j Ψ j (X) (3.4)
in which X is the random vector defining the random space, Ψ j are the orthogonal basis function, and y j are the deterministic unknown coefficients. Once the coefficients are determined, the statistics of Y can be immediately exploited. Without too much details, we explore a short list literature on the development and application of the spectral method in the UQ. Reference [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] is wildly recognized as the original work for spreading the spectral method in UQ. This early work starts from solving the linear mechanical problem with spatial varying elastic material properties and loading uncertainties [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Ghanem | Spectral stochastic finite-element formulation for reliability analysis[END_REF]. Afterwards, it is developed to deal with UQ in physical systems of fluid dynamic [START_REF] Le Maître | Spectral methods for uncertainty quantification: with applications to computational fluid dynamics[END_REF], bio-mechanics [101], heat transfer problem [START_REF] Carnevale | Uncertainty quantification: A stochastic method for heat transfer prediction using LES[END_REF] and topology optimization [111]. However, the essential ingredients to the spectral method rely on the determination of the coefficients y j to their corresponding basis function Ψ j .

In the following, we introduce a polynomial chaos expansion base spectral method. The selected orthogonal polynomials basis is used to construct the basis function according to the PDF to the uncertainty inputs. Typically, two approaches, non-intrusive and intrusive solutions, to obtain coefficients y j are briefly introduced.

Polynomial Chaos Expansion (PCE) based Spectral Methods

The PCE is one of the spectral representation methods, the random output Y is expressed as an infinite series expansion in equation 3.4, where, {Ψ j } +∞ j=0 forms basis of Hilbertian space H ⊃ Y , {y i } +∞ j=0 is the set of coordinates of Y with respect to this basis. Hilbertian analysis guarantees the existence of such bases and representation, however many choices are possible.

The equation 3.4 can be developed with different series expansion methods, such as, Karhunen-Loève (K-L) expansion, Orthogonal Series Expansion (OSE), and Expansion Optimal Linear Estimation (EOLE) method, [104]. Herein, we consider Polynomial Chaos Ex-pansion [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]110] as an effective tool for constructing Hilbertian basis, {Ψ j } inf j=0 , in which the basic terms are multivariable orthogonal polynomials in the inputs vector X, and y j are the coefficients that can be determined accordingly. Such a tool is called PCE stochastic model [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF] which is qualifying the uncertainties on the Quantity of Interest (QoIs) and obtaining surrogate model to original numerical model within acceptable accuracy.

Key ingredients of PCE-based spectral method

The key ingredients to develop such stochastic model in equation 3.4 rely on: (1) truncate after a finite series term, (2) decide the orthogonal function Ψ j , (3) determine the coefficients y i

• Univariate orthogonal polynomials [103]

The multivariable orthogonal polynomials basis Ψ j is built startung from a set of univariate orthogonal polynomials ψ

(i) k (x i ).
For the sake of simplicity, one assumes the random input vector has independent components denoted as {X i , i = 1, ..., M }, therefore the joint distribution is given by a simple product:

f X (x) = M i=1 f X i (x i ) , x i ∈ D X i (3.5)
Where the D X i is the support of X i , for each X i and given two functions ψ 1 , ψ 2 : x ∈ D X i , provided an exist function:

ψ 1 , ψ 2 i = D X i ψ 1 (x)ψ 2 (x)f X i (x)dx (3.6)
The right hand side is the expectation E [ψ 1 (X i ) ψ 2 (X i )] with respect to the f X i . When evaluated to be 0, one determines that such two functions are orthogonal with respect to the probability P(dx) = f X i (x)dx. Using this notion and applying the Gram-Schmidt orthogonalization procedure, one allows to build a family of orthogonal polynomials. For some parametric pdf distributions, table 3.1 reports the well-known families of orthogonal polynomials. For example of uniform distribution, the Legendre Polynomials up to 4th order are depicted in figure 3.1.

• multivariable orthogonal polynomials basis Ψ α [103]:

One intends to construct the basis in equation 3.4 by the tensor product of univariate orthogonal polynomials. From the start, one defines the multi-index α = (α 1 , α , ...α M ), which are corresponding to multivariate polynomial Ψ α : k are defined according to the i-th marginal distribution. One admits Ψ α , α = (α 1 , α , ...α M ) in the input vector X are orthonormal:

Ψ α (x) def = M i=1 ψ (i) α i (x i ) (3.7)
E [Ψ α (X)Ψ β (X)] def = D X Ψ α (x)Ψ β (x)f X (x)dx (3.8)
if α = β, the right hand side is equal to 1, otherwise 0. For ease of implementation, the equation is truncated to a certain term, let one defines the total degree of Ψ α by:

|α| def = M i=1 α i (3.9)
where the M is the truncation threshold. One standard truncation approach is selecting all the polynomials that |α| is smaller than a provide value p, finally, the number of terms is [START_REF] Huang | Collocation-based stochastic finite element analysis for random field problems[END_REF]127]:

M = N + p p = (N + p)! N !p! (3.10)
where N is the uncertainty dimension.

In such context, one builds the PCE stochastic model upon the outline above, expressed as:

Y = M α=0 y α Ψ α (3.11)
The calculation of the polynomial coefficients for uncertain parameters uses various procedures, e.g., Galerkin projection, collocation method, and moment method.

Intrusive solution: stochastic finite element method

We briefly described the intrusive solution to determine the coefficients exemplified by the elastic mechanical problem with material uncertainties. The general details may refer to [START_REF] Herzog | From Nano to Space[END_REF]. In the case the deterministic linear equation system is described as:

KU = F (3.12)
where K is the stiffness matrix, U is the displacement solution vector and the F is the force loading vector. Assumed the force loading is deterministic and material properties are random, it leads the two stochastic representations of K(X) and U (X), to be expanded by equation 3.11:

K(X) ≈ M j=0 K j Ψ j (X) U (X) ≈ M i=0 U i Ψ i (X) (3.13)
where the K j and U i are the corresponding coefficients of the selected basis function Ψ j and Ψ i , respectively.

Submit them into the deterministic elastic mechanical equation 3.12, consequently leads to:

M i=0 M j=0 K j U i Ψ i Ψ j ≈ F (3.14)
Next, one performs a Galerkin projection onto the orthogonal basis polynomials of the Chaos space:

M i=0 M j=0 K j U i Ψ i Ψ j Ψ k = F Ψ k = F δ 0k k = 0, 1, . . . , M (3.15)
The stochastic linear system can be rewritten as:

   M j=0 K j Ψ 0 Ψ j Ψ 0 • • • M j=0 K j Ψ M Ψ j Ψ 0 . . . . . . M j=0 K j Ψ 0 Ψ j Ψ M • • • M j=0 K j Ψ M Ψ j Ψ M    Kpc    U 0 . . . U M    Upc =      F 0 . . . 0      Fpc (3.16)
Building the K pc is also called the matrices second level of assembling. One may refer to references [START_REF] Herzog | From Nano to Space[END_REF]104] for the detail on assembling the matrix K pc . The vector U pc is the obtained PCE coefficient, in which the unknown dimension is N ×M . The M is the size of polynomial basis, and N is the degree of freedom of the structure (e.g. twice the node number for twodimension elastic-mechanics problem). It obviously shows the disadvantages of the intrusive method in the perspectives of: (1) high memory usage; (2) require rewriting/adaptation of source code to the deterministic model; (3) difficult to implement for some interesting values (e.g. von Mises stress) [START_REF] Herzog | From Nano to Space[END_REF]. However, the reduced-order computational model approach is the alternative solution to overcome the high memory usage issue; it has been applied in solving the stochastic linear dynamic problem [START_REF] Capillon | Uncertainty quantification in computational linear structural dynamics for viscoelastic composite structures[END_REF].

Non-intrusive solution: stochastic collocation points

With the non-intrusive solution, we are not interested in a specific deterministic numerical model, for we take the deterministic model as back box model. Different samplers, e.g. MC LHS and QMC (using Halton's sequence) can be applied. However, in the PCE stochastic model context, we used stochastic collocation points strategy based on the root of polynomial basis [START_REF] Huang | Collocation-based stochastic finite element analysis for random field problems[END_REF][START_REF] Sharma | Acoustic performance of a periodically voided viscoelastic medium with uncertainty in design parameters[END_REF]127], which significantly saves the computation requirement of the stochastic process. The key procedures of PCE-based stochastic collocation points are described as follows:

• Design of experiments X (stochastic points method):

The design of experiments is the key ingredient; it is fundamental that the experimental points set pursue a good cubature rule; it is simple for a single parameter problem (N = 1), an optimal choice is the Gaussian Quadratures [125]. In addition, the classical rule-of-thumb for selecting 1D points is to consider the roots of the polynomial of order p + 1 [START_REF] Berveiller | Stochastic finite element: a non intrusive approach by regression[END_REF]. However, it is challenging in higher dimension space (N >> 1), One choice is to use the tensor product of the 1-D points, e.g., Gaussian Quadratures or roots of the polynomial. The approach has been validated with an acceptable error. However, if one uses q points in the 1D problem, it turns q N in N dimension problem. The problem is that the number of collocation points increases rapidly in high dimension space.

To overcome the dimension increase issue, sparse grids constructed by the Smolyak algorithm were proposed in high order stochastic collocation points design [126]. Fig. 3.2 shows the 2D grids using a different strategy based on the same 1D nodes. It is clear that full tensor grids consist of much more nodes than sparse grids; since then, one can design higher dimension uncertainty space efficiently. • Determine coefficients (or coordinates) of basis y α :

Once the truncated basis has been selected, the coefficients y α can be calculated. Introducing Ŷ = {y 0 , . . . , y P -1 } and Ψ(x) = {Ψ 0 (x), . . . , Ψ P -1 (x)} Minimizing the mean-square error of the residue ε P leads to determine the unknown coefficients Y:

Ŷ = arg min Y∈R P E M(X) -Y Ψ(X) 2 (3.17)
• Post-processing Due to the orthogonality of the PCE basis, one can explicitly gives the mean and standard deviation without sampling, they read:

µ PC Y = y 0 (3.18) σ PC Y = M -1 α=1 y 2 α Ψ α 2 (3.19)

Comments

Similar to stochastic Galerkin methods (intrusive methods), the collocation methods(nonintrusive methods) take advantage of an assumption of smoothness of the solution in random space to achieve fast convergence. However, the numerical implementation of PCE based stochastic collocation is trivial, as it requires only repetitive runs of an existing deterministic solver, similar to Monte Carlo (MC) methods. However, when one introduces a large number of uncertainties into a deterministic system, the problem called dimension disaster arises due to the number of collocations points being exponential to the uncertainty dimension to the power of PCE order. The intrusive PCE-based stochastic Galerkin method is a memory storage-consuming method; it requires to rewrite and adapt the original deterministic model. However, compared to the non-intrusive method, it still has intense competition in some particular aspects:

(1) it achieves fast convergence compared to MC; (2) the size of a linear system of equations grows linearly with the number of basis polynomials, can be competitive compared to exponentially increasing the dimension system of the non-intrusive method (a dimension disaster for the non-intrusive method); (3) better handling of the transient approach compared to non-intrusive method (temporal discretization is necessary for a non-intrusive method before PC expansion).

However for our multiscale problem (in chapter 2, a transient thermal problem via heterogeneous medium with and without thermal memory effect), we proposed a scale separation approach to solve the macro problem at coarse mesh with applied homogenized thermal capacity, conductivity, and thermal memory function (obtained offline at RVE refined mesh mesh). Modifying these equations will lead to a complexity to our models. Therefore, we will use non-intrusive method.

In the following sections (and chapters), the QoIs (effective material properties, structure performance, optimization design targets) will vary in each presented numerical simulation system. Therefore, developing the corresponding non-intrusive stochastic model is elaborate. By using PCE-stochastic points method for its non-intrusive implementation way lead computation acceptable, and it will be validated by MC approach.

D X Distribution PDF f X ≡ w Family [-1,1] Uniform 1/2 Legendre R Gaussian e -x 2 /2 / √ 2π Hermite R + Exponential e -x Laguerre [-1,1] Beta (1-x) α (1+x) β B(α+1,β+1)
Jacobi Table 3.1: Families of classical orthogonal polynomials

Numerical implementation: micro uncertainties propagate through different scales

The deterministic homogenization model has been fully understood in chapter one; further, the uncertainties of material properties and geometry at micro-scale propagation through macro scale will be studied using stochastic method with PCE and also will be validated using brute force MC method. The section presents a computational framework for quantifying the uncertainties in microscale propagation to a full-scale transient thermal problem with and without thermal memory effects. The memory effect is derived from large contrast in the thermal conductivity of two phases. We combine homogenization model and stochastic collocation method to reduce the computation cost. First, 3 macro thermal properties, thermal capacity, conductivity, and memory function, are evaluated via the simple analytical expression, FEM-based homogenization on unit cell, and transient thermal simulation on low-conductivity phase, respectively. In section 2.4.4, we have proved the accuracy of our macro model with memory function to capture the thermal delay phenomenon. It is the deterministic model combined with three offline homogenized models (homogenized thermal capacity, conductivity, and memory function) we consider here. Secondly, the uncertainties of material properties and microarchitecture propagation through the macro thermal characteristics and have been studied using Polynomial Chaos Expansion (PCE) method. We apply a dimension reduction analysis and Probability Transform Method (PTM) to reduce the number of uncertainties (overcoming the dimension disaster issue), efficiently construct the PCE-surrogate model. Lastly, we proposed a PCE 2 approach to provide a globally efficient predictive model to assess the micro uncertainties propagate to the full-macroscopic transient thermal model with and without thermal memory effect.

Macroscopic thermal capacity

When the architectured material is made with two homogeneous phases, obtaining the macroscopic thermal capacity as a quantity of interest (QoI) relies on a simple analytical expression of the microscopic input parameters as

C M = f (C F , C S , n) = nC F + (1 -n)C S .
As a stochastic model of the input, we consider herein the material parameters and the geometric parameter n as uniform independent random variables, whose characteristics are given in Table 3.2.

To predict and quantify the macroscopic capacity, at least two approaches can be considered:

• obtain the stochastic model of the output, i.e. the QoI C M , as the sum of two products • derive the empirical cumulative density function (ecdf) numerically using a sampling technique.

This last approach, as more general, is selected herein: from a sampling of the 3 concerned inputs (with 400 samples obtained with a Sobol sequence [START_REF] Sobol | Construction and comparison of high-dimensional Sobol generators[END_REF]), the set of computed outputs give a plot of the ecdf. Eventually, a resampling strategy, e.g. bootstrapping [START_REF] Efron | An introduction to the Bootstrap[END_REF], can be used to check the level of confidence in the obtained ecdf. Figure 3.3 plots the empirical cdf associated with samples of C M as well as its 95% confidence level.

Macroscopic thermal conductivity

The case of the macroscopic thermal conductivity tensor is more challenging, since it requires a steady-state thermal problem to be solved on a micro-structure to get one value of k M =

M f (k F , k S , n).
With isotropic component material properties k F and k S , and a symmetric geometry parametered with a single input n, the effective thermal conductivity tensor is isotropic, and one value for k M is observed. The polynomial chaos expansion (PCE) [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] is therefore herein selected as a tool for obtaining a surrogate model suited for uncertainty propagation. Indeed, since 3 random input parameters are involved, and since each evaluation of the model now relies on a FE computation on a micro-structure, it is less costly than a Sobol sampling to get the ecdf as done previously for the macroscopic thermal capacity.

The stochastic analysis of k M , taking consideration of the uncertainties of parameters (k S , k F , n) given in Table 3.2, is performed with PCE, and also with Monte Carlo (MC) and Sobol sampling [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF][START_REF] Sobol | Construction and comparison of high-dimensional Sobol generators[END_REF] for validation purpose, and evaluation of the surrogate modeling error.

In the case of large contrast in conductivities, the situation is simplified since

k F k S , so k M ≈ k SM = f (k S , n).
In the following, we therefore test the case without a strong conductivity contrast.

The precision of the different approaches in term of convergence for the mean and the standard deviation values requires different numbers of resolutions of the deterministic problem. For MC and Sobol methods, Figure 3.4 shows the (non-smooth) evolutions of mean and standard deviation errors when increasing the sample size (to compute an error, a reference value has been obtained by the over-killing approach: a much higher number of runs is used, 8000 for MC and 1000 for Sobol). A convergence is decided for MC with 3000 samples, and for Sobol with 400 samples. Asymptotic theoretical behaviors are known for these approaches, see e.g. [START_REF] Leobacher | Introduction to Quasi-Monte Carlo Integration and Applications[END_REF]: for MC, the error evolves as O(N -0.5 ), where N is the number of samples, while for Sobol, it evolves as O((log N ) p N -1 ), where p is the number of input parameters. The trends observed on Figure 3.4 show that Sobol sampling is more efficient.

Once an order p for the PCE is selected, r collocation points are considered, where ξ denotes inputs in [-1, 1] by normalization of physical inputs X. The coefficients of the PCE are then determined, to build the stochastic model M PCE A relative true error e T can be built on a set of points ξ j as

e T = mean j |M(ξ j ) -M PCE (ξ j )| mean j |M(ξ j )| (3.20)
If the set of points is the previous collocation set, this is a fit error. A reference error is nevertheless obtained with a larger (and possibly costly) set obtained for instance with MC or Sobol samplings. Once a reference solution has been obtained with the previous samplings, the precision of the PCE can be settled, depending of the order p. Increasing p indeed increases precision, but also the cost of building the surrogate model. The core of the cost leads in the FE nonintrusive computations performed at each collocation point. Therefore, the cost-precision curve of Figure 3.5a is depicted as reference error vs number of collocation points. The results concern the '3D' curves, where errors with respect to a reference obtained either with MC or Sobol are almost identical. '3D' is the size of the parametric space 3.

As for the capacity, using the PCE for generating a costless sample of conductivity allows to derive its ecdf, which is depicted in Figure 3.5b.

Since the cost is strongly influenced by the number of input parameters, the effectiveness of using dimensional analysis is also under concern. The dimensional analysis allows to reduce the number of parameters but usually complexifies the stochastic answers of the nondimensional involved parameters, as well as for inputs and outputs. The non-dimensional deterministic problem now reads: k M = g(a, n) (for the strong contrast case, it would be k SM = g(n)), involving ratios of random variables: k M = k M /k S and a = k F /k S . Uncertain inputs are a and n. Once an ecdf F a of a is obtained at low cost by sampling, the Probability Transform Method (PTM) [START_REF] Papoulis | Probability, random variables, and stochastic processes Fourth[END_REF] is used herein to transform the non-uniform distribution of a into an uniform one ξ: a = F -1 a (ξ), useful as input for standard PCE stochastic model using Legendre polynomial basis.

The same approach as for the dimensional case is performed and compared on Figure 3.5a, as '2D' curves. We can see the advantage of the dimensional analysis, the cost-precision curve being below the previous one.

Concerning the uncertainty propagation, as for the thermal capacity, we can see its influence at the macroscale as well, Table 3.2, that cannot be neglected and should be taken into account even for the macroscopic answer.

Macroscopic memory function

Since we are concerned also with the case with a large contrast in conductivities, the memory function has also to be dealt with. This case is the most challenging since it requires a full transient thermal problem to be solved with finite elements on the F phase of the microstructure.

We consider a test case with the scale ratio = 0.1, with k F = 2 k F . Now the stochastic model for k F is selected as the same as the previous one for k F in Table 3.2. Now the quantity of interest is a function of time:

β M (t) = f (t; C F , k F , n).
As this will be justified with the dimensional analysis, the time scale along which the vanishing function β M (t) is computed is set to t ∈ [0, 3τ F ] and discretized in 75 time steps.

τ F = C F l 2 /k F = C F L 2 /k F
is the characteristic time of the F phase; due to the strong contrast in conductivities, it exhibits the same order of magnitude as τ M .

Considering the uncertainties onto C F , k F and n, the macro memory function is described by a PCE such as

β M (t, ξ) ≈ M -1 i=0 β i (t)Ψ i (ξ) (3.21)
in which β i (t) are functions of time. Indeed, dimensional analysis can be used as well, and leads to: β M ( t) = f ( t; n). It boils down the number of input random variables to 1, that is n, and with a rescaled time line

t = t/τ F for β M ( t) = β M (t)/C F .
The same framework is used for the conductivity issue: the reference case is obtained for validation using a Sobol sequence of 400 samples, conforming to the convergence results of Figure 3.6.

Table 3.3 gives the PCE approximated mean and standard deviation of β M (t = 0) and β0 = βM ( t = 0) with respect to different order p and different models ('1D' and '3D'). Figure 3.7 gives the accuracy of β M (t) obtained with PCE meta model ('3D' for the dimensional model, and '1D' for the non-dimensional one) with respect to Sobol simulation, with increasing order p. For the dimensional case, the targeted precision is obtained with order p = 3 and 64 collocation points, for which the relative error is 0.017%. For the non-dimensional case, this is drastically improved to order p = 1, 2 collocation points, and a relative error 0.099%.

The mean, standard deviation and envelop of β M plotted with respect to the rescaled time line (t/τ F ) are post-treated and depicted in Figures 3.8 and 3.9.

We can note that in '3D' case, with p = 1, convergence seems attained on the mean and with p = 2 for the standard deviation, respectively though it is still not the case for the true error with respect to the Sobol sequence (still 17.94 % error) so this last one is a valuable (but costly) tool for assessing the accuracy, especially for reliability analysis for which not only the mean and the standard deviation are desired, but also the tails of the ecdf, see e.g. [START_REF] El Hami | Overview of Structural Reliability Analysis Methods -Part II: Sampling Methods. Incertitudes et fiabilité des systèmes multiphysiques 1. Optimisation et Fiabilité[END_REF]

Full macro-scale problem

The aim being to be able to perform a stochastic analysis, and therefore a many query problem, not only the macroscopic problem solving is required to be performed with a cost reduction, but also the evaluation of the macroscopic parameters. We therefore wish to rely on the PCE models of these macroscopic parameters previously built, together with a reduced model of the macroscopic resolution, in a form embedding PCE sub-models within the determination of a surrogate PCE of the macroscopic problem, leading to a kind of PCE 2 approach (inspired with the name of the FE 2 approach [START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fiber SiC/Ti composite materials[END_REF]), to cope with both multiscale and stochastic analysis.

Taking into account the memory effects leads to an increase in the complexity of the analysis, leading to a problem coupling 4 random inputs, when compared to a single random input for a standard macroscopic problem without memory effect. The QoI is field θ M (x, t; ζ) in which ζ is the random vector representing the uncertain parameters (C F , C S , n, k F , k S ) (described in detail in Section 3.5.4).

Macro thermal problem without memory effect

The classical macro transient thermal problem has been illustrated in Figure 2.6(c) with an assumed given (deterministic) thermal loading θ d (t) with an amplitude θ a and a characteristic time T c , as in (2.42). As previously, using dimensional analysis may be valuable. The dimensional macroscopic temperature is the solution of a problem whose inputs are:

θ M (x, t) = f (x, t; L, T c , C M , k M ) (3.22)
where an additional parameter T is the simulation time of transient thermal problem for t ∈ [0, T ]. We nevertheless enforce here T = 3τ M , where τ M = C M L 2 /k M is the macroscopic characteristic time. The simulation time is therefore no more an independent variable, and is discarded from the dimensional analysis. Let us introduce dimensionless variables as

ĈM = C M /C M = 1, kM = k M /k M = 1, L = L/L = 1, (3.23) t = t/τ M , x = x/L, θM = θ M /θ a , Tc = T c /τ M (3.24)
so that the non-dimensional problem reads:

θM (x, t) = f (x, t; Tc ) (3.25)
yielding to a single random input parameter Tc . One can then get back to θ M (x, t) using ratios of length, time and temperature: L, τ M and θ a respectively. Note also that the initial input random variables C M and k M are neither uniform, as can be seen in Figures 3.3 and 3.5b, nor independent, but since a single random input arises in the non-dimensional problem, all the stochastic aspects are concentrated into the parameter Tc . In this case, a stochastic analysis of the full problem can easily been performed as previously:

• a sampling of independent random microscopic inputs (C S , C F , k S , k F , n) allows to get a dependent sampling of (C M , k M ) with a low cost using the PCE model of Section 3.5.2;

• this last sampling allows to get a sampling of the non-dimensional random macroscopic input Tc = T c /τ M , hence its ecdf F Tc ;

• with this ecdf, a PTM allows to re-parameter this input with a uniform input single parameter ξ, as Tc = F -1 Tc (ξ);

• now the macroscopic black-box problem takes as input ξ and as output the field process θM (x, t);

• this black-box can be used to build a stochastic model with a non-intrusive PCE for the full stochastic macroscopic problem:

θM (x, t; ξ) ≈ M -1 i=0 θi (x, t)Ψ i (ξ) (3.26)
In order to compare this model with the other macroscopic model, the results are postponed to the next Section.

Macro thermal problem with memory effect

Carrying out a stochastic analysis of a macro problem with a memory effect (2.31) is more challenging, since (i) C M , k M and β M (t) are cross-dependent inputs, and (ii) k M and β M (t) are obtained from numerical simulations and therefore are not explicit functions of random inputs; they will therefore be replaced with their respective PCE non-dimensional models.

An intermediate step leads in building a surrogate model to the non-dimensional problem associated to (2.31), still with T = 3τ SM . With a dimensional analysis, this problem leads to

∂ θSM ∂ t = ∂ 2 θSM ∂ x2 + τ = t τ =-∞ βM ( t -τ ; n) ∂ 2 θSM ∂ t2 (τ )dτ (3.27) with θSM = θ SM /θ a , t = t/τ SM ∈ [0, 3], x = x/L ∈ [0, 1]
, and βM ( t; n) = 1 C βM ( t; n) is obtained with a time line rescaling as t = τ SM τ F t, and C = C M /C F . Additionally, the boundary condition, here a prescribed deterministic temperature, also exhibits a characteristic time

T c , whose non-dimensional equivalent is Tc = T c /τ SM = k(n) C × k S L 2 C F T c .
Therefore, a possible choice of random inputs are the non-dimensional quantities:

ζ 1 = k S Tc L 2 C F , ζ 2 = n and ζ 4 = C S /C F so that C = n + (1 -n)C S /C F = ζ 2 + (1 -ζ 2 )ζ 4 , and ζ 3 = k F /k S so that τ SM τ F = ζ 3 ζ 4 . ζ = (ζ 1 , ζ 2 , ζ 3 , ζ 4
) is needed to deal with the non-dimensional problem, but is not a set of random uniform independent variables. A surrogate model based of these, though using PCE approximation, is therefore not a stochastic model, but a simple surrogate polynomial model used only for cost reduction, similarly to existing reduced modelings, as in [START_REF] Casenave | A nonintrusive reduced order model for nonlinear transient thermal problems with nonparametrized variability[END_REF][START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF], as

θSM (x, t; ζ) ≈ M -1 i=0 θi (x, t)Ψ i (ζ) (3.28)
To do so, the following steps are performed:

• from uniform distribution of (C F , C S , n, k F , k S ) evaluate the ranges of ζ i , i = 1 . . . 4 for a rescaling in [-1, 1];

• select the associated Legendre polynomial basis;

• using surrogate collocation points, solve the macroscopic FE black-box problems: for each value of the set ζ,

evaluate C M analytically, k M with its PCE sub-model, β M ( t) with its PCE submodel

rescale β M ( t) to βM ( t) by ratios of their true values and time scales (with interpolation due to the time line potentially different discretizations):

βM = βM /C, t = t × ζ 3 ζ 4 .
with Tc and βM ( t) solve the full order macroscopic problem (3.27) to get the solution field θSM (x, t)

• fit macroscopic surrogate PCE model on the non-dimensional outputs at collocation points to get: θSM (x, t; ζ)

• rescale to obtain θ SM (x, t; ζ)

Even using a dimensional analysis, the thermal problem with memory effects still leads to a stochastic problem with 4 input random variables, but allows to embed as sub-models for the homogenized parameters the previously built costless PCE models of k M and β M (t). To assess the performance of the surrogate model, a comparison with the full order macroscopic problem is also set with a Sobol sampling sequence, for different orders p of the macroscopic PCE polynomials. Since embedding PCE sub-models, the comparison relies on the full solution strategy for obtaining the evolution of the temperature field θ SM (x, t; ζ) from the micro random parameters (C F , C S , n, k F , k S ).

Herein, the macroscopic QoI is the full evolution of the thermal field; if a micro temperature at specific point, or time is the QoI, a relocalization is needed to go back to the micro fields. The relocalization operators also could be provided with their respective PCE sub-models; nevertheless, if they are only required at a post-processing step (as for linear problems), this is not the computational bottleneck of the approach.

The accuracy of our proposed model in equation 2.31 to deal with the thermal memory problem has been studied and proven in section 2.4. Nevertheless, we check this deterministic model with inputs of mean variables in table 2, with applying the same boundary loading as shown in case of Figure 2 in Figure 3.12. For the macro thermal memory problem, it is computationally unaffordable to build reference case, due to the facts that thermal memory effect arises from thermal conductivity large contrast, and from the micro structure parameter n. Thus the accuracy of the surrogate model relies on the convergence of PCE mean and standard deviation. The mean and standard deviation of θ SM (x, t) are reported in Figure 3.14, together with the concentration on central location that are reported in Figure 3.15. We observe more uncertainty compared to the case without thermal memory effect. These results exemplify the feasibility to use reduced order models for parametric deterministic problems, such as those embedded into stochastic models where input parameters are the random quantities especially for uncertainty propagation, but could be used also for optimization problems. For the problem with memory function, the additional complexity of the model induces a large dispersion of the quantity of interest, emphasizing the interest in using the uncertainty propagation, for model validation and design issues.

Conclusion

This chapter carries out the PCE based stochastic analysis of the problem built upon deterministic numerical models in chapter one, being the homogenization, macro transient thermal problem, macro transient thermal problem with thermal memory effect, when considering the uncertainties of the material properties and geometry of underlying unit cell.

The cost of the stochastic analysis on C M with three input uncertainties (C F , C S , n) is negligible since the C M is analytically expressed.

The case of UQ analysis on k M exhibits more challenges, since it requires the steadystate thermal problem to be solved on a micro-structure to evaluate the effective thermal conductivity. Meanwhile, it involves three uncertain inputs k S , k F , n. We reduce the uncertainty dimension to two by using the dimensional analysis. It reduced computation cost compared to the 3-D problem ( original three uncertainties) while maintaining some accuracy order.

The thermal memory function is the time dependent macro property, requires solving the transient thermal problem on the micro scale, and three uncertainties C F , k F , n are involved. Before the PCE expansion, time discretization is necessary, the dimensional analysis is also applied to reduce the uncertainty dimension.

Uncertainty Quantification of two full macro-scale problems, with and without memory effects are performed. Due to five uncertainties (C S , C F , k S , k F ) existence, and their propagation to different scales, it is rather burden in computation requirement. We propose dimensionreduction and PCE 2 numerical strategies to release the cost. However, it is the random variable we consider in stochastic processing, the random field can be represented based on polynomial chaos decomposition [START_REF] Desceliers | Identification of chaos representations of elastic properties of random media using experimental vibration tests[END_REF]. p = PCE order ζ = uncertainty vector θ M (x, t) = macro temperature field as function of space x and time t θM (x, t) = non-dimensional macro temperature field θ M S (x, t) = macro temperature field considering thermal delay effect as function of space x and time θMS (x, t) = non-dimensional macro filed temperature field with considering thermal delay effect

Introduction

The optimization has been broadly applied in all domains. One of the most significant application is the topology optimization (TOP), to design the material layout to satisfy design targets. Since pioneer work [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] introduced the concept of topology optimization, it has taken time to develop in algorithms and computation hardware to be wildly adopted by industry. By now, the concept has been developed into several directions, element 'density' [START_REF] Sigmund | Topology optimization approaches[END_REF], nodal 'level set' [117] and 'basic components movement' [START_REF] Guo | Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[END_REF]; the fundamental difference relying on implicit/explicit representations of macro topology [START_REF] Guirguis | Evolutionary Black-Box Topology Optimization: Challenges and Promises[END_REF].

In this dissertation, we contribute to extending the classic TOP into systemic multiscale, multi-physics topology optimization. The multiscale TOP often produces geocentrically complex topology, even more, a design with complex micro architecture consisting of different phases, which are hardly fabricated by traditional tailoring. However, additive manufacturing (AM) is arising to overcome this recently. Thus, TOP is a key foundations of AM.

The extension of TOP to muli-scale cases is made possible by embedding homogenization model (refer to chapter 2), which is used to optimize the macroscopic layout as well as the micro configuration. One excellent use is for design of a structure consisting of periodically (or quasi-periodically) distributed unit cells, mesostructures, meta materials, architectured materials, lattice structures, or cellular structures. Designing optimized structure satisfying the objectives involving multiple physics inherits more challenges as: (1) multiphysics coupled finite element modeling, (2) multiple constraints required for multiphysical system modeling, (3) sensitivity analysis of objectivity to design parameters (when gradient based update scheme is used).

In this chapter, we extend the density-based topology optimization (refer to chapter 1) into multiscale, multiphysics, multiphase topology optimization (3M-Top), by embedding homogenization model (refer to chapter 2)). The outline of this chapter is organized as follows:

• Multi-scale. In section 4.2, we apply the non-gradient and gradient methods to solve multiscale topology optimization, taking into consideration fiber-reinforced composite material, which effective properties are analytically expressed by mixture-law and rotation matrix.

• Multi-scale, Multi-phases. In section 4.3, with considering two candidate composite materials, one investigates the effects of global mass constraint and candidate material property ratios of the elastic modulus to density on the phenomenon of two-phases competition in designed topology space. Note that, standard macro mechanical loading and a simple micro laminate-rank architecture are considered. For the purpose of numerical efficiency and achieving a good optimum, gradient optimization is adopted.

• Multi-scale, Multi-phases, Multi-physics. The multi-physics design problem (in section 4.4), involving thermo-mechanical coupled loadings, a global mass constraint, and considering locally two candidates, is built upon the multi-phases optimization mentioned above, it requires detailed derivation information for the gradient optimization solution. We will prove that in particular for the multi-physics loading problem, the multi-scale multiphase topology optimization will improve designed structural performance compared to monoscale topology optimization.

• 3M-Top complex architecture material. Finally, One embeds the complex micro architecture material into the 3M-Top framework, in which FEM based homogenization indeed provided an efficient numerical tool for obtaining effective thermal-mechanical properties of complex architecture material, and PCE is adopted to achieve a cheap semi-analytical surrogate model to the homogenized model.

Multiscale topology optimization

In this chapter, we employ a parametric microstructure, density-based macro topology representation, gradient-based and/or non-gradient-based optimization algorithms to simultaneously optimize the global macro topology and the local parametric cell structure.

To mathematically formulate this problem, two sets of design parameters

P = {ρ 1 , ρ 2 , ..ρ n , η 1 , η 2 , ..., η n }
are employed, including the pseudo-density parameters ρ to represent the material layout on the design space, and micro scale parameters η constructing the underlying multi-phases unit cell at different locations. The choice of η depends on the geometry characteristics of micro basic cell structure. For example, to fiber-reinforced composites material, η determines the orientation of fiber. To a porous material, η relates the micro volume fraction of the solid phase.

We apply both gradient and non-gradient optimization approaches to solve the optimization problem. The gradient-based approach requires firstly modeling the problem to calculate the objective function. Then, obtaining derivative/sensitivity information of objective function w.r.t. the design parameters is required to search the design variable along a direction. Finally, the optimizer updates the design parameters iteratively until it converges to a user-defined criterion.

In a non-gradient topology optimization approach, various heuristic optimizers are proposed to address the problem by taking an objective function as a black-box model but without sensitivity information. Therefore, the optimizer requires global search information of the objective space and the design variable space, and additional constraints are applied to drive the search towards a target design.

Another challenge relies on the cost of modeling the objective function. Therefore, instead of modeling the full-scale problem of heterogeneous material/structure, we adopted the homogenization surrogate model to assist multi-scale modeling and sensitivity analysis (when gradient-optimizer is adopted).

For the example of a thermo-mechanical loading design system, the parametric configuration at every macro point is associated with η. Given material properties of constituent and micro architecture, one can apply the homogenization method to obtain the corresponding effective elastic tensor D H , thermal stress tensor d H and thermal conductivity tensor C H . Afterwards, one can construct the surrogate homogenization model using PCE we discussed in previous section. Therefore, it avoids carrying out numerical homogenization at each optimization iteration to solve the effective thermal and mechanical properties.

Similar to standard density-based topology optimization, which use a most well-known SIMP material interpolation scheme [START_REF] Sigmund | On the design of compliant mechanisms using topology optimization[END_REF], we consider the same order of formulations to calculate the D M A i , d M A i and C M A i at ith macro element, which read:

D M A i = ρ p i D H (η i ) (4.1) d M A i = ρ p i d H (η i ) (4.2) C M A i = ρ p i C H (η i ) (4.3)
where ρ is strictly in range of zero and one, the p penalty factor of pseudo-density ρ, choosen to be 5, attempt to drive ρ (almost) binary. Starting with a purely mechanical problem, the general formulation of maximizing the stiffness of structural is expressed as:

min c(X) = F T U s.t. : KU = F V (X) ≤ V * (4.4)
where c is the structural compliance, the X is the set of the design variables, it consists of macro and micro geometry information. V (X) is the constraint problem, it will be specific in each case. We solve the problem with level set genetic algorithm and gradient method. We first take the example of a simple but most commonly used micro architectured material, the fiber reinforced composite material, that has the great advantages of high strength ratio and light mass, applied to aerospace vehicle engineering. We are interested in designing such composites, where the fiber is continuously distributed in a single pile. In such context, micro design variable η, is the orientation in range of [0, π]. The homogenized elastic module can be analytically expressed as:

D H (η i ) = T(η i )D 0 T(η i ) T (4.5) D 0 = 1 1 -ν xy ν yx   E x ν yx E x 0 ν xy E y E y 0 0 0 G xy   (4.6) T =   cos 2 η i sin 2 η i -2 cos η i sin η i sin 2 η i cos 2 η i 2 cos η i sin η i cos η i sin η i -cos η i sin η i cos 2 η i -sin 2 η i   (4.7)
where the T(η i ) is the orientation matrix and, D 0 is the original elastic matrix of fiber reinforced composite material without orientation. E, ν and G are Young modulus, Poisson ratio and shear modulus, with each of the subscripts x, y denoting different directions. The 2d geometry of such material is depicted in figure 4.1, within back inclusion fiber and surrounding white matrix material, assumed they are isotropic material of Young modulus are

E f = 1, E m = 0.5, identical Poisson ratio ν f = ν m = 0.3, G f = E f 2(1+ν f )
, volume fraction of fiber v f = 50%. For the first simple case, a mixture law can be used to predict the

E x = v f E f + (1 -v f )E m E y = (E f E m )/(vf E m + (1 -v f )E f ) G xy = (G f G m )/(v f G m + (1 -v f )G f ) (4.8)
It the following, we employ such a simple micro architectured material, exploring the multiscale topology optimizaiton using non-gradient based genetic algorithm and gradient based optimization algorithm.

Genetic Algorithm with multi-resolution transition scheme

Non-gradient based optimization (or derivative-free) requires no derivative information in the classical sense to find the optimum, when the derivative of the objective function is unavailable, unreliable or impractical to obtain.

Typically a non-gradient optimizer that uses global search techniques may have a high potential to converge towards a global optimum. As a wild heuristic optimization approach, an evolution Genetic Algorithm (GA) is common and wildly adopted [START_REF] Boichot | A genetic algorithm for topology optimization of area-to-point heat conduction problem[END_REF][START_REF] Prasad | Synthesis of bistable periodic structures using topology optimization and a genetic algorithm[END_REF]118]. The general principle of GA to assess the potential global optimum among the starting random of population of design variable X, is to keep the best populations (in terms of the objective function), and cross and mutate to get new generation children. The limitation of GA is evident because of the large population size, making the optimization a large query problem. Thus, we propose a two-resultant transitional scheme. It requires less decision variables in low resolution as inputs for GA optimization, afterwards obtained interpolated ones on the high resolution scheme as inputs for FEM to evaluate the cost function. The following one particularly demonstrates the framework of multi-resolution transition schemes, and the tuning parameters in genetic algorithm. Note that, we only consider the multi-scale optimization in the context of mechanical loading, in which the orientation of the unit cell rules the micro material properties.

Decision variable transition Firstly, a level-set interpolation scheme is used to describe the macro topology and micro orientation distributions. To reduce the number of decision variables, we proposed a two scale level set function translation approach, which defines the decision variables LSF c = {LSF 1 n1 i , LSF 2 n1 i } on coarse mesh (see in figure 4.2 (a)). LSF 1 ∈ [-3, 3] and LSF 2 ∈ [-π/2, π/2] are used to represent the macro topology and micro orientation distribution, respectively. Their superscript n1 denotes the mesh number i.e. the size of their discretization space. Afterward, processing interpolation, (see in figure 4.2 (b)), it allows one to determine the level set function LSF f = {LSF 1 n2 i , LSF 2 n2 i } on the finer n2 dimension mesh grid . Herein, a radial basis function is employed to obtain a smooth level set function on fine mesh grid, with given node position on coarse mesh {x} n1 i and its corresponding LSF value {f } n1 i

LSF f = N i=1 λ i ϕ i ( x -x i ) (4.9)
where the x is the fine mesh nodes location where we want to approximate the level set value.

x i is node locations on coarse resolution mesh grid. N is the mesh size at low resolution. λ i and ϕ i are the weights and the corresponding basis function.

Different common basis functions ϕ can be used, a basic multiquadric is used:

ϕ i ( x -x i ) = 1 + (ε x -x i ) 2 (4.10)
where, ε = 1 in our case. Then the weighs λ i can be estimated with linear least square method.

After the interpolation process, the macro topology has to be selected for the following objective function evaluation. The threshold T s is given to truncate the LSF f , consequently, an interaction is determined, which represent the boundary of solid and void material (see figure 4.3a). In uniform mesh domains, a single mesh is completely in or out of the intersection is solid or void material, respectively (see figure 4.3b). However, an element crossover in the intersection is regarded as intermediate density material; refining this element into 20 × 20 sub mesh structure (shown in figure 4.4) allows to select a sub-mesh within or crossing the level set as a solid material. The intermediate density of such element is measured by volume fraction of solid material in sub-mesh.

For such macro topology optimization design problem, normally, a structural volume fraction constraint V (LSF ) = V * should be satisfied by a line-search method to decide a threshold truncation T s ∈ [-3 , 3]. Therefore, at each iteration of GA optimization, evaluating the objective function consists of:

• interpolate decision variable {LSF n1 1 , LSF n1 2 } from coarse mesh into fine mesh {LSF n2 1 , LSF n2 2 } .
• with a prescribed initial T s , transform the level set function LSF n2 1 into density representation ρ and decide the volume fraction V .

• use a line search method to find a T s that guarantees the volume fraction constraint V * .

• Formulate the FEM to calculate the objective function c, with the macro density ρ n2 and micro orientation distribution LSF n2 2 . Note that, the density of the void material is assigned to a small positive value for avoiding the singularity problem.

In the context of GA optimization, the advantage of level-set interpolation would make it possible to construct the geometry with less design variables, meanwhile satisfying the quasi-periodic condition. Moreover, it avoids material non-connectivity issues when using 0 -1 binary represented structure [118].

GA parameters After the transition from decision variables to corresponding configuration via two resolution schemes, the FEM is performed to evaluate the cost/objective function, subsequently updated decision variables is performed with a black-box GA optimizer [START_REF] Chipperfield | A genetic algorithm toolbox for MATLAB[END_REF], where the solver only requires the objective value to proceed. • Select np initial diversity populations i.e. np configurations.

• Evaluate objective function of each configuration.

• Use genetic operators to update next generation populations from previous parent generations by: (1). elitism operation, keeping the best (in term of fitness) f e of total np parents. [START_REF] Afshar | Additive manufacturing and mechanical characterization of graded porosity scaffolds designed based on triply periodic minimal surface architectures[END_REF]. crossover operation, randomly crossing two parents among f c of the best parents, other than the elite ones. (3). mutation operation, adding f m chosen from a Gaussian distribution, to each entry of the parent vector.

• Stop if convergence or if the maximum generation number ng is reached.

Gradient based optimization algorithm

For two scale topology optimization subjected to mechanical loading F, the formulation is generally expressed as: min c(X) = F T U s.t. :

KU = F V (ρ) = ne i=1 ρ i ≤ V * η c1 ≤ η i ≤ η c2 0 ≤ ρ i ≤ 1 (4.11)
where, c is the compliance U and F are the nodal displacement and mechanical loading, X = {ρ, η} are the two scale design parameters. V is the volume constraint, η c1 and η c2 are lower value and upper value design parameters describing the micro architecture, ne is the total number of elements. K is the global stiffness matrix, assembled by elementary stiffness K i , that can be written as:

K = ne i K i = ne i Ω i B T D MA i BdΩ i (4.12)
where B is the strain displacement matrix, D MA i is defined in equation 4.1. The differentiation of c(X) with respect to the design variable X is expressed as:

∂c ∂X = ∂F T ∂X U + F T ∂U ∂X (4.13)
Considering the governing equation 4.11, we have:

∂K ∂X U + K ∂U ∂X = ∂F ∂X (4.14)
Then the derivative of the compliance c can be expressed as:

∂c ∂X = -U T ∂K ∂X U (4.15) = - ne i=1 U T i ∂K i ∂X i U i (4.16) = - ne i=1 U T i Ω i B T ∂D M A i ∂X i BdΩ i U i (4.17) finally , ∂c ∂X i = -U T i Ω i B T ∂D M A i ∂X i BdΩ i U i (4.18)
With the offline solution of homogenized surrogate model, D i H and ∂D H i ∂η i can be obtained with affordable computation cost.

Features control

The features control, as the filter approach discussed in introduction, is essential to gradient based approach to address the checker-board, mesh dependency problem. Meanwhile, it is employed to control features of the design structure. To obtain (almost) binary macro topology and micro spatially continuous architecture, a filter method is used to help the gradient-based optimizer to ensure global and local features control. In this study, a linear density filter is used as discussed in the introduction:

X i = j w ij X j (4.19)
where X is the design variable set, subscript i denotes the central element, subscript j denotes the j element within a certain distance to element i, w ij is the weighting factor. Conic weights w ij are defined as:

w ij = R-d(i,j) k∈N i (R-d(i,k)) j ∈ N i 0 j / ∈ N i (4.20)
where the R is the prescribed filter radius. However, in the context of multi-scale optimization, the design structure expects binary density variables ρ, but continuous micro architecture variable η therefore ready for additive manufacturing. Thus we propose a two-scales filter scheme, R ρ = 3 and R η = 3 in the beginning iterations, when to the end of convergence, turns R ρ = 1, driving a binary macro design structure.

Embedding the two-scale filters and modifying the sensitivity ∂ c/∂X i and ∂ M /∂X i lead to:

∂ c ∂ρ i = j∈N i ∂c ∂ρ j w ρ ij ρ j ρ i (4.21) ∂ c ∂η i = j∈N i ∂c ∂η j w η ij η j η i (4.22) ∂ M ∂ρ i = j∈N i ∂M ∂ρ j w ρ ij ρ j ρ i (4.23) ∂ M ∂η i = j∈N i ∂M ∂ρ j w η ij ρ j ρ i (4.24)
where the w ρ i,j and w η i,j are weight factors in Equation 4.20 with R = R ρ and R = R η , respectively.

Comparison of two approaches

Two approaches are used to solve the optimization of macro topology and micro fiber orientation distribution subjected to a half Messerschmitt-Bolkow-Blohm (MMB) Beam working condition (figure 4.5, where F = 1). The material properties of fiber f and polymer matrix m and the micro volume fraction of fiber material are

E f = 1.8, E m = 0.2, ν f = ν m = 0.3, v f = 45%
, respectively. The solution from gradient and non-gradient approaches will be presented and discussed.

In level set based GA approach, we define two level set function values LSF 1 ∈ [-3, 3] and LSF 2 ∈ [-π/2, π/2] on 11×6 on coarse mesh, use the proposed multi-resolution scheme and linear search method to construct FEM model at fine resolution mesh with a guaranteed global volume fraction constraint. GA parameter settings are population size, elite fraction, crossover fraction and mutation fraction, as: np = 200, f e = 5%, f c = 80% and f m = 20%. The maximum iteration number is ng = 600. In each GA iteration, the many query problem requires to solve the FEM problem with different inputs of populations, consequently, the cost of GA optimization is proportional to the np × ng. In mutliscale gradient optimization scheme, the design space is 40 × 20 rectangular mesh, the design variables are η ∈ [ -π 2 , π 2 ] and ρ ∈ [0, 1]. The filter approach is applied for the purposes of: [START_REF] Abueidda | Mechanical properties of 3D printed polymeric cellular materials with triply periodic minimal surface architectures[END_REF]. avoid mesh-dependent issue, [START_REF] Afshar | Additive manufacturing and mechanical characterization of graded porosity scaffolds designed based on triply periodic minimal surface architectures[END_REF]. obtain spatially continuously fiber orientation, and their filter distance are given as R η = 3, R ρ = 3. Derivatives of objective w.r.t. design variables have been implemented in equation 4.18, to iteratively guide towards an optimized design. We performed such design with different initial uniform orientations, which are η = 0, and η = π/2, and the corresponding final designs are presented in figure 4.8. The topology is related to the initial guess, however, we can obtain similar target design in terms of objective function, which are indeed local optimum. We perform designs on 100 × 40 fine mesh, have achieved a better design compared to low resolution designs in terms of smooth macro topology and design targets, because the intermediate density may under-evaluate the structural stiffness (or over-evaluate the compliance).

We observe a truss-like structure in figure 4.8, where fiber orientation becomes parallel to the direction of the dominant truss component. The smoothly varying fibers are present on the connections of each truss members. We obtain a similar macro topology from GA optimizer, and the figure 4.7 (b) is the best solution (but still not guaranteed as a global solution) to this problem with 40 × 20 mesh structures. GA may slightly win in achieving a better solution; however, it is more computational costly than the gradient approach. That is why in the future, we will use gradient-based optimization method. However, gradient-free is the only optimization solution to cases, with non-differentiable objective function (discrete bistable compliant structure) [START_REF] Ohsaki | Enumeration of optimal pin-jointed bistable compliant mechanisms with non-crossing members[END_REF] or non-connected design space (optimization of photonic or phononic band gap materials) [START_REF] Sigmund | Geometric properties of optimal photonic crystals[END_REF].

Mass constraint in multi-scale multi-phases topology optimization

In this section, we study a multi-scale topology optimization, where there are parametric cells spatially varying in design region. A unit cell equivalent to the original one in the fiber-reinforced composite optimization case is adopted. However, the micro design variable η is no longer the orientation but the microvolume fraction. Mass constraint plays an essential role in multi-phase topology optimization due to the material difference densities and lightweight design topic. Thus the interpolation scheme is not only for equivalent material properties (e.g., Young modulus, thermal conductivity, and thermal expansion) involved in the objective function but also for the material density involved in structural mass constraint. Finally, we apply a linear scheme to obtain the mass function; the mass constraint M (ρ, η) is:

M = ne i ρ i ρ η i ≤ M * (4.25)
where

ρ η i = η i ρ S + (1 -η i )ρ F
is the density of i-th element. ρ S and ρ F are densities of S and F phases. The optimization scheme has been illustrated in detail previously. The optimized outputs ρ η and the parametric unit cell are able to reconstruct the optimized multi-scale architecture.

It is the most basic case, where a micro parallel architecture is adopted here, thereby homogenized properties can be analytically expressed by mixture rules [START_REF] Lielens | Prediction of thermomechanical properties for compression moulded composites[END_REF]. This case works on (show figure 4.9a with applied Q = 0), consisting of 80 × 40 elements and three candidate materials (two non-zero f and m phases and a void phase) with a constant mass fraction constraint 50%. The Young's modulus of void material will be given by a small value E v = 1e -3 , to avoid singularity problem as wildly applied in SIMP like material interpolation scheme. Its corresponding density is set to be ρ

v = 1.e -3 , leading to E v /ρ v = 1.
Effect of different ratios E/ρ of non-zero phases on design structure is presented in table 4.1, where the first column is the candidate material mechanical properties. The second and third columns present the corresponding binary macro topology design and post-treatment of local density representation. Its color bar shows that deep blue, light blue, and red represent the void, lower density material, and higher density material, respectively. Continuous color regions between light blue and red indicate, at micro-scale, the two candidates material mixture with varying volume fraction. And the last column shows the results of objective function c, volume fraction v f , mass fraction of f and m phases mf f and mf m with respect to the corresponding cases.

In case 1, where E f > E m and

E f ρ f > Em ρm
, the approach designs a structure composed of only f material, which is indeed the benchmark result obtained by classical SIMP method. In case 2, compared to case 1, one decreases the ρ m , leading to increase the ratio Em ρm , however still maintaining the relation

E f ρ f > Em ρm .
Due to the mass constraint, the design is no longer single phase structure (only the stiffer f phase), however, consisting of both materials. It reduces the c but increases the volume fraction v f = 56.96%; the mass percentage of f and m materials are 44.3% and 55.64%, respectively. In case 3, compared to case 2, one increases ρ f = 1.5, leading to

E f ρ f = Em
ρm , the design increases c = 50.82, and reduces the consumption of f material to mf f = 25.14%. Last case considers

E f ρ f = 1 and Em ρm = 1.2, E f ρ f < Em ρm .
Compared to last case, it gives a similar macro topology design, but quite different in microscopic material consumption, c = 49.15, vf = 48.9% and mf f .

From cases 2 to 4, in the context of the maximum structural stiffness problem, the f phase is no longer the only choice for design (even if it is stiffer compared to another one) because of the mass constraint. An unexpected result was observed by comparing case 1 and case 4, indicating that the design structure mainly composed of material with stiffer and higher ratio E/ρ exhibits worse compliance, due to the mass constraint. One demonstrates a challenge of mechanical optimization of mass constraint problem, and even great a challenge for multi-physics loading case, due to the hybrid influence of E, k, E/ρ and Eα, to the best knowledge of the author, there is no previously case for such a result. Multi-scale topology optimization subjected to thermo-mechanical loading under the mass constraint is mathematically expressed as:

Materials

min c(X) = F T U s.t. : K t T = Q F t = A(T -T ref ) K m U = F = F m + F t M (ρ, η) ≤ M * η c1 ≤ η i ≤ η c2 0 ≤ ρ i ≤ 1 (4.26)
where, X is the design variables, K t , K m and A are the elastic stiffness matrix, thermal conductivity matrix, and 'thermal stress stiffness', respectively. F m is the design independent mechanical loading, F t temperature dependent thermal stress, it is design dependent. For simplification, T ref = 0. we adopt efficient gradient based optimizer, thus requiring to compute iteratively design sensitivities.

For efficiency of the optimization solution, it requires to explicitly express the derivatives of the objective function w.r.t to the design variables X. Such derivatives are carried out by adjoin sensitivity analysis; one express the objective function c in the Lagrangian form:

c * = (F m + F t )U -λ m (K m U -F m -F t ) -λ t (K t T -Q) -λ mt (A(T -T ref ) -F t ) (4.27)
Applying the chain rules, the derivative of the objective function with respect to X is expressed as:

∂c * ∂X = ∂c ∂X = (F m + F t -λ m F m ) ∂U ∂X -(λ mt A + λ t F t ) ∂T ∂X + (λ mt + λ t ) ∂F t ∂X -λ m ∂F m ∂X U -λ mt ∂A ∂X T -λ t ∂F t ∂X T (4.28)
This adjoin solution yields auxiliary equations:

U T + λ T m = λ T mt (4.29) λ T m K m = F m + F t (4.30) λ T mt A = λ T t K t (4.31)
where λ m T is a pseudo displacement solution, and λ mt T is a pseudo temperature solution, yielding from heat flux loading -K t λ mt T . Substituing such adjoint solutions into Equation 4.28, the total sensitivity is reformulated into:

∂c * ∂X = ∂c ∂X = -λ m ∂K m ∂X U + λ mt ∂A ∂X T -λ t ∂K t ∂X T (4.32)
Finally, the derivative of global c with respect to X i can be written as:

∂c ∂X i = n i=1 -λ (m,i) ∂K m i ∂X i U i + λ (mt,i) ∂A i ∂X i T i -λ (t,i) ∂K t i ∂X i T i = -λ (m,i) ∂( Ω i [B m ] T D M A i B m dΩ i ∂X i )U i + λ (mt,i) ∂( Ω i [B mt ] T d M A i B mt dΩ i ) ∂X i T i - λ (t,i) ∂( Ω i [B t ] T C M A i B t dΩ i ) ∂X i T i (4.33)
where the

D M A i , d M A i and C M A i
are the macro i-th element elastic tensor, thermal stress tensor and thermal conductivity, respectively. B m , B mt and B t are the strain-displacement , strain-temperature and flux-temperature, respectively.

For the homogenized properties estimation, it is low cost while using an analytical mixture law expression. However, with complex micro architectures, the mixture law is no longer precise, a numerical homogenization model is proposed. The details on homogenization and surrogate modeling technologies has been discussed previously. We could adopt such a PCE surrogate homogenized model to facilitate the sensitivity analysis in optimization. Method of Moving Asymptotes (MMA) optimizer [106] is used to interactively update the design parameters, until it is fulfilled with stopping criterions.

Numerical implementation

Numerical examples are implemented, firstly, a simple parallel two-phase microarchitecture is adopted, but the proposed approach is not limited to a certain parametric micro architecture.

When using candidate material shown in table 4.2, we study this case under the working condition sketched in figure 4.9a. The design domain is 80 × 40 mesh elements, the mass constraint M * = 80 × 40 × 0.5. We start with uniform initial guess ρ = 0.5, η = 0.5. The filter radius are R ρ = 2 and R η = 2, when it converges, R ρ = 1 is assigned to drive a almost binary macro density design.

Firstly, when applied Q = 0 and F = 1, it is exactly the mechanical problem that has been presented in the first row of the table 4.1, showing no phenomenon of material competition. Secondly, heat flux Q = 1 and Q = 2 are applied, respectively. The corresponding designs and temperature fields are represented in figures 4.11 and 4.12, respectively. Comparing the two designs, we note that a 'tail-like' component is sticking to the left side boundary where the heat source is located, consisting mainly of the m-phase. This component increases in size as Q increases. It is assumed to efficiently conduct the heat from the boundary through the rigid body, meanwhile expecting a lower thermal deformation.

However, notice that temperature near the heat flux region is exceptionally high, leading to the structure being non-working. Our proposed approach aims to minimize structural compliance, ignoring the risk of extreme temperature. To address this problem, the literature proposes temperature constraint topology optimization [134], or multi-objective formulation (minimize the thermal and mechanical compliance) [START_REF] Deng | Multi-objective concurrent topology optimization of thermoelastic structures composed of homogeneous porous material[END_REF],

We herein propose a simple but practical approach to overcome this issue by adding a non-design rigid domain material on the heat source boundary, as depicted in figure 4.9b, to efficiently transform the heat to the design domain Ω. The non-design domain consisted of f phase (material with higher thermal conductivity), which has same width and 1/80 length to the design domain. Figure 4.13 presented the designs obtained from this proposed nondesign approach, compared to design in figure 4.12, the objective function value is increased, meanwhile, the maximum temperature decreases dramatically, a smooth temperature field has been observed. Local optimum is wildly an issue in gradient-like topology optimization, especially when multiple candidate materials and mass constraint are introduced. Reference [START_REF] Gao | A mass constraint formulation for structural topology optimization with multiphase materials[END_REF] reports a uniform initial weighting in purely mechanical multiple-phase topology optimization problems with mass constraint. However, in thermo-mechanical coupling optimization problem, it is challenging to decide an unbiased initial guess, because temperature field, thermal stress are high related to the material distribution and macro topology. Therefore, for each numerical design examples, one re-run it considering different initial guesses to demonstrate that the proposed optimization approach can achieve a 'good' optimum. The sensitivity of the final solution to the initial guesses is one of the interest.

In the following, we implement the multiscale optimization subjected to different thermomechanical cases shown in figure 4.9, the different initial guesses to the target designs will be investigated. In addition, We show the corresponding monoscale optimization on f phase (the stiffer candidate material ) for each case, to demonstrate the advantage of the multiscale topology optimization.

Firstly, in working conditions shown in figure 4.9b, where the horizontal F m = 1 and Q = 2, test cases are implemented using candidate materials in Table 4.2 with fixed density distribution ρ = 0.5 and different initial values η. Note that, it is not necessary to satisfy the mass constraint at the beginning iteration. Table 4.3 reports their concerning designs and quantity of interest compliance c, volume fraction v f and f phase mass fraction corresponding to different uniform and non-uniform initial guesses.

In post-treatment of density representations, the deep blue, light green, and red colors indicate the void, f, and m material, respectively. Among 5 cases with a different initial guess, one produces five completely different macro topologies, but similar objective function values c, denoting that the design is highly related to the initial guess. However, it can always obtain a feasible design structure. With η 1 and η 2 , one achieves non-symmetric design structure, while obtained symmetric design with η 3 η 4 , and η 5 . In fact, one can explore the diversity of such designs with different initial guesses. However, there is still no guarantee that the global optimum has been found.

The last column shows the mono-scale topology design, when applying the same working and mass constraint as to the multiscale case, the objective value is c = 69.10, higher than any of the multiscale designs. Note that, if one executes the design with varying uniform initial density, however, they converge to the same target macro topology .

Secondly, we consider the L-bracket structure shown in figure 4.9c discretized with 80×80 elements. Applied mechanical and thermal loadings are F m = 1 Q = 2, respectively. A micro parallel architecture material consisting of two distinct materials shown in table 4.2 is considered here. The post-treated physical density is shown in figure 4.14, corresponding to different initial uniform ones η = 0.3, η = 0.5 and non-uniform one. There is no significant difference among the designed topologies, indicating that, in such working condition, the effect of initial guess has negligible effect.

Thirdly, for the workings shown in table 4.9d, the approach designs macro topology structure corresponding to different initial guesses are reported in table 4.4. There is no clear guidance on how to design the initial material distribution to avoid a local optimum in thermo-mechanical problem.

3MTop with TPMS micro structure

We consider one of the Triply Periodic Minimum surface (TPMS) family, Schwarz Primitive minimal surfaces [START_REF] Rajagopalan | Schwarz meets Schwann: design and fabrication of biomorphic and durataxic tissue engineering scaffolds[END_REF], to produce the parametric micro architecture. It is one of the minimal surfaces, where any points with a mean curvature of zero. They can be explicitly expressed with level set function φ(x, y, z) + iso = 0, in which φ(x, y, z) describes a surface calculated at the isovalue -iso that will produce a minimum surface-like topology. A Schwarz Primitive is considered herein, it reads: The mechanical and thermal properties of such structures with evolving iso ∈ [-2, 2] have been reported in literature [START_REF] Abueidda | Mechanical properties of 3D printed polymeric cellular materials with triply periodic minimal surface architectures[END_REF][START_REF] Cheng | Morphology, flow and heat transfer in triply periodic minimal surface based porous structures[END_REF][START_REF] Maskery | Insights into the mechanical properties of several triply periodic minimal surface lattice structures made by polymer additive manufacturing[END_REF]130], and the design spatially graded periodic minimum surface structure has been investigated [START_REF] Afshar | Additive manufacturing and mechanical characterization of graded porosity scaffolds designed based on triply periodic minimal surface architectures[END_REF]. Also, [START_REF] Challis | Design of three dimensional isotropic microstructures for maximized stiffness and conductivity[END_REF] reported that composites with the Schwartz primitive and diamond minimal surfaces as the phase interface have been shown to have maximal bulk modulus and conductivity. With adopted homogenization and PCE methods, the homogenized material properties with respect to cell geometry parameters are reported in appendix. Accordingly, the established PCE meta homogenization model will be applied to calculate the derivative information of targets to the design variables on both macro and micro scales (shown in figure 4.17).

f = cos(X1) + cos(Y 1) + cos(Z1) + iso = 0 (4.
We normalized material properties shown in table 4.6, found some interesting clues when solving design problem subjected to working condition in figure 4.15a:

• For optimization of the mechanical loading problem, the designed outputs are distinguished to applied different constraints, which are volume fraction, or mass fraction constraints.

• Since E M 2 > E M 1 , in volume fraction constraint context, one predicts that the design consists of one unique cell with iso = -1 (where the maximum volume fraction of M2 phase in microscopic ).

• However, due to

E M 2 /ρ M 2 ≈ E M 1 /ρ M 1
, we find that the two materials are in the same level of competition in context of mechanical loading and mass fraction constraints, which increases the potential to explore the diversity of the design space. (the result is shown in 4.15b)

• It is an engineering requirement to apply both volume and mass fraction constraints for fabrication convenience and lightweight designs (the result is shown in figure 4.15c)

• It is challenging to predict the level of completion of two materials in thermo-mechanical loading considering both constraints, for involving the heat loading may change the macro topology and the microphase distribution. (result is shown in figure 4.15d).

• Finally the corresponding 3D geometries (shown in figure 4

F th = 0, F m = 2, Mf * = 50% (c) F th = 0, F m = 2, Mf * = 50%, V f * = 60% (d) F th = 40, F m = 2, Mf * = 50%, V f * = 60% (e) F th = 0, F m = 2, Mf * = 50%, V f * = 60% (f) F th = 40, F m = 2, Mf * = 50%, V f * = 60%

Conclusion

In this chapter, firstly, we explored the performance of different optimization algorithms, the GA (non-gradient) and MMA (gradient), for solving the multiscale topology optimization. While using GA, we proposed a two-resolution transition approach and linear search method for reducing the dimension of design variables but maintaining the accuracy of the evaluated objective value. Furthermore, the linear search method guarantees the constraints at each iteration. In a multiscale topology optimization context, we are able to mathematically derivate the sensitivity information in mechanical (and thermo-mechanical coupling) problems as a result of the offline-built PCE surrogate model. The GA is a black box model for tackling almost all the optimization problems. the gradient method is efficient but requires plenty of offline works for having sensitivity information during the online optimization iteration. Secondly, we analyze the multiscale multiphase topology optimization in mechanical loading problem with considering the mass constraint, finding out that material competition can be due to the different ratios of stiffness to density.

Thirdly, we explore the multiscale multiphase multiphysics topology optimization, where we elaborately create two candidates material, one of them is stiffer with higher thermal expansion coefficient but with lower density, to analyze the effect of the additional thermal loading on the macro topology and micro material distribution.

Finally, we embedded the TPMS micro architecture consisting of real engineered material in the 3M-Top framework. Moreover, we applied both volume fraction and mass constraint to design a fabrication feasibility and lightweight structure.

However, uncertainty is always present in topology optimization for material properties, manufacturing process, so its propagation across the scale will be investigated in following charter.

Nomenclature

Introduction

The power of topology optimization (TOP) has been discussed previously. The previous chapter has been developed to explore the width of TOP, including the multiphysics, multiscale, multi-phases design system problem. In such context, deterministic conditions are assumed, disregarding various uncertainties that may arise from loading conditions, material properties, design variables. This leads the deterministic designs to be less robust (or more sensitive) to the present uncertainties. Thus, the formulation of the uncertainties is of tremendous significance under real-world engineering.

Both the uncertainty quantification and the optimization are many query problems; the computation affordability is the concerning aspect. In this chapter, we applied the PCE and gradient topology optimization to deal with the uncertainty propagation through the optimization system. The key ingredients rely on the (1) robust formulation and (2) sensitivity analysis.

Robust topology optimization: illustration case

In the deterministic TOP, we have a solution of design points X 1 to minimize the design target c; however, a real engineering problem always exhibits uncertainties ξ that will propagate to design system lead a varying c(X 1 , ξ). For demonstration purpose, figure 5.1 gives two solution points of c with respect to different designs X 1 and X 2 with same range of uncertainties ξ, It turns out that design X 1 is more sensitive to the uncertainty, while X 2 is the least sensitive one but with a less target performance compared to X 1 when absent uncertainties. It indicates that, in the context of uncertainties and topology optimization, one does not select the design point that minimizes the design target performance but instead selects the design points that reduce the fluctuation range of the target performance. However, it still needs to maintain a good target performance. Therefore, there is a trade-off between the best target performance and less sensitivity to uncertainty.

Robust topology optimization (RTO) is proposed, the objective function is often given as a weighted sum of the mean, and standard deviation of deterministic design target performance, a robust objective function c R is expressed as: min

X c R (X) = µ c (X) + λσ c (X) (5.1)
where the X is the design variables, λ is the weight factor λ, µ c and σ c are the mean and standard deviation of deterministic objective function with respect to uncertainties ξ. µ c and σ c are expressed as:

µ c (X) = E(c(X, ξ)) (5.2) σ c (X) = Var(c(X, ξ)) (5.3)
where the E and Var are operators for obtaining the mean and variance. With a given λ, applied constraint on design variables X, and the derivatives of µ c (X) and σ c (X) to the design variable X, one can adopt the gradient-based optimizaotn approach to iteratively drive to the design target. Note that, the probabilistic approach does not directly provide standard deviation but the variance. The derivative of ∂c R /∂X is expressed as:

∂c R ∂X = ∂µ c ∂X + λ 2 σ 2 c ∂σ 2 c ∂X (5.4)
There are alternative ways to formulate such robust optimization; one can minimize the mean [START_REF] Dunning | Robust topology optimization: minimization of expected and variance of compliance[END_REF], or one may imagine the way of minimizing the mean while constraining the variance. Many other statistical objectives do not include the statistics moments, but employ the confidence interval [131] or reliability metric to minimize the probability of exceeding some critical value [START_REF] Asadpoure | Robust topology optimization of structures with uncertainties in stiffness-Application to truss structures[END_REF].

Solving the RTO problem requires to iteratively express the derivations of ∂µc ∂X and ∂σ 2 c ∂X along the optimization process. The simplest solution, a brute Monte Carlo sampling [100], is employed to approximate the mean and variance, their derivations are expressed:

∂µ c ∂X = 1 nc nc i=1 ∂c i ∂X (5.5) ∂σ 2 c ∂X = 1 nc -1 nc i=1 2(c i -µ c )( ∂c i ∂X - ∂µ c ∂X ) (5.6)
where nc is the number of samples, c i is a deterministic objective function corresponding to sample i. Consequently, the solution cost of RTO is nc times the deterministic one, which is unacceptable when increasing the dimension of design variable X and uncertainty dimension ξ.

Herein, we employ stochastic response surface method , based on PCE and collocation points approach, to explicitly express the mean and variance with a moderate number of 'samples' [START_REF] Chu | Weakly intrusive time homogenization technique to deal with pseudo-cyclic coupled thermomechanical problems with uncertainties[END_REF]. The PCE method and standard simulation have been discussed in the chapter 3, herein we recall some key ingredients from the previous chapters that will help us realize the RTO.

Given a random ξ obeying to a certain probability density function, one selects the corresponding univariate polynomials. For example, with uniform ξ, Legendre Polynomials are selected. The objective c(ξ) is assumed to be the weighted sum of finite polynomial basis Ψ i (ξ):

c(ξ) ≈ M i=1 y i Ψ i (ξ) (5.7)
where M is the truncated threshold, and y i are the PCE coefficients. One may consider the roots of the polynomial of order p + 1 for selecting 1D experiment points, further design a higher dimension N experiments ξ k by tensor product, leading to r = (p + 1) N experiments.

Evaluate the unknown weight coefficients ŷ (or PCE coefficients) by linear square fit [START_REF] Schöbi | Polynomial-chaos-based Kriging[END_REF], reads: ŷ = Wĉ (5.8)

where ĉ being the vector of c(ξ k ), k = 1...r, W = (V T V) -1 V T is a M × r matrix with w j , j = 0, ...M -1 being each row vector. V is calculated from the evaluation of the polynomials basis onto each of the design experiments:

V = V ij def = Ψ j (ξ i ) , i = 1, . . . , r, j = 0, . . . , M -1 (5.9)
where V is M × r matrix, Due to orthogonality of the basis, one can explicitly give the mean µ c and variance σ 2 c , they are given as:

µ c = y 0 = w 0 ĉ (5.10) σ 2 c = M -1 j=1 ŷ2 j Ψ 2 j = M -1 j=1 Ψ 2 j (w j ĉ) 2 (5.11)
Consequently, the gradients of the µ c and σ 2 c can be given as: ∂µ c ∂X = w 0 ∂ĉ ∂X (5.12)

∂σ 2 c ∂X = 2σ c ∂σ c ∂X = M -1 j=1 2 Ψ 2 j w j ĉw j ∂ĉ ∂X (5.13)
The ∂σc ∂X is expressed:

∂σ c ∂X = 1 σ c M -1 j=1 Ψ 2 j w j ĉw j ∂ĉ ∂X (5.14)
Note that, the obtained expansion is only related to the PCE order p and polynomials basis. It allows approximating the mean and variance with ξ k collocation points. The precision of such an approximation relies on PCE order p. The derivative of ∂ĉ/∂X on collocation points is exactly the sensitivity in DTO problem, therefore, given a RTO problem considering N uncertainty variables, the PCE-RTO computation cost is roughly r = (p+1) N times DTO one.

Numerical validation on mechanical topology optimization

To summarize, the critical ingredients of RTO are: first benefits from using PCE and stochastic collocation method to approximate the mean and standard deviation with acceptable computation cost, second the gradient-based optimizer is used to update the design variables. Moreover, since PCE with collocation points is a non-intrusive method, it allows to express semi-analytically the stochastic moments, the derivation in robust optimization is no longer a trouble compared to the intrusive stochastic finite element method.

For validation, we implemented a macro case based on working conditions of figure 5.3a under a 40% volume fraction constraint. Material properties are E = 1, ν = 0.3, the 2D design domain is discretized by 80 × 40 elements. With applied deterministic mechanical loading F θ , the DTO case can be solved firstly to minimize the structure compliance.

In RTO context, F θ is mechanical loading with deterministic amplitude F = 1, however, the orientation of the loading is random, obeys a uniform probability low within range of

θ = [-π 2 -π 6 , -π 2 + π 6 ]
. With given different λ, one will obtain different macro topology designs.

First, convergence of PCE approximated mean and standard deviation when increasing p has been checked in figure 5.2 where the weight factor is λ = 0.1. Therefore, p = 2 is the best choice for a good moment approximation and a computation cost reduction. Other MMA optimizer parameters are the same as for a deterministic TOP case.

Two DTO cases are implemented where, θ 1 = -π/2 and θ 2 = -π/2 + π/6, finally producing 2 designs X d1 and X d2 depicted in figure 5.3b and figure 5.3c, respectively. One takes the θ 1 = -π/2 case as a 'worst' loading case, because the vertical loading to the supported structure will produce the maximum bending moments [START_REF] Hoang | Robust adaptive topology optimization of porous infills under loading uncertainties[END_REF]. With λ = 0.1, RTO design X r is produced (in figure 5.3d). In figure 5.4b, we perform the structural analysis of designs by DTO and RTO with respect to varying θ to calculate the corresponding compliance c. As one can observe, the RTO design and DTO design with θ = -π/2 are less sensitive compared to the DTO when θ = -π/2 + π/6. It indicates that, in this simple mechanical loading problem, the DTO to the worst loading condition and our proposed RTO can both have a more robust design. We then created the 8000 θ mc samples to measure the c(X, θ mc ) with respect to the 3 designs of figure 5.3, their histograms are plotted in figure 5.4a. We find the RTO design exhibits better robust performance compared to DTO with the worst loading condition, for the former one has the slight narrow data bar.

With λ varying, the RTO designs are reported in figure 5.5, and corresponding µ c and σ c are reported in table 5.1. Figure 5.6 gives the structural performance with those designs with respect to varying θ. It concludes that, increasing λ, improves the 'robust' performance of the structure, however reduces the structural performance in term of compliance. One may choose a suitable λ value depending on engineering requirement. We consider the different RTOs, where uncertainties being Young's modulus of m phase E m , loading orientations θ F and heat flux magnitudes Q, respectively. Table 5.3 reports the uncertain and the deterministic variables, where these uncertainties are independent and varying within given bounds, obeying uniform probability distribution. Note that, Q = [0, 8] is a severe changing environment, more than small uncertainties. We first consider individual uncertainties and then study problems combining them.

λ = 0.1 λ = 0.2 λ = 1 λ = 3 λ = 5
In a multiscale RTO problem, the same concept as in DTO is used. One decides design variables, macro pseudo-density x 1 , and micro-volume fraction of f phase x 2 to minimize a given target performance. The objective function is the weighted sum of the statistic mean and standard deviation of compliance in our case. For simplification, a parallel microarchitecture is adopted, therefore, homogenized material properties with varying x 2 are analytically expressed by mixture law. However, considering more complex micro-architecture in RTO may not bring too much difficulty; it causes no explosion of time cost compared to adopting a simple-architecture-like material, since the PCE based response surface method is adopted to map the micro geometry features (indeed the micro design variables in optimization system) of unit cell into the material space. It is an offline technology to obtain a semi-analytic surrogate homogenized model.

The objective/cost function c R is the weighted sum of µ and σ with a weighted factor λ to the σ. An uncertain quantification technology, PCE with collocation points approach is used to approximate the µ and σ. SIMP-like material interpolation is adopted and detailed in section 4.2.

Adaptive algorithms

One solves the RTO with a gradient optimizer, so with the possible drawback of getting a local optimum design. Generally, a multi-start with different initial points will somehow avoid this problem. Herein, we apply three different strategies to compare the performance of the convergence, the capability of achieving a 'good' optimum, etc. The first algorithm is the common one, where applying uniform initial guess on design mesh. The second one is more elaborate; the coarse mesh design problem generates the design as an initial guess to the refined mesh design problem. Note that this generated design on coarse mesh, as an input does not necessarily need to converge to a binary density topology. The initial guess for the coarse mesh problem is still uniform but can be userdefined. Unlike direct uniform starting, it has excellent advance in convergence speed to a feasible design. This non-uniform starting approach is indeed a progressive refine-mesh method (PRM). Specific in our case, one performs RTO on 20 × 10 mesh; constructs a rough design, where intermediate pseudo-density (where design between 0 and 1) is not an issue. Typically, it corresponds to a solid component structure on a finer mesh domain. Afterward, one continually performed RTO in 40 × 20, and 80 × 40 mesh domain along with inputs generated by the previous solutions. For the optimization period of 80 mesh, we adopted p = 5, and r = 1, since the high penalty factor, and low density filter radius will accelerate the convergence speed and eliminate intermediate density, respectively. Applying higher p at the beginning stage of optimization will lead to a non-convex problem, making it more difficult to find a good optimum. This PRM strategy is similar to the multiresolution strategy proposed in [START_REF] Boucard | A multiscale strategy for structural optimization[END_REF]. Both of them reuse results coming from a 'coarse' design process.

The last strategy is called a Continuation Penalty (CP) scheme, where the applied penalty factor p = 1.5 at the beginning continuously increases to 3 during the optimization process. The structure compliance objective function c shows convexity while the stiffness is linear to the design variable (in our case, the macro density and micro volume fraction) [105]. It indicates that the choice of smaller p leads to an undesirable situation, where intermediate design material existing, however, achieves a good local optimum. So, our proposed CP approach is presumably a good initial guess at primary stage with p = 1.5, and from there to increase p gradually to achieve the solution without intermediate density. Finally, for an obtained binary macro topology, we assign filter radius r = 1 at the last step.

In another aspect, both the CP and PRM are intended to provide a good initial guess (or diverse initial guesses) to the final optimization stage; they are named adaptive strategies used to progressively change parameters (mesh size to PRM and p to CP). introduces more computation burden. We applied the 3 different strategies to explore the diversity of generated designs. Table 5.2 lists the topology optimization input parameters with respect to different strategies.

Uncertain heat flux magnitude

The DTO have been solved in preceeding section, we are now tackling the corresponding RTO case taking into consideration the uncertainty of heat flux. Typically, an uncertain range is twenty percent to its mean. In our RTO case, the heat flux is in the range of [0,[START_REF] Auriault | Upscaling Heterogeneous Media by Asymptotic Expansions[END_REF] for building an enormously varied heat flux case. On one hand, we are particular interested in considering zero heat flux in the range of uncertainty. Since it is a pure mechanical loading case, where we will achieve a monophase design structure consisting of a stiffer material between two candidate materials. On the other hand, one can physically consider it as a time-evolving thermal work condition (imagine longer machine works, the more heat it generates). We implement RTO to reduce the sensitivity of compliance to the heat flux loading changing. Firstly, with λ = 1, using three different strategies with their concerning topology optimization inputs listed in the table 5.2, RTO is performed with different initial guesses as shown in the first line of table 5.4. It is noted that, we apply the same x 1 (the macro density initial guess, x 1 = 0.3), but different x 2 (micro design parameter in each element, two uniform, one non-uniform) for different cases. Moreover, we play a trick here, the non-uniform case, where x 2 is a random field on the mesh domain, obeying uniform probability low. Indeed, the effect of macro topology initial guess is interesting and well studied in topology optimization [START_REF] Bendsoe | Topology optimization: theory, methods, and applications[END_REF]. However, in our case, considering effects of different initial guesses at both macro and micro scales is more complicated. The other lines in table 5.4 show the corresponding designs, means, standard deviations. By observing the designs 5.4, we have several concluded points:

• All the proposed strategies perform well for the converging ability with the given λ = 1, even with the random field initial guess x 2 .

• The RTO-UNI is highly biased by the initial guess, while different starting points guide towards different macro topology designs.

• RTO-PRM and RTO-CP are unbiased by the initial guess point. Even with a random initial x 2 guess, the design will converge to a topology similar to others achieved by applying uniform initial guesses.

• These presented solutions show the feasibility to improve the structural robustness, nevertheless, it may yield local optimum designs, which is a common issue of gradientlike optimization approach Tables 5.5-5.7 report the effect of weight factor λ to the RTO-(PRM,UNI,CP), where different cases start with same initial uniform guess x 1 = 0.3, x 2 = 0.3.

• Firstly, Table 5.7 reports RTO-UNI designs concerning λ = 1, 4, 8, 12. One admits that the approach fails to converge with λ = 2 resulting in the missing displayable designs between λ = 1 and λ = 4. The local optimum issues have been observed among these discrete points. For instance, in the context of λ = 8, the design corresponding to λ = 4, however, has a better performance in terms of objective function value (c R = 10.852). Nevertheless, increasing λ somewhat increases µ, while minimizing σ.

• Secondly, Table 5.5 reports RTO-PRM designs concerning λ = 0.01, 0.1, 0.5, 1. One observes a similar macro topology pattern with increasing λ but decreasing mass fraction of f phase (the stiffer phase of two candidate materials), and σ rises while µ drops. One observes a similar macro topology pattern with increasing λ but decreasing mass fraction of f phase (the stiffer phase of two candidate materials), and σ rises while µ drops. When λ is over 1, converging to an exhibitable design is difficult, possibly due to m f nearly reaching to 100%. Thus, this pattern macro topology is no longer suitable, but RTO-PRM fails in jumping to a new macro topology pattern.

• Lastly, Table 5.6 reports RTO-CP designs concerning λ = 0.1, 1, 2, 4. This method is unbiased by initial guess, it has good behaviors in convergence ability with different λ compared to others methods.

Uncertain mechanical loading orientation

To demonstrate the performance of RTO with the mechanical loading uncertainty, the boundary conditions in figure 4.9d, is employed, where the load fluctuated with an angle θ = µ θ ±σ θ . µ θ is used for a deterministic case, and the variance σ θ corresponds to a RTO case.

It is a remarkable result in our RTO for a pure mechanical uncertainty problem that a worst-case DTO design corresponds to a solution that the structure sustained the most significant bending moment, which improves the structural 'robustness' performance subjected to the uncertain environment. However, it is not easy to decide the worst case in thermomechanical coupled loading problem, because the design-dependent thermal deformation may counterbalance the mechanical one. This highlights the contributions of our proposed method in Robust Topology Optimization with mechanical uncertainty problem.

Currently, RTO is employed to minimize the µ c and σ c considering µ θ = 0 and σ θ = π/12. RTO-UNI and adaptive RTO-PRM and RTO-CP are employed. The effects of λ are reported in tables 5. [START_REF] Auriault | Upscaling Heterogeneous Media by Asymptotic Expansions[END_REF] • Table 5.9 reports the different RTO-CP designs w.r.t. varying λ = 0.1, 1, 2, 8. This method shows a significant advantage in exploring diverse macro topologies. To the first three designs, as for the RTO-UNI and RTO-PRM, RTO-CP performs somewhat shape optimization to a baseline topology to slightly adjust the σ and µ. However, in the last case, one increases the λ a lot, it optimizes the macro topology to dramatically reduce the σ.

Uncertain material properties

Taking into account the uncertain material properties is one of the interest in the RTO problem. Apart from the thermo-mechanical properties reported in table 5.3, elastic modulus is one of the most interesting uncertainties, because it is related to both the mechanical and thermal stress. Besides, we do not have any knowledge of the nature of variation, one considers first a uniform random material variable. Moreover, a random field is expanded in several random variables using PCE, therefore, our proposed RTO solves the system of equations of higher uncertain dimension. Herein, we considered that the elastic modulus E m obeys to the uniform probability, varying between [0.6, 1]. Tables 5.11 -5.13 report mean µ and standard deviation σ to the

Conclusion on individual uncertainties

We perform a robust multiscale topology optimization considering the independent individual uncertainties of heat flux magnitude, mechanical loading orientation, and material elastic modulus. The objective function defined as the weighted sum of mean µ and standard deviation σ is our choice to formulate the robust optimization. Nevertheless, alternatives, e.g., minimize µ while constraining σ and probability described µ and σ, are of interest.

Local optimum is a wildly known common issue. Thus we applied RTO-UNI, RTO-PRM, RTO-CP to explore the possible 'good' local optimum. In addition, the effect of λ on the RTO has been investigated. By combined all the solutions of different optimization search algorithms (UNI, PRM, CP) and different weight factors λ, we give a pseudo-pareto-front of RTO design with heat flux uncertainty ( Figure 5.10), and mechanical loading orientation uncertainty (Figure 5.11). Performing RTO with applied different λ is a brute force way to construct a completed pareto-front line. We could apply high and low values of λ first, then refine the range progressively. However, it is still a heavy burden task, for (1) our RTO is a costly computation problem (compared to classical macro topology optimization problem), considering two-scale design variables and uncertainties; (2) we do not know how the standard deviation is sensitive to the varying λ, it requires an adequate test to choose a suitable range of λ to design the Pareto front; (3) with a given λ, a convergence solution sometimes is not guarantee. One may refer to a Pareto Optimal Tracing (POT) approach that has been applied to efficiently trace the Pareto front from a single solution already had [START_REF] Nathan | A Multi-objective Approach for Robust Structural Topology Optimization[END_REF]. However, this method may help but not guarantee to reduce the cost in our problem, for a solution obtained with a given λ as an input for a new RTO problem (different λ) may have a convergence issue.

Combined uncertainties problem

We have already designed the RTO corresponding to a single uncertainty. However, in this section, in the context of the multi-physics loading problem, we will exemplify the robust topology optimization taking into account the combination of loading uncertainties. As discussed before, the weighted sum of the mean and standard deviation of thermo-mechanical compliance is the chosen objective function to minimize. Combined uncertainties in RTO are none of challenging to implement. Nevertheless, it is a dramatic computation cost since the cost of RTO is in proportion to the DTO to the given number of experiments used to evaluate these statistical moments precisely. We are not interested in validating the precision of the PCE approximated µ and σ for every optimization iteration, due to the linearity of the deterministic objective function. Indeed, one can decide an adequate PCE order p = 7 to enforce the accuracy of the estimated mean and variance.

RTO-PRM has the faster convergence speed. However, even with different λ, it often converges to a similar pattern macro topology. Note herein, when one evaluates the structural compliance of each design subjected to the samples of uncertainties, the penalty factor p = 1 is applied because we have the binary (almost) topology.

RTO-UNI is the most direct way. For comparison purposes in this case, we only produce two RTO-UNI-(Q, θ) designs shown in figures 5.13 (c) (d) with λ = 3 and λ = 5, respectively, and then compared to RTO-UNI design considering sole uncertainty(either Q or θ, as shown in 5.13 (a),(b)). Figure 5.13 (e)-(h) are the corresponding structures' response surface subjected to different loadings (samples of heat flux and orientation uncertainties). Figure 5.13 (e),(f) indicates that the RTO design shows stronger sensitivity to absent uncertainty (not considerated in a RTO problem) compared to the presented one (considerated in a RTO problem), illustrating an achievement of our proposed RTO design in improving the robustness to considerated uncertainties. Figure 5.13(g) shows the sensitivity of RTO-UNI design subjected to both uncertainties with applied λ = 3. The trend of c(X, Q, θ) is similar to the 5.13(e), however, the former one has more minor variance (compared the range of colorbar of two plots). When considering two uncertainties and applying a higher λ = 5, the c(X, Q, θ) is a almost flat surface if one applies same level color bar as in figure 5.13(e).

We implement RTO-CP-(Q, θ) cases, and table 5.14 depicts the several designs with a given λ, and the corresponding µ and σ, and figure 5.12 plots the response of compliance to the uncertainties concerning different designs.

Firstly, with applied smaller λ and considering either individual or combined uncertainties (compared design cases of RTO-UNI with different uncertainties recourse), one obtains a similar pattern of typologies, indicating the robustness of the RTO-CP strategy. Secondly, with the figure 5.12, one can conclude how sensitive the design with respect to different uncertainties.

Thirdly, we observe an intermediate density region (the lighter blue close to the green region) by a RTO design with given λ = 1.88. Note herein, using SIMP-like material interpolation, for example, the homogenized elastic property is expressed as x p 1 E H (x 2 ), where x 1 ∈ [0, 1]. With p > 1, one relaxes the material space to achieve a differentiable objective function to the x 1 . Normally, for topology optimization problem p = 3. For a binary solution (x 1 = {0, 1}), different p have no influence on evaluating the c(x 1 , x 2 , Q, θ) (such as cases in figure 5.12 (a),(b),(d),(e)). However, for a structure with intermediate density region, applying

x p 1 E H (x 2 ) (x 1 ∈ [0, 1], p > 1) to evaluate c(x 1 , x 2 , Q, θ) lacks of physical interpretation.
Therefore, we perform c(x 1 , x 2 , Q, θ) with p = 1 (the intermediate density material interpretation), p = 3 (the material relaxation interpretation) and p = 5 (an experiment test), shown in figure 5.12 (c), (f) and (g), respectively. They show difference in sensitivity to uncertainties, however, have the almost same mean value (depicted with the color of surface).

We have several inspiration points by the RTO hybrid uncertainties case:

• The presented uncertainties, in fact, could be interpreted as range of use to the structure, therefore, we can manipulate the structural target performance (as same concept for the sensitivity in context of robust optimization) by topology optimization.

• An intermediate density material can be interpreted as porous microstructure materials in a particular case; indeed, x p 1 , x 1 ∈ [0, 1] with varying p (p ≥ 1) can be physically interpreted as different porosity density. Accordingly, microscopic material design may provide a significant potentiality for optimization and a more feasible space for manipulating material/structure performance.

We implement the RTO considering the complex architecture which corresponds to DTO design considered in section 4.6. The results are shown in figure 5.14 It is not managed to discuss details, for example effects of different λ, single or multiple uncertainties, or the performance of the different RTO adaptive strategies, however, only for demonstration purpose. Note that, volume fraction and mass fraction constraints are both guaranteed.

It denotes that, with our proposed multiscale topology optimization (and robust optimization), we can consider arbitrary parametric cell (as shown in figure 1.6), and then reconstruct the full-scale structure from optimized macro topology and micro geometry features.

Conclusions

This chapter extends the gradient 3M-Top framework in the previous chapter by considering the uncertainties to deal with a robust topology optimization problem. In such context, the new design target is the weighted sum of the mean and standard deviation of the structure compliance subjected to the uncertainties. We apply the non-intrusive PCE stochastic analysis (details are illustrated in chapter 3), to approximate the mean and standard deviation, consequently, the sensitivities of the robust objective function to the design variables (micro and macro) are easily calculated. Gradient optimizer MMA is used to solve the RTO problem.

For illustration purposes, we start with a macro topology optimization to consider the loading orientation uncertainties. The accuracy of mean and standard deviation with different PCE orders p are validated. In addition, the effect of λ on the designed topology and PCE obtained mean, and standard deviation have been analyzed. To the end, we find that: (1) increase λ typically increases the mean and decreases the standard deviation (i.e., improve the robustness of design); (2) the deterministic optimization subjected to the worst loading case (usually, the loading corresponding to the maximum bending moment) improves the structure's robustness. However, RTO can supply designs with varying mean and standard deviation by changing λ, which is not possible in a deterministic optimization (the worst loading case). Furthermore, one may not be able to decide the worst case when suffering from complex thermo-mechanical coupled loading.

We apply the proposed RTO solution for solving the independent uncertainty of heat flux magnitude, mechanical orientation, and elastic modulus in 3M-Top cases, on a hybrid uncertainty problem. One direct uniform initial guess (UNI), and two adaptive strategies, the progressive refined mesh (PRM) and the continuously increase of the penalty factor (CP), are proposed. We have validated that, RTO-PRM and RTO-CP approaches are less biased to the initial guesses.

Our proposed PCE based RTO with adaptive strategies (UNI, PRM, and the CP) and the PCE-surrogate homogenization model are appealing in solving the multi-scale topology optimization problem with uncertainties in the following aspects: (1) non-intrusive implementation. (2) computation affordable: for using the gradient optimization, stochastic collocation method, scale separation, PCE surrogate homogenization model to the FEM based homogenization model. (3) Compatibility. It is ready for RTO design with different micro parametric architecture materials and different uncertainty sources (random variable, field, or described by an empirical probability density function).

Conclusion

This thesis explores and contributes to several areas and topics, including multiscale multiphysics multiphases topology optimization, homogenization, stochastic analysis, gradient and non-gradient optimization algorithm, robust optimization. Specifically, it mainly contributes to developing a systemic robust topology optimization that takes into consideration of (1) multiphysics loading, (2) multi-phase micro architectured material, (3) multiple independent uncertainties, aims to reduce the design sensitivities to uncertainty. The target of robust optimization is quite clear. However, the most challenging task is the computation requirements, since it is a many query problem to evaluate the performance of structure that embeds different macro topologies and micro material distribution and on uncertainties. Therefore, some problems should be pre-solved: (1) heterogeneous structure performance analysis, (2) uncertainties propagation, (3) multiscale topology optimization algorithm developments. They are concerned with several essential technologies: (1) homogenization (chapter 2), (2) uncertainty quantification (chapter 3), (3) topology optimization and multiscale topology optimization (chapter 4). Chapter 1: Chapter 1 provides a primer of density-based topology optimization, which consists of material interpolation, penalty, optimization algorithm, to satisfy some design targets by the layout of the material. It is well-developed in overcoming common issues in TOP, like mesh dependence and checkerboard pattern, leading to a manufacturing feasibility design structure. Literature related to the trend of development on multiscale optimization and robust optimization has been surveyed. To extend the width and depth of TOP, we are interested in developing a systemic robust multiscale multiphysics multiphase (3M) topology optimization in an affordable computation manner.

Chapter 2: In chapter 2, scale separation of heterogeneous structures is achieved by homogenization. Instead of performing analysis on a full-scale structure, one applies the multiscale-homogenization method based on the periodicity (or quasi-periodicity) assumption to perform the analysis on a macro scale.

We take the transient thermal problem as an example, developing the FEM-based homogenization method on an underlying unit cell to obtain the homogenized macro properties. As a result, two classical thermal conductivity, thermal capacity, and one non-classical thermal memory effect in the context of large contrast of thermal conductivity of two phases are obtained. Details of numerical developments i.e. applying boundary conditions, solutions of each property component, to obtain these homogenized material properties have been demonstrated. In addition, it is ready for other physics problems of heterogeneous material, mechanical, and coupled thermo-mechanical (in Appendix), and can be extended to others as acoustic, magnetism, etc.

For validation purposes, we perform three transient thermal cases on structure scale, i.e., the full-scale, the macro scale with and without thermal memory effects, to validate our obtained homogenized three thermal properties. As a result, the proposed thermal memory function and homogenized thermal conductivity and capacity show good accuracy and efficiency in predicting the thermal delay effect compared to the full-scale case.

Chapter 3: We start from a deterministic heterogeneous transient thermal problem, for the physical quantities are both spatially and temporally related. As a result, the present uncertainty analysis to this system is computationally costly. Thus it requires developing an easily implemented and affordable computational method for stochastic analysis.

For this purpose, the chapter 3 proposed a non-intrusive PCE-based stochastic method, for analysis of the uncertainties in unit cell built upon deterministic numerical models in chapter 2, being the homogenization, macro transient thermal problem, macro transient thermal problem with thermal memory effect, and considering the uncertainties of material properties and geometry of underlying unit cell. In stochastic analysis processing, one applied dimensional analysis to reduce the number of uncertainty parameters. Consequently, three stochastic models are formulated, i.e., the Monte Carlo simulation, PCE-based stochastic model, PCE-based reduction model (with applied dimensional analysis), and compared in terms of stochastic moments convergence and computational cost. Besides the full-scale problem, we applied the dimension-reduction and PCE 2 numerical strategies to release the cost.

This chapter has a significant foundation on the robust multiscale optimization for: (1) a non-intrusive stochastic method, and (2) a semi-analytically surrogate homogenization model development (indeed the PCE 2 method).

Chapter 4: The chapter 4 is the exploration of the width of density-based topology optimization in the developments of:

• multiscale and multiphase topology optimization;

• non-gradient and gradient-based optimization algorithm;

• applied in multiphysics coupled loading problem. Firstly, we applied parametric cells, as a formula of material, into the density-based topology to develop a multiscale topology representation. For each integration point of the FEM model, the pseudo-density and micro parameters are the design variables for (1) evaluating the full scale structure performance and (2) updating by an optimization algorithm to satisfy design targets.

Secondly, we develop the non-gradient and gradient based optimization i.e., genetic algorithm (GA) and Method of Moving Asymptotes (MMA), respectively, to solve the problem. In GA, we develop a two-resolution and linear search scheme to reduce the design variable and guarantee the constraints. In MMA, the sensitivity analysis is derived with respect to the macro and micro design variables. It requires (1) the material interpolation scheme for a relaxation density ∈ [0, 1] (2) offline homogenized material property by an analytical formulation (in our case, either the mixture law and rotation matrix described simple ranked parallel micro structure material with rotation angular, or the PCE surrogate homogenization model for arbitrary parametrilized complex architecture material). We prove that both optimization algorithms are feasible to solve the multiscale topology optimization problem through a numerical case of a continuously oriented fiber-reinforced composite material. GA is an ideal non-gradient processing applied to a non-differentiable objective function. Gradient-based MMA is efficient compared to GA. However, it may be stuck into local optimum problems, especially in a complex design system with multiple design variables.

In the following, we apply MMA to solve our different cases,

• two material competition under mass constraint and mechanical loading (applied simple ranked parallel microarchitecture);

• two material competition under thermo-mechanical coupled problem, (applied simple ranked parallel micro architecture), we achieve a better design compared to the monoscale optimization;

• the multiscale optimization case considers a Triply Periodic Minimum Surface (TPMS) micro architectured material.

Chapter 5: With the development of the previous method, the robust multiscale optimization is naturally solved. The stochastic moments (mean and standard deviation) are approximated by PCE stochastic method, via a weighted sum of the collocation points on the uncertainties. We apply MMA to solve the robust optimization. Nevertheless, the local optima and the numerical convergence are still the main issues. Two adaptive strategies, i.e., the progressively refined mesh (PRM) and continuously increasing penalty value (CP), have been developed. Compared to the simple uniform initial guess (UNI), both the PRM and CP show a better convergence ability. PRM is the most efficient. However, it converges to a similar pattern of macro topology, reducing the possibility of exploring the global optimum. The CP is the most desired method, for its convergence ability, diversity of explored topology, and unbias to initial guess. We have applied this method for dealing with individual uncertainties of heat flux magnitude, mechanical orientation and material properties variable, and dealing with combined uncertainties in the robust optimization process.

Outlooks and future direction

We have achieved our academic targets, however, some aspects remain for further improvements to be applied to solve real engineering problems:

• A 3D multiscale optimization development. This is merely an extension to be tested on problems of a larger scale.

• In a thermo-mechanical topology optimization problem, one chooses structural compliance as design target performance. Mechanical and thermal compliance are interesting aspects of achieving a stiff and efficient thermal transfer structure component in real engineering problems. A weighted sum of the mechanical and thermal compliance is one of interest. It is effortless to extent such design targets in our proposed multiscale topology optimization and robust optimization.

• The uncertain variable is considered in robust optimization. However, in a real problem, the uncertainty features with the random field, for example, the uncertainty surface loading, spatially uncertainty material properties. It may increase the uncertain dimension dramatically, and need to consider further reduction development.

• Nonlinear material models have to be also tacked. In this case, the fully non-intrusive GA is straightforward but not computationally efficient. Therefore, an advanced gradient-like approach needs to be tested, as well as advanced homogenization and multiscale technologies.

We have implemented the robust multiscale optimization and robust optimization embedded with TPMS micro architectured material and modeled the full-scale structure for demonstration purposes. In such cases, the filter method smooths the spatially distributed micro parameters at the macro scale. However, the microgeometry mismatching occurs because of the functionally graded isotropic architectured material. Therefore, it has significant side effects on structure performance, i.e., stress and thermal concentration problems. Typically, damage happens on the geometry miss-match region. We can handle this problem by either post-treatment modeling the full-scale structure or a fundamental approach by parametrically describing the functional graded anisotropic architectured material. In addition, material space is extended by the latter method (isotropic and anisotropic material space).

Non-gradient is a potential solution to topology optimization, it is a black box mode optimization approach but time exhausting. However, it is a hot point in design material with applied machine learning. Where databases are lacking, unbiased data generation can be achieved (for example, populations generated by GA ).

As we mentioned in robust optimization, the presented uncertainties could be interpreted as varying evolving stimuli to structure. A topology optimization approach can program the structural to a target performance subjected to given stimuli and apply it in designing unique microarchitectured material. where σ is the Cauchy stress tensor defined by the Duhamel-Neumann law for thermoelasticity, being the strain tensor. D is the elastic module tensor, and the d is the thermal stress increment per unit temperature. T and T ref are the structural temperature and reference temepearture while zero strain, for simplification, assumed T ref = 0. It is supposed to be satisfied to the boundary condition, it is written as:

σ • n = t bc in Γ t u = u bc in Γ u
where, u bc and t bc are displacement and traction tensor on boundary Γ u and Γ t , respectively. Displacement u can be expressed by an asymptotic expansion series in:

u(x, y, t) = u 0 (x, y, t) + u 1 (x, y, t) + 2 u 2 (x, y, t) + . . .

Submit it into govern equation and by applying the terms with expansion series -2 , -1 and 0 equal to 0, it allow to solve the elastic and thermal characteristic displacement w.r.t. different orders. First, given terms to expansion power of -2 to zero, one obtained:

∇ y (D∇ y u) = 0 (7.4) lead to u 0 is the macroscopic displacement, and the macro scopic strain can be expressed as

M = ∇ x u 0 .
Then, terms to -1 gives:

∇ y (D∇ x u 0 ) + ∇ x (D∇ y u 0 ) + ∇ y D(∇ y u 1 ) -∇ y dT 0 = 0 (7.5)

Considering ∇ y u 0 = 0, we can have relationship of ∇ y u 1 to ∇ x u 0 and T 0 .

∇ y u 1 = -Ly ∇ x u 0 + Ly T 0 (7.6) where the Ly , Ly are liner operator, for build relationship between macro displacement and micro displacement and temperature filed, respectively. Terms to expansion power of 0 gives:

PCE stochastic analysis on homogenized thermomechanical properties

The close approximation to the structure level surface consists of trigonometric functions F at points of (x, y, z), and they are formulated as:

F (x, y, z) = cos(x) + cos(y) + cos(z) -t (7.12)

Parameter t is the variable that determines the interphase region, separating the unit cell into two regions of different phase material. Moreover, if the parameter t is sufficiently positive and negative, the inclusion will be globally surrounded by the matrix, which turns an unconnected inclusion into a full heterogeneous scale. As a result, one defines t ∈ [-0.9, 0.9] to guarantee the connectivity of inclusion trough the macro scale, figure 7.

Random field uncertainties

The main issue is an increase in cost when the number of parameters increases. Indeed using a dedicated parametrization of the microarchitecture, and eventually of the materials of the phases, the number of variables may be large, especially with the versatile manufacturing process as for 3D printing. In any case, parametrization is a potentially useful approach to ensure manufacturability (and as a second hand, be compatible with quasi-periodicity assumption). Therefore, taking into account a random field material parameter is to be taken into consideration concerning the computational cost. With random parameter fields, p(M, ξ), we can approximate these as a Karhunen-Loeve expansion with a single product: p(M, ξ) ≈ p 0 (M ) * (1 + f (ξ)). Indeed we mainly need an estimate of the standard deviation of the cost function c, and this approximation will probably lead to an overestimate of it. For instance, f (ξ) = α p ξ when ξ ∈ [-1, 1] and α p is a deterministic amplitude (so we got a mean and a uniform variation, PCE order 1). Nevertheless, we may consider different random variables ξ for the different parameters. Note that p 0 (M ) are deterministic fields, so the context is similar to the one when only the given load is random.

Obtaining a macro characteristic q M from p makes use of an hybrid between PCE and ROM: we therefore now denote p(M, ξ) = g(Ξ), where Ξ is a parameter allowing to describe both the nominal value range, and the variations around a nominal value. We now have q M = i Q I ψ(Ξ) as previously, but with this different interpretation of variable Ξ. A nominal value of q M is therefore obtained from a the corresponding nominal value of the micro parameters, so the nominal value Ξ 0 of Ξ.

The previous objective was c(p), now with p(Ξ). Since Ξ has both a deterministic part and a stochastic one, we write Ξ 0 its nominal (deterministic) value and ∆Ξ its interval of variation so that Ξ = Ξ 0 ± ∆Ξ.

The new objective is µ + λσ where µ(Ξ 0 ) is the mean and σ(Ξ 0 ) the standard deviation of c(Ξ); λ is a weight.

We first make the following approximation: µ ≈ c(Ξ 0 ) (mean and nominal value). For gradient-like optimization, we need the gradient of µ, which is easy to obtain: ∂c/∂p or ∂c ∂Ξ (Ξ 0 ) = ∂c ∂Ξ 0 (Ξ 0 ), so it is the same as for the deterministic case. We now also need the gradient of σ. With a first-order sensitivity analysis, we got ∆c = | ∂c ∂Ξ |∆Ξ. To get rid of the absolute value, let's consider (∆c) 2 = ∂c Since ∆Ξ is directly related to the standard deviation of Ξ, ∆c also allows to get σ. The last step will be to express the gradient of σ. Les nouveaux procédés de fabrication, comme la fabrication additive, permettent de réaliser des géométries de pièces mécaniques complexes, mais aussi de microachitectures multiphasiques. Cette flexibilité peut augmenter les possibilités d'optimisation topologique et paramétriques de ces pièces, en particulier quand elles sont soumises à des sollicitations multiphysiques. Néanmoins, des incertitudes de fabrication peuvent aussi intervenir. L'objectif de cette thèse est de combiner les méthodes d'homogénéisation, de propagation d'incertitudes et d'optimisation de la topologie pour formuler une approche robuste d'optimisation multiéchelle, pour déterminer à la fois la forme macroscopique de la structure, et l'évolution microscopique de l'architecture du matériau, de façon robuste, pour satisfaire les performances cibles. Les applications présentées visent donc à pouvoir concevoir la macrotopologie et la distribution des micromatériaux pour une rigidité maximale de la structure, et à réduire la sensibilité de la structure aux incertitudes présentes (par exemple, le chargement thermomécanique et les propriétés des matériaux) en vue d'augmenter la robustesse de l'approche. Ce cadre général proposé a la capacité de résoudre les problèmes d'optimisation en considérant des cellules architecturées paramétrées multiples, des problèmes de chargement complexe, des incertitudes hybrides, etc., avec un coût de calcul abordable.

MOTS-CLÉS

  (a) Functionally graded bamboo material, left bamboo culm, right functionally-graded structure[START_REF] Shah | Mapping thermal conductivity across bamboo cell walls with scanning thermal microscopy[END_REF] Solution by a multiscale topology optimization in section 4.6. Q and F are the heat flux and mechanical force loading, respectively

Figure 1 . 1 :

 11 Figure 1.1: Multiscale functionally graded structure (a) by nature observation (b) by multiscale topology optimization design

Figure 1 . 2 :

 12 Figure 1.2: density-based topology optimization

Figure 1 . 3 :Figure 1 . 4 :

 1314 Figure 1.3: The flowchart of density-based TOP method

Figure 1 . 5 :

 15 Figure 1.5: Macro and micro scale topology optimization technique maps, red box is our approaches

  (a) Triply periodic minimal surface unit cells[START_REF] Maconachie | SLM lattice structures: Properties, performance, applications and challenges[END_REF] (b) Strut-based lattice structures[START_REF] Maconachie | SLM lattice structures: Properties, performance, applications and challenges[END_REF] (c) topology optimized microstructure[START_REF] Andreassen | Design of manufacturable 3D extremal elastic microstructure[END_REF]

Figure 1 . 6 :

 16 Figure 1.6: Three types of 3D unit cells

Figure 2 . 1 :

 21 Figure 2.1: Illustration of scales separation

Figure 2 . 3 :

 23 Figure 2.3: Mesh of unit cell

Figure 2 . 4 :

 24 Figure 2.4: Thermal flux continnuty check

Figure 2 .

 2 Figure 2.5 shows the homogenized time dependent β M with respect to different time resolutions (nt = 15, 45, 75, respectively).With applied nt = 45, one achieves a good convergent β M compared to coarse time discretization. Since the accuracy of this proposed numerical solution to solve β M has been validated compared to a analytically solution[START_REF] Dureisseix | Numerization of a memory effect for an homogenized composite material with a large contrast in the phase thermal conductivities[END_REF]. In next section, we will validate it from a different perspective: comparing the solution of equation 2.31 to the solution obtained by full-scale modeling.

Figure 2 . 5 :

 25 Figure 2.5: β M /τ F w.r.t. different time steps nt

Figure 2 . 6 :

 26 Figure 2.6: (a) details of unit cell (F phase in red, S phase/matrix in blue), (b) transient thermal boundary-value problem in a 3D composite and finite element mesh, (c) same problem using effective media

Figure 2 . 7 :

 27 Figure 2.7: Comparison of the macroscopic part of the temperature for 3 loading cases: (a) T c /τ M = 20, (b) T c /τ M = 2, (c) T c /τ M = 0.2, (d) T c /τ M = 2 while the classical macroscopic thermal characteristic time is τ M = C M L 2 /k M = 200.

Figure 3 . 1 :

 31 Figure 3.1: Univariate Legendre polynomials up to 4th order

Figure 3 . 2 :

 32 Figure 3.2: 2D nodes based on the same 1D grids. Left: Sparse grids. The total number of points is 145. Right: Tensor product grids. The total number of nodes is 1,089 [126]

Figure 3 . 3 :

 33 Figure 3.3: Histogram and ecdf of the homogenized capacity C M

Figure 3 . 4 :

 34 Figure 3.4: MC and Sobol convergence properties for the homogenized conductivity k M

  Figure 3.5: (a) PCE (2D for non-dimensional problem, 3D for dimensional one) meta-models cost-accuracy performance, (b) use of the PCE to derive the histogram and ecdf, for the macroscopic conductivity k M

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Convergence of Sobol sampling for β M function

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Mean and standard deviation (std) of β M

  (b). The Figure3.10 reports the norm errors of |θ SM -θ S |/θ a and |θ M -θ S |/θ a with respect to time evolution and intersections of adjacent cells. The accuracy of surrogate model with PCE order p to the macroscopic thermal model is reported in Figure3.13, and p = 6 is settle down to construct surrogate model. The maximum minimum and mean, together with standard deviation of θ M (x, t) are post-treated and depicted is Figure3.11. In particular, the envelope and statistic characteristics of θ(x = L/2, t) are illustrated

Figure 3 . 10 :Figure 3 . 11 :Figure 3 .Figure 3 . 13 :Figure 3 . 14 :Figure 3 . 15 :

 3103113313314315 Figure 3.10: Absolute errors of temperature fields along x direction

NomenclatureM

  = numerical model X = uncertain inputs Y = uncertain outputs µ y = mean value σ y = standard deviation f Y (y) = probability density function {Ψ j } ∞ j=0 = orthogonal basis function {y j } ∞ j=0 = coefficients to corresponding orthogonal basis function E = mean operator N = uncertainty dimension M = truncation threshold p = order of polynomial µ P C Y = mean approximated with PCE-stochastic collocation method σ P C Y = standard deviation approximated with PCE-stochastic collocation method C = thermal capacity k = thermal conductivity C M = homogenized thermal capacity k M = homogenized thermal conductivity kM = non-dimensional homogenized thermal conductivity n = micro volume fraction M P CE = PCE surrogate model τ = characteristic time β M (t) = homogenized thermal memory function βM ( t) = non-dimensional homogenized thermal memory function β i (t) = PCE coefficient to homogenized thermal memory function

Figure 4 . 1 :

 41 Figure 4.1: Fiber reinforced composites micro 2D structure

Figure 4 . 2 :

 42 Figure 4.2: Multi-resolution transition scheme, (a) Decision variable at nodes of coarse uniform grid LSF , (b) Interpolate on Finite element mesh grid LSF f

Figure 4 . 3 :

 43 Figure 4.3: (a) level set function on finite element mesh grids, and a threshold boundary with respect to a given T s . (b) Level set function to density representation

Figure 4 . 4 :

 44 Figure 4.4: Measure the intermediate density

Figure 4 .

 4 6 presents the cost function of the best population in each iteration. The optimized macro topology and micro orientation design structure with multiple restarts are illustrated in figure 4.7 as well.

Figure 4 . 5 :

 45 Figure 4.5: MMB working condition

Figure 4 . 6 :

 46 Figure 4.6: GA: iteration vs fitness value

  (a) GA c = 242 (b) GA c = 102.22

Figure 4 . 7 :

 47 Figure 4.7: GA generates different solutions by restarting

Figure 4 . 8 :

 48 Figure 4.8: Optimization solutions subjected to working condition shown in figure 4.5. Effects of different initial guesses and mesh sizes on the design targets are demonstrated.

  Topology Local densityResultsE f = 1.5, ρ f = 1; E m = 0.8,ρ m = 2 c = 50.44 vf = 50% m f = 100% E f = 1.5, ρ f = 1; E m = 0.8,ρ m = 0.5 c = 44.07 vf = 56.96% m f = 44.36% , E f = 1.5, ρ f = 1.5; E m = 0.8,ρ m = 0.8 c = 50.82 vf = 55.17% m f = 74.86% E f = 1.5, ρ f = 1.5; E m = 1.2,ρ m = 1 c = 49.15 vf = 48.9% m f = 93.42% Table 4.1: Mechanical loading case with different materials properties 4.4 Topology optimization with thermo-mechanical coupled loading 4.4.1 Compliance minimization

  phase 1.5 10.e -4 2 0.3 1 1.5

Figure 4 . 9 :

 49 Figure 4.9: Different boundary conditions, (b) (c) (d) are cases with the grey non-design region

Figure 4 . 10 :

 410 Figure 4.10: Optimization convergence proof

Figure 4 .

 4 Figure 4.11: Q = 1, c = 70.13, vf = 48.89%, mf f = 95.56%, mf m = 4.44%

Figure 4 . 12 :

 412 Figure 4.12: Q = 2, c = 78.32, vf = 47.37%, mf f = 89.46%, mf m = 10.54%, may not be a feasible design if f bearing high temperature

Figure 4 . 13 :

 413 Figure 4.13: Q = 2, c = 108.11, vf = 46.62%, mf f = 90.22%, mf m = 9.78%

Figure 4 . 14 :

 414 Figure 4.14: F = 1, Q = 18. (a), (b) and (c) are local optimum designs w.r.t. different initial η to working condition 4.9c. (d) is the mono-scale topology optimization.

3 : 4 :

 34 Multiscale designs. Red and green represent the m and f phase respectively. c, v f , m f are the compliance, volume fraction, mass fraction of m material, respectively. η 1η 3 are uniform initial guesses. η 4 and η 5 are non-uniform ones. Multiscale topology optimization subjected to working condition 4.9d. Different designs w.r.t. different initial guesses. The last row is the corresponding mono-scale topology optimization where X1 = 2απx, Y 1 = 2βπy and Z1 = 2γπz, α, β, γ are constant values, related to the unit cell size in x, y, z coordinate, and we prescribed α = β = γ = 1 to consider an isotropic material.

Figure 4 . 15 :

 415 Figure 4.15: Solutions to different parameter settings subjected to boundary condition (a). Mf * , V f * , c are mass fraction, volume fraction, and objective value, respectively. (b) Pure mechanical loading case with mass fraction constraint, c = 38.89. (c) Pure mechanical loading case with mass and volume fraction constraint, c = 41.33. (d) Thermal mechanical loading case with mass and volume fraction constraint c = 61.37. (e) convergence of V f, Mf constraints corresponding to (c) solution. (f) convergence of V f, Mf constraints corresponding to (d) solution.

Figure 4 .

 4 Figure 4.16: (a), (b) and (c) are the corresponding full scale 3D geometries (only M2 phase) to the solutions shown in figures 4.15b, 4.15c, 4.15d. (d) shows how to manage M1 and M2 phase in micro structure.

Figure 4 . 17 :

 417 Figure 4.17: Thermo-mechanical properties evolving with respect to iso

P

  = design variable sets ρ = density design variable η = micro geometry features design variables D H = homogenized elastic modulus D M A i = homogenized elastic modulus of ith element d H = homogenized thermal expansion coefficient d M A i = homogenized thermal expansion coefficient of ith element C H = homogenized thermal conductivity C M A i = homogenized thermal conductivity of ith element F = mechanical loading U = displacement field K = global stiffness matrix V (X) = volume fraction T = rotation matrix LSF c = level sets on coarse mesh LSF f = level sets on fine mesh N = mesh size at low resolution ψ = basis function λ = weight factor to basis function T s = truncate threshold np = number of initial population for GA f e = elitism fraction to initial population for GA f c = crossover fraction to initial population for GA f m = mutation fraction to initial population for GA ng = number maximum generation for GA η c1 , η c2 = lower and upper bound for micro design variable B = shape function X = filtered design variable c = filtered compliance M = filtered mass fraction Q = heat flux vector T ref = reference temperature A = global thermal stress stiffness matrix F t = equivalent thermal force loading F m = mechanical force loading K m = global stiffness matrix c * = augmented Lagrangian compliance {λ m , λ t , λ mt } = adjoint vector R ρ = filter radius of density ρ R η = filter radius of micro geometry η B m = strain-displacement shape function B mt = strain-temperature shape function B t = flux-temperature shape function 94

Figure 5 . 1 :

 51 Figure 5.1: Deterministic optimization and robust optimization, c , X and ξ are the objective function, design variable and the uncertainty respectively.

Figure 5 . 2 :

 52 Figure 5.2: µ c and σ c versus RTO iterations w.r.t. different PCE orders p, λ = 0.1

Figure 5 . 3 :

 53 Figure 5.3: Compared designs of DTO and RTO. (a) The loading condition. F θ 1 and F θ 2 refer to two different DTO cases. The red lines show the range of uncertain loading orientation. (b) The worst loading DTO case, θ 1 = -π/2. (c)DTO case with θ 2 = -π/2 + π/6. (d) RTO case with loading orientation uncertainty ranged in [θ 2 = -π/2 -π/6, -π/2 + π/6], weight factor λ = 0.1, PCE order p = 2.

Figure 5 . 4 : 5 Figure 5 . 5 :Figure 5 . 6 :

 5455556 Figure 5.4: Visualize the 'robustness' of designs. (a) x-axis indicates the range of compliance, and y-axis represents the corresponding frequency distribution

Figure 5 . 7 :

 57 Figure 5.7: Working conditions to the design problem

Figure 5

 5 Figure 5.8: RTO-PRM-θ, difference

Figure 5 . 10 :

 510 Figure 5.10: Pseudo-Pareto front plot for RTO heat flux uncertainty

Figure 5 .

 5 Figure 5.11: Pseudo-Pareto front plot for RTO mechanic orientation uncertainty

  (a) λ = 0.5 (b) λ = 1 (c) λ = 1.88 (d) λ = 2 (e) λ = 3 (f) λ = 1.88 (g) λ = 1.88

Figure 5 . 12 :Figure 5 . 13 :

 512513 Figure 5.12: Sensitivity of RTO-CP design w.r.t different λ subjected to combined uncertainties (Q, θ). Applied unified colorbar. (a)-(e) with applied p = 1. (f) with applied p = 3. (g) with applied p = 5

  λ

Figure 5 .

 5 Figure 5.14: (a) and (b) are RTO designs with different λ. (c) and (d) are the corresponding full scale geometry.

NomenclatureX

  = design variable vector ξ = uncertain variable vector c R = robust objective function µ c = mean of compliance w.r.t. uncertainty ξ σ c = standard deviation of compliance w.r.t uncertainty ξ E = mean operator Var = variance operator ∂µ c /∂X = derivative of mean to the design variable ∂σ 2 c /∂X = derivative of standard deviation to the design variable ξ k = design experiments based on collocation points p = penalty factor, PCE order 1 N = numbers of uncertain dimension Q = heat flux amplitude uncertainty θ = mechanical loading orientation uncertainty θ mc = random samples on mechanical loading orientation uncertainty r = filter radius nelx × nely = mesh size x 1 = macro design variable x 2 = micro design variable

  -mechanical homogenizationThermal mechanical equilibrium equation is written as:div σ + f = 0 (7.1) σ = Dd(T -T ref ) (7.2)

∇

  x (D∇ x u 0 ) + ∇ y (D∇ x u 1 ) + ∇ x (D∇ y u 1 ) -∇ x (dT 0 ) -∇ y (dT 1 ) = -f (7.7) submit equation 7.6 into the equation above, and apply average operator x ,∇ x D 1 -Ly x M + D Lyd x T 0 = -f (7.8)note that, u 1 x = 0 and T 1 x = 0σ M = D M M -d M T (7.9) This leading to homogenized elastic tensor D M and thermal stress tensor d M : D M = D 1 -Ly x (7.10) d M = d -D Ly x (7.11)

2 × 10 - 5 [ 1 , 3 ] × 10 Figure 7 . 1 :

 2105131071 Figure 7.1: Structure of inclusion phase with t ∈ [-0.9, 0.9]

∂Ξ 2 (

 2 ∆Ξ) 2 .

  FOLIO ADMINISTRATIF THESE DE L'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : CHU DATE de SOUTENANCE : Prénoms : Chenchen TITRE : Multiscale and multiphysics robust design of a complex microstructure with uncertainties, and driven by target performances NATURE : Doctorat Numéro d'ordre : 2022LYSEI022 École doctorale : MEGA (Mécanique, Énergétique, Génie civil et Acoustique) Spécialité : Génie Mécanique RESUME :

  : matériau architecturé paramétré, homogénéisation, propagation d'incertitudes, optimisation topologique multi-échelle, optimisation robuste Laboratoire (s) de recherche : Laboratoire de Mécanique des Contacts et des Structures (LaMCoS), UMR5259 Directeur de thèse : David Dureisseix Président de jury : Composition du jury : Pierre-Alain Boucard, Christophe Desceliers, Boris Desmorat, David Dureisseix, Béatrice Faverjon, Nawfal Blal

Table 1 .

 1

	remarks		Multiobjective evolutionary	algorithms; NSGA II non-	gradietn optimization		25 tensor product quat-	drature points to 10000MC	points		-		-		a intrusive method; precon-	ditioned Conjugate Gradi-	ent (PCG) solver to large	linear system	Multiscale TOP without	the homogenization and	local volume constraints	Variance is not considered		Design dependent uncer-	tainties; intrusive method	Multi-Resolution Finite	Element Models within an	inexpensive low-resolution	model	Multiscale topology optim-	ization; Hybrid uncertain-	ties
	UQ		Tauguchi L16 ar-	ray design experi-	ments		MC; Gauss-type	quature sampling			Gauss-type quad-	rature	MC and quadrat-	ure points	PCE and	stochastic Galer-	kin		MC			Analytical formu-	lation	-		MC			hybrid univariate	dimension reduc-	tion
	Uncertainty rep-	resentation	-				K-L				Uniform variable		EOLE		PCE				principle of linear	superposition		-		Taylor series		K-L			Uniform interval
	Ref TOP Recour		[49] DTO load orient-	ation; mag-	nitude;material	properties	[25] LST geometric shapes	and load-	ing/material	uncertainties	[63] SIMP properties of can-	didate materials	[67] SIMP material and geo-	metry	[111] SIMP material property	uncertainty manu-	facturing errors		[57] MMB loading orienta-	tion		[37] LST Loading uncer-	tainties	[66] SIMP Material proper-	ties	[64] SIMP Loading and geo-	metry uncertain-	ties	[132] SIMP Material prop-	FE 2 erties Loading	uncertianties

• To develop a density-based multiscale topology optimization. It enables one to embed the homogenization model for optimizing the macro topology and micro material distribution. 1: Literature of robust topology optimization

•

  To develop a multiscale robust topology optimization. Once we have tackled scale separation, uncertainty quantification, and multiscale topology optimization, the robust optimization is naturally solved.

  the Equations of 2.21, 2.27, 2.28, 2.29 and applying the average operators, the macroscopic thermal problem finally reads:

Table 3 .

 3 2: Uniform distribution of input micro parameters, and some characteristics of the output macroscopic quantities. k S , k F , n, C F , C S are random inputs, l f , L are deterministic inputs, C M , k M are outputs of uniform random variables, for which some analytical expressions are available,

	variables	uniform distribution interval	mean	coefficient of variation
	k S / (daW/m/K)	[0.5, 2]	1.25	0.346
	k F / (daW/m/K)	[1, 3]	2	0.289
	n	[0.3539, 0.3981]	0.376	0.034
	C F / (MJ/m 3 /K)	[1, 2]	1.5	0.192
	C S / (MJ/m 3 /K)	[2.5, 2.8]	2.65	0.033
	l F / m L / m	--	0.8 × 2 1	0 0
	C M / (MJ/m 3 /K)	-	2.22	0.055
	k M / (daW/m/K)	-	1.49	0.236

Table 3 .

 3 . Indeed, building costless reliable error indicators is a current field of interest to assess the accuracy of such fitted models. 3: Mean and std of β 0

		p=1	p=2	p=3
	mean	PCE 1D 0.30405 0.30405 0.3045 PCE 3D 0.45608 0.45608 0.45608
	std	PCE 1D 0.0188 PCE 3D 0.0895	0.0188 0.0896	0.0188 0.0896

Table 4 .

 4 2: Thermal mechanical material properties of candidate materials

Table 4 .

 4

	34)

Table 4 .

 4

	.16) assembled via graded

5: Material properties provided in

[START_REF] Gao | Topology optimization of thermo-elastic structures with multiple materials under mass constraint[END_REF]

Table 5 .

 5 1: µ c and σ c w.r.t. different λ

	Mean, µ c	79.69	78.244 79.88 80.61 119.94
	Standard deviation, σ c 14.80	15.36	13.81 11.66 0.6089

Table 5 .

 5 2: Optimization input parameters used (brackets indicate continuation )The prescribed deterministic material properties and loading are reported in table 5.3. Note that, they not real physical materials. For comparison, individual uncertainty is considered in RTO first. Then, combination of the 3 uncertainties on the RTO is trivial, but

		mesh size nelx × nely	penalty p	filter radius	max. iter. nb.
	UNI	80 × 40	5	[2,1]	[ 200,20]
	PRM [20 × 10, 40 × 20, 80 × 40, 80 × 40]	[3, 3, 3, 5]	[2,2,2,1]	[50,50,100,30]
	CP	80 × 40	[1.5,2,2.5,3,3] [2,2,2,2,1] [100,100,100,100,20]

Table 5 .

 5 3: Uniform distribution of uncertainty sources

Table 5 . 5

 55

	-5.10

Table 5 .

 5 6: RTO-CP-Q.

	λ	1	4	8	12
	Solutions
	µ	9.2	10.34	12.61	14.99
	σ	1.2	0.064	0.02	0.009

Table 5 . 7

 57 Table 5.8 reports different RTO-PRM designs w.r.t. varying λ = 1, 4, 5, 10. One has achieved identical pattern of macro topologies. Nevertheless, referring to figure 5.8, the different designs, show shape variation with different λ. It indicates that, in context of mechanical uncertainty, RTO-PRM with λ increasing, in fact, performs a shape optimization on a baseline macro topology, and allows to reduce the σ. • Table 5.10 reports the different RTO-UNI designs w.r.t. varying λ = 1, 2, 4, 8. Figure 5.9 displays the outlines of different designs. We have a similar pattern macro topology, however, in the same table, a RTO case with * mark, where one increases the uncertainty range σθ = [-π/9, π/9] shows a different topology. It indicates that both the λ and the interval range of the uncertainties affect the RTO design in macro topology, micro material distribution, mean, and standard deviation.

: RTO-UNI-Q 108 •

Table 5 . 8

 58

	λ	1	4	6	10
	Solutions
	µ	10.81	10.014	10.226	10.228
	σ	0.232	0.11	0.076	0.052

: RTO-PRM-θ

Table 5 . 9

 59

	λ	1	2	4	8	8 *
	Solutions
	µ	11.195	10.713	10.572	10.72	13.93
	σ	0.4308	0.388	0.303	0.2	0.0033

: RTO-CP-θ

Table 5 .

 5

	λ	0.1	1	5
	Solutions
	µ	11.103	10.352	9.517
	σ	1.272	1.165	1.04
	m f	68.59%	76.69%	81.87%

10: RTO-UNI-θ. * represents a RTO case with increasing uncertainty range of σ θ = [-π/9, π/9] 110 Figure 5.9: RTO-UNI-θ compliance c with different λ to RTO(UNI, CP, PRM); the standard deviation seems less sensitive to the increase of weight factor λ. And with given λ = 6, RTO-PRM does not manage to obtain a converged solution. Ref [25, 85, 131] have analyzed the effect of the random material field to the robust optimization (however, they are only macro topology optimization problem), which is still feasible to be introduced into our proposed RTO solution framework. It is only required to discretize the random field by series expansion method [104] (PCE and other methods), the remaining parts of stochastic process and optimization in RTO are consistent with our current implementation.

Table 5 .

 5

11: RTO-PRM-E m

Table 5 .

 5

		0.5	1	1.88	2	3
	Solutions
	µ	10.279	8.205	9.2	10.069	10.97
	σ	1.081	1.494	0.061	0.017	0.013

14: RTO design with hybrid uncertainties and different λ, solved by RTO-CP

  1 depicted geometry of inclusion with varying t. With given phases material proprietress in table 7.2, and parametric cells with varying t (shown in fig 7.1), tables report the fit error of independent elastic tensor components, E 11 (table 7.1), E 12 (7.2), E 44 (7.3), and thermal expansion α H (7.4) with respect to different PCE order p. With p = 7, the PCE can construct an accurate thermo-mechanical surrogate model to the FEM-based one, suitable for the multiscale thermo-mechanical coupled topology optimization.

	Parameters	mean	uniform interval
	E s (Pa)	2 × 10 11	[1.25, 2.75] × 10 10
	E f (Pa)	5 × 10 10	[2, 8] × 10 10

Table 7 .

 7 9, 0.9] PCE order p mean of E 11 /E f std ofE 11 /E f fit error % 1: QoI,E 11 /E f 129 PCE order p mean of E 12 /E f std of E 12 /E f fit error(%)

	1	1.9067	0.85101	11.613
	2	1.8943	1.0095	5.191
	3	1.898	1.0399	3.148
	4	1.9054	1.0628	1.030
	5	1.9004	1.0398	0.66

Table 7 .

 7 2: QoI,E 12 /E f PCE order p mean of E 44 /E f std of E 44 /E f fit error(%)

	1	0.45169	0.23419	13.64
	2	0.44517	0.27892	6.093
	3	0.44554	0.28982	4.1552
	4	0.44883	0.29609	1.4705
	5	0.44591	0.28925	0.8552

Table 7 .

 7 3: QoI,E 44 /E f PCE order p mean of α H /α f std of α H /α f fit error(%)

	2	2.371	0.90229	30.429
	3	2.3278	1.0936	15.526
	4	2.3323	1.0821	7.5374
	5	2.3328	1.0598	4.4332
	6	2.3379	1.0681	2.6365

Table 7 .

 7 4: QoI,α H /α f

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

Note herein, p indicates two different variables, and it is the penalty factor in the context of topology optimization, and the PCE order in which of PCE-based stochastic collocation points method. They are the classical nomenclature in their contexts, respectively. Thus the author wants to keep both, and the reader can tell the difference dependent on the context.

Acknowledgements

. Thanks to Mainak Bhattacharyya, Nawfal Blal for sharing a very enjoyable collaboration. Thanks to the committee members who played an integral role in my doctoral research. Thanks also to the CSC-UT/ISNA funding.

Thanks to my friends Ting SHI,