.2, the diffusion tensor is applied to obtain the heat flux at the central

+(a 2 *D i+ 1 2 ;j+ 1 2

confusions over the next three years of research. Secondly, I would like to thank my cosupervisor Gaeton Hello for his hard work over the past three years. I have also learnt a lot from him during the countless discussions. Not only scientific knowledge, but also a positive attitude towards research and life. Then I would like to thank my senior, Mu He, and my classmate, Xudong Wang, for our joint efforts to complete the analysis of the thermal conductivity of sintered silver materials and the comparison of different methods.

I would also like to thank the other teachers and colleagues in my lab, such as Cong Yu and Hu Libang, for your valuable advice when I was confused. In addition, thanks to Professors Mabrouk BEN TAHAR and Yuming LI who kindly agreed to review this thesis.

Heterogeneous materials

A heterogeneous material is one that is composed of sub domains with different materials (phases), such as a composite, or the same material in different states, such as a polycrystal. According to the review article [START_REF] Bargmann | Generation of 3d representative volume elements for heterogeneous materials: A review[END_REF], heterogeneous materials can be divided into porous and non-porous materials.

Porous materials at their turn can be divided into three main groups: fabrics, agglomerates, and aggregates. Aggregates (see Fig. 1.1-(a)) are composed of stacked individual particles that are surrounded by a network of inter-particle spaces. Aggregates encompass a diverse range of coarse particle materials, such as sand, rice, gravel, crushed stone, and slag. Agglomerates are collections of solid primary particles of any size or shape that are kept together by adhesive and/or cohesive forces. Agglomerates are by definition multiphase materials, containing at least one fluid phase in the interstitial volume between main particles [START_REF] Bika | Mechanical properties of agglomerates[END_REF]. Cellular materials can be described as a type of agglomerate. Cellular materials consist mainly of polyhedral cells with solid edges or faces [4]. Cellular materials are further divided into open-cell structures (see Fig. depending on the direction of processing. 2D weave structures can be divided into plain, twill and satin weave. 3D weave structures can be divided into interlaced or interlocked fabrics. Ansar et al. [START_REF] Ansar | Modeling strategies of 3d woven composites: A review[END_REF] gives an review of modelling of 3D woven composites. Non-woven structures are made by mechanically, chemically or thermally bonding their fibres into a random web structure. By changing the material, the bonding mechanism or the bonding point and fibre density, the non-wovens can be tailored to different requirements and thus used in a variety of applications.

Non-porous materials are further classified as polycrystals, bicontinuous composites, and matrix-inclusion materials. A polycrystal (see Fig. 1.1-(f)) is made up of numeros crystallites with varying crystal orientations separated by grain boundaries. Metals, their alloys, ceramic materials, rocks, and ice are mostly polycrystalline. Bicontinuous composites (see Fig. 1.1-(e)) comprise two solid phases, each spatially continuous and forming an interpenetrating network, such as carbon/silica composite [START_REF] Yang | Preparation and lithiumstorage performance of carbon/silica composite with a unique porous bicontinuous nanostructure[END_REF] or nanoporous titanium-based composites [START_REF] Lee | Periodic bicontinuous composites for high specific energy absorption[END_REF]. Matrix-inclusion composites (see Fig. ) ceramics after thermal etching [START_REF] Ikesue | Polycrystalline nd: Yag ceramics lasers[END_REF] (g) Carbon epoxy bicontinuous composite [START_REF] Lee | Periodic bicontinuous composites for high specific energy absorption[END_REF].

(h) Fiber-reinforced composite [START_REF] Bargmann | Generation of 3d representative volume elements for heterogeneous materials: A review[END_REF]

Homogenization methods

Materials exhibit different properties for different physical contexts, including mechanical [START_REF] Pelleg | Mechanical properties of materials[END_REF], thermal [START_REF] Buck | Thermal Properties[END_REF][START_REF] Cverna | ASM ready reference: thermal properties of metals[END_REF], and electrical [START_REF] Solymar | Electrical properties of materials[END_REF] properties.

Material's mechanical properties define the material's behavior under the action of a loading external force, generally including elasticity, strength, plasticity, toughness.

Thermal characteristics refer to a material's response to temperature variations and heat application. This response may take the form of an increase in temperature, a phase change, a change in length or volume, the commencement of a chemical reaction, or some other change in physical or chemical properties. They include thermal conductivity, heat capacity, and thermal expansion (thermo mechanical property).

Electrical properties reflect materials ability to conduct electricity. Resistivity, electrical conductivity are all examples of electrical properties.

The main purpose of homogenization is to determine the <<equivalent>> properties of heterogeneous materials based on the knowledge of their components properties and spacial arrangement. This makes it possible to replace the composite medium with an equivalent homogeneous medium at the macroscopic level.

In general, there are two classes of methods to estimate the equivalent properties of materials: analytical and numerical ones. The analytical methods were the first homogenization approaches to be proposed, leading to several expressions to estimate equivalent properties of a heterogeneous medium. Mean field methods are one type of analytical methods.

These methods consider the immersion of inclusions with simple geometries into an infinite matrix or effective medium and require simple microstructural informations such as the volume fraction and mean morphology of the inclusions [START_REF] Bargmann | Generation of 3d representative volume elements for heterogeneous materials: A review[END_REF].

Among the mean-field homogenization techniques, a distinction can be made between the bounding method and the estimated mean-field method [START_REF] Müller | Homogenization of linear elastic properties of short-fiber reinforced composites -a comparison of mean field and voxel-based methods[END_REF]. The bounding method generally gives an acceptable range of valid attributes. For example: Voigt [START_REF] Voigt | Ueber die beziehung zwischen den beiden elasticitätsconstanten isotroper körper[END_REF] and Reuss [START_REF] Reuss | Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle[END_REF] estimates. They determine the upper and lower bounds of a material's equivalent properties. Hashin and Shtrikman [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF] used the variational theorem to derive bounds on the equivalent properties of macroscopically homogeneous and isotropic two-phase materials. For other related models, a review is presented by Pietrak et al. [START_REF] Pietrak | A review of models for effective thermal conductivity of composite materials[END_REF]. The second class of method gives approximate valid attributes. For example: Maxwell model [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF] or Maxwell-Eucken model [START_REF] Eucken | Die wärmeleitfähigkeit keramischer feuerfester stoffe : ihre berechnung aus der wärmeleitfähigkeit der bestandteile[END_REF], the Mori-Tanaka mean field [START_REF] Mori | Average stress in matrix and average elastic energy of materials with misfitting inclusions[END_REF][START_REF] Benveniste | A new approach to the application of mori-tanaka's theory in composite materials[END_REF], the self-consistent method [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF][START_REF] Kröner | Berechnung der elastischen konstanten des vielkristalls aus den konstanten des einkristalls[END_REF][START_REF] Hill | A self-consistent mechanics of composite materials[END_REF]. The Mori-Tanaka mean field is achieved by averaging the effective fields acting on each heterogeneity as a field on a matrix. The self-consistent method was originally proposed for use with crystalline aggregates [START_REF] Hershey | The elasticity of an isotropic aggregate of anisotropic cubic crystals[END_REF][START_REF] Kröner | Berechnung der elastischen konstanten des vielkristalls aus den konstanten des einkristalls[END_REF]. It was later extended by Hill [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] for multiphase materials. Many researchers have also successfully applied mean-field methods to heterogeneous materials such as porous materials [START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF][START_REF] Sharma | Size-dependent eshelby's tensor for embedded nanoinclusions incorporating surface/interface energies[END_REF], particulate composites [START_REF] Hashin | Assessment of the self consistent scheme approximation: conductivity of particulate composites[END_REF][START_REF] Garboczi | Elastic moduli of a material containing composite inclusions: effective medium theory and finite element computations[END_REF], fibre-reinforced composites [START_REF] Duschlbauer | Computational simulation of composites reinforced by planar random fibers: homogenization and localization by unit cell and mean field approaches[END_REF][START_REF] Chou | A self-consistent approach to the elastic stiffness of short-fiber composites[END_REF][START_REF] Pipes | Self-consistent properties of carbon nanotubes and hexagonal arrays as composite reinforcements[END_REF] and polycrystals [START_REF] Lebensohn | A self-consistent anisotropic approach for the simulation of plastic deformation and texture development of polycrystals: application to zirconium alloys[END_REF][START_REF] Molinari | A self consistent approach of the large deformation polycrystal viscoplasticity[END_REF].

Typically, Mean field methods are only accurate for assessing equivalent properties when the volume fraction of the inclusion phase is small. For inclusions, phases are limited to simple mathematical morphologies, i.e., spheres, ellipsoids, etc [START_REF] Bargmann | Generation of 3d representative volume elements for heterogeneous materials: A review[END_REF]. As a result, analytical methods are insufficient for estimating the equivalent properties of materials with complex microstructures and morphologies, such as nanostructured materials or advanced composites. Numerical methods are therefore very helpful in practically obtaining the equivalent properties. With the computational homogenization methods, we need to consider a characteristic sample of the material whose homogenized properties are expected to coincide with the material's ones. Since computation of large volumes is often prohibitive, we need to create a model sufficiently large to provide accurate results and small enough so that computations can be carried out with reasonable time and memory, such a model is called Representative Volume Element (RVE) [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF]. In order to satisfy the assumption that the average properties of the RVE are equal to the average properties of a given heterogeneous material, the RVE must contain sufficient microstructural information. Furthermore the RVE is the smallest equivalent material volume element in the composite [START_REF] Drugan | A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites[END_REF]. Currently, computational homogenization methods have become powerful approaches to determine the properties of heterogeneous materials due to the relative speed of its calculation, the accuracy of its results, and its application to potentially any heterogeneous material. Such computational methods include Finite Element Method, Fast Fourier Transform, Boundary Element Method, Finite Different Method, Phase Field, etc.

Computational homogenization of heterogeneous materials 1.2.1 Representative models for the material architecture

The structure of a material is generally obtained from design assumptions, the observation of manufactured samples or sample computations. Once the material structure has been determined, a digital representative model needs to be built for further study.

CAD models

There are many available computational tools to generate CAD models for various kinds of heterogeneous materials. Popular in the composites world are WiseTex [START_REF] Verpoest | Virtual textile composites software wisetex: Integration with micro-mechanical, permeability and structural analysis[END_REF], TexGen [START_REF] Stig | Spatial modelling of 3d-woven textiles[END_REF], and Texmind braider [START_REF] Kyosev | Geometrical modeling and computational mechanics tools for braided structures[END_REF]. These are powerful and efficient softwares that can handle textile reinforcement of any structure, including unidirectional (UD), non-crimped fabrics (NCF) and 3D, and all feature cell meshing [START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF]. In addition to the tools mentioned above, several researchers have proposed methods for CAD modelling. For example, Wedling et al. [START_REF] Wendling | Consistent geometrical modelling of interlock fabrics[END_REF] proposed a novel CAD modeling method of the geometry of interlock fabrics unit cells.

Chen et al. [START_REF] Chen | Cad/cam for complex woven fabrics. part ii: Multi-layer fabrics[END_REF] also propose a CAD modelling approach for Complex Woven Fabrics. For CAD modelling of polycrystalline materials, Neper [START_REF] Quey | Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF] is a powerful tool for polycrystal generation and meshing. Some researchers have also proposed some CAD modelling for polycrystalline materials, such as Ghosh et al. [START_REF] Ghosh | Cad-based reconstruction of 3d polycrystalline alloy microstructures from fib generated serial sections[END_REF] provided a robust CAD-based approach for modelling the three-dimensional microstructures of polycrystalline metals using crys-tallographic data acquired using a focused ion beam (FIB)-scanning electron microscopy (SEM) system. For porous materials, Kou et al. [START_REF] Kou | A simple and effective geometric representation for irregular porous structure modeling[END_REF][START_REF] St Kou | An approach of irregular porous structure modeling based on subdivision and nurbs[END_REF] proposed a CAD-based approach for irregular porous structures. There are other software such as Surface Evolver [START_REF] Brakke | The surface evolver[END_REF] for the study of surfaces formed by surface tension and other energies, which can handle arbitrary topologies, volume constraints, boundary constraints, boundary contact angles, crystal quilt functions, gravity, etc. As can be seen from the brief review above, most constructive CAD modelling methods are able to represent only some of the various architectures that can be found in real world materials.

Models based on conformal meshes

In a volume mesh-based model, a heterogeneous material is represented as a collection of polyhedra. They are conformal in the sense that the boundaries of all elements are expected to represent faithfully the interfaces between the different phases involved. For meshing, the following widely used methods have been proposed by several researchers: the model meshed by the Delaunay triangulation method [START_REF] Ruppert | A delaunay refinement algorithm for quality 2-dimensional mesh generation[END_REF][START_REF] Shewchuk | Triangle: Engineering a 2d quality mesh generator and delaunay triangulator[END_REF] shown in Fig. 1.2-(a). [START_REF] Yerry | A modified quadtree approach to finite element mesh generation[END_REF][START_REF] Yerry | Automatic three-dimensional mesh generation by the modified-octree technique[END_REF].

Lo [START_REF] Lo | Finite element mesh generation and adaptive meshing[END_REF] provided a comprehensive review of the above common mesh generation techniques for general applications. These most frequently used approaches prioritize geometric and topological compatibility, with geometric compatibility guaranteeing that the mesh roughly matches the geometry of the target model and topological compatibility ensuring that the elements have the proper adjacency [START_REF] Sullivan | A three-dimensional mesh generator for arbitrary multiple material domains[END_REF]. However, for some heterogeneous materials, traditional meshing methods are often difficult to implement, such as 3D woven fiber structures, where it is often not feasible to avoid distorted elements when using traditional meshing methods, and a lot of time has to be invested in handcrafting high quality meshes [START_REF] Potter | Mesh generation and geometrical modelling of 3d woven composites with variable tow cross-sections[END_REF]. Therefore some researchers have proposed different meshing methods for different heterogeneous materials, such as Grail et al. [START_REF] Grail | Consistent finite element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements[END_REF] who created a new automated method for the mesoscale representation of textile composites unit cells to generate smooth Finite and conformal discretization of 3D woven composites RVEs. Rassineux [START_REF] Rassineux | Robust conformal adaptive meshing of complex textile composites unit cells[END_REF] presented a methodology to create conformal tetrahedral meshes of realistic complex woven unit cells.

In the field of polycrystalline materials, Quey [START_REF] Quey | Large-scale 3d random polycrystals for the finite element method: Generation, meshing and remeshing[END_REF] proposed a method for the generation and meshing of large-scale three-dimensional random polycrystals. Resk [START_REF] Resk | Adaptive mesh refinement in crystal plasticity finite element simulations of large deformations in polycrystalline aggregates[END_REF] developed an adaptive mesh refinement and automatic reshaping scheme to study crystal plasticity with an anisotropic mesh refinement strategy based on the level set description. S2M is a well-established software for generating conformal meshes for polycrystalline materials.

Voxel based models

Voxel based modeling is a way of representing three-dimensional objects on regular hexahedral cells. Due to their inherent discrete nature, they are naturally suited to scientific analysis and visualization involving volumes. Voxel modeling becomes more accurate as the size of voxels gets smaller and often requires relatively large amount of memory to handle the discretized data. With the increased affordability of computers with large memory, this defect tends to become less penalizing. More importantly, voxel models can effectively represent the complex internal geometrical details of solid models and allow to easily handle internal properties that are anisotropic or heterogeneous. Models that are difficult to realize in traditional modeling schemes, such as the human body [START_REF] Caon | Voxel-based computational models of real human anatomy: a review[END_REF] and internal organs, are easily realized with voxel models [START_REF] Patil | Voxel-based solid models: representation, display and geometric analysis[END_REF][START_REF] Patil | Voxel-based representation, display and thickness analysis of intricate shapes[END_REF]. Generally, x-ray tomography and magnetic resonance imaging equipment can be used to obtain the data. Overall, the advantages of voxel modeling include simplicity, intuitive representation, identical complexity for all objects, model heterogeneity that can be easily incorporated into the analysis. Disadvan-

tages

Computational homogenization methods

Once the model has been obtained by the above method, it can be calculated using a computational homogenisation method, such as Finite Element Method, Fast Fourier Transform, Boundary Element Method, Finite Different Method, Phase Field, etc.

Finite Element Method

The Finite Element Method (FEM) is a highly effective and commonly used numerical calculation method. FEM is based on the variational principle and the weighted residue method. The primary solution idea is to divide the computational domain into a finite number of non-overlapping cells. Within each cell, select some suitable nodes as interpolation points for the solution function. The PDE is solved discretely with the help of the variational principle or the weighted residual method.

The advantage of the Finite Element Method is that it is highly versatile, popular and well documented. There is a lot of very good available open source as well as commercial softwares. The disadvantage is, for example, that the model must be meshed in order to perform finite element calculations. And the cost of meshing is too high for some complex structures. The implicit form of the FEM, depending on the particular conditions and the solver used, has a time complexity of O(n 2 -n 3 ). This complexity also limits the use of very large meshes [START_REF] Lucarini | On the accuracy of spectral solvers for micromechanics based fatigue modeling[END_REF].

The Finite Element Method has been used by many researchers to calculate the properties of heterogeneous materials. For instance, Hollister and Kikuchi [START_REF] Hollister | A comparison of homogenization and standard mechanics analyses for periodic porous composites[END_REF] used Finite Element Method to compare homogenization and standard mechanics RVE based analyses for periodic porous composites. Hollister and Kikuchi [START_REF] Hollister | Homogenization theory and digital imaging: a basis for studying the mechanics and design principles of bone tissue[END_REF] provided a finite element homogenization method that can be adapted to analyse bone tissue. Smite et al. [START_REF] Smit | Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling[END_REF] proposed a Finite Element homogenization method to predict the mechanical behaviour of non-linear heterogeneous systems. Kanit et al. [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] determined the size of representative volume elements of the random composite and applies FEM to calculate the effective physical properties of the material. Fritzen et al. [START_REF] Fritzen | Computational homogenization of elasto-plastic porous metals[END_REF] investigated the effective material response of ductile metals with spherical pores with volume fractions between 0.1% and 30%. Dirrenberger et al. [START_REF] Dirrenberger | Effective elastic properties of auxetic microstructures: Anisotropy and structural applications[END_REF] calculated the elastic modulus of three auxetic periodic microstructures using Finite Element combined with periodic homogenization technique.

Khdir et al. [START_REF] Khdir | Computational homogenization of elastic-plastic composites[END_REF] discussed a method for estimating the effective elastic-plastic response of random two-phase composite media using a computer homogenization technique. Florez et al. [START_REF] Florez | Effective thermal conductivity of sintered porous media: Model and experimental validation[END_REF] investigated sintered porous media's equivalent thermal conductivity. They demonstrated that the geometry of the solid matrix in a porous medium has a relevant effect. El Moumen et al. [START_REF] Moumen | Computational thermal conductivity in porous materials using homogenization techniques: Numerical and statistical approaches[END_REF] calculated the equivalent thermal conductivity of porous materials at the microscopic scale using numerical homogenization techniques and morphological analysis. Vel et al. [START_REF] Vel | Computational homogenization and micromechanical analysis of textured polycrystalline materials[END_REF] presented a FEM based homogenization method for the thermomechanical analysis of real polycrystalline material microstructures obtained using electron backscatter diffraction techniques. Kaddouri et al. [START_REF] Kaddouri | On the effect of inclusion shape on effective thermal conductivity of heterogeneous materials[END_REF] estimated the effective thermal conductivity of random two-dimensional two-phase heterogeneous materials using a numerical homogenization technique. Sukiman et al. [START_REF] Sukiman | Effective thermal and mechanical properties of randomly oriented short and long fiber composites[END_REF] calculated the effective thermal and mechanical properties of randomly oriented short and long fiber composites. Signor et al. [START_REF] Signor | Evolution of the Thermal Conductivity of Sintered Silver Joints with their Porosity Predicted by the Finite Element Analysis of Real 3D Microstructures[END_REF] studied the thermal conductivity with the change of porosity for sintered silver by using Finite Element analysis in a realistic 3D microstructure. Torquato [START_REF] Torquato | Predicting transport characteristics of hyperuniform porous media via rigorous microstructure-property relations[END_REF] used Finite Element homogenization method for predicting transport properties in hyper homogeneous porous media. Qing et al. [START_REF] Qin | Crack Effect on the Equivalent Thermal Conductivity of Porously Sintered Silver[END_REF][START_REF] Qin | Evaluation of thermal conductivity for sintered silver considering aging effect with microstructure based model[END_REF] proposed a quantitative technique for evaluating the influence of crack evolution on the equivalent thermal conductivity of porous sintered silver, and also proposed a numerical approach based on specific microstructural features to calculate the equivalent thermal conductivity by taking into account the aging response.

Fast Fourier Transform

The Fast Fourier Transform (FFT) is a fast method for computing the discrete Fourier transform (DFT) using a computer [START_REF] Rao | Fast Fourier Transform -Algorithms and Applications[END_REF]. In particular, the more sampling points N that is transformed, the more significant the savings in the computational effort of the FFT algorithm, i.e., O(Nlog N) of time instead of O(N 2 ) of time [START_REF] Rao | Fast Fourier Transform -Algorithms and Applications[END_REF]. The principle of the Fast Fourier Transform algorithm is to implement large-scale transformations by means of many small and more easily performed transformations, reducing the operational requirements and increasing the speed with which they can be performed. The advantages of Fast Fourier Transform are its excellent numerical performance (the algorithm scales to n log n), its reduced computational cost compared to FE and its ability to increase the speed of simulations by several orders of magnitude [START_REF] Eisenlohr | A spectral method solution to crystal elasto-viscoplasticity at finite strains[END_REF]. The disadvantage of the FFT method is that its convergence rate and accuracy are highly dependent on the contrast between the phases represented in the domain [START_REF] Lucarini | Adaptation and validation of fft methods for homogenization of lattice based materials[END_REF]. It requires that the signal must be smooth.

In the real world, most of the signals we collect are non-stationary and non-linear. In the field of homogenization, the Fast Fourier Transform homogenization approach, first introduced by Moulinec and Suquet [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] for the purpose of studying composite materials, is a handy method for determining the effective properties of heterogeneous materials having periodic microstructures. Lebensohn et al. [START_REF] Lebensohn | N-site modeling of a 3d viscoplastic polycrystal using fast fourier transform[END_REF] first extended this FFT-based approach to crystal viscoplasticity. In a subsequent study, this FFT-based method was extended to predicting the elastoviscoplasticity of polycrystalline materials in two and three dimensions [START_REF] Lebensohn | Orientation imagebased micromechanical modelling of subgrain texture evolution in polycrystalline copper[END_REF][START_REF] Lebensohn | Study of the antiplane deformation of linear 2-d polycrystals with different microstructures[END_REF][START_REF] Lebensohn | Dilatational viscoplasticity of polycrystalline solids with intergranular cavities[END_REF][START_REF] Lebensohn | An elasto-viscoplastic formulation based on fast fourier transforms for the prediction of micromechanical fields in polycrystalline materials[END_REF][START_REF] Lebensohn | Modeling viscoplastic behavior and heterogeneous intracrystalline deformation of columnar ice polycrystals[END_REF]. Kabel et al. [START_REF] Kabel | Use of composite voxels in fft-based homogenization[END_REF] proposed a method which use of composite voxels in FFT-based homogenization. Li et al. [START_REF] Li | A non-local fracture model for composite laminates and numerical simulations by using the fft method[END_REF] applied the FFT-based homogenization method to composite laminates. Monchiet et al. [START_REF] Monchiet | Numerical homogenization of nonlinear composites with a polarization-based fft iterative scheme[END_REF] applied FFT-based homogenization methods to nonlinear composites.

Boundary Element Method

The Boundary Element Method (BEM) is a more accurate and efficient method developed after the Finite Element method. The most significant feature of boundary elements is that, in contrast to the Finite Difference or Finite Element methods, the methodology of formulating boundary value problems as boundary integral equations describes problems only by equations with known and unknown boundary states [START_REF] Antes | A short course on boundary element methods[END_REF]. The Boundary Elements Method has become a powerful alternative to the Finite Element method, particularly where the domain extends to infinity. The most important feature of this method is that it only requires discretization of surfaces rather than volumes. As a result, the code of the Boundary Elements Method is easier to use with existing solid modelers and mesh generators. This advantage is particularly important for design, as the process usually involves a series of modifications, Boundary Elements Method meshes can be easily generated and design changes do not require complete re-division, which is more difficult to do with Finite Elements [START_REF] Brebbia | Boundary Element Techniques: Theory and Applications in Engineering[END_REF]. If the boundary is not smooth but contains corners and edges, or when the boundary conditions are discontinuous then the solution to the boundary value problem contains singularities at the boundary. In addition, pure BEM cannot solve problems with inhomogeneities or non-linear differential equations. These are the main difficulties with the BEM. More information on the principles, advantages, and disadvantages of the BEM can be found in this article [START_REF] Costabel | Principles of boundary element methods[END_REF].

On the other hand, there has been much research activity in the use of BEM to assess the equivalent properties of heterogeneous materials, such as the BEM-based homogenization method proposed by Okada et al. [START_REF] Okada | Homogenization method for heterogeneous material based on boundary element method[END_REF][START_REF] Okada | Homogenization analysis for particulate composite materials using the boundary element method[END_REF] for heterogeneous elastic materials with periodic microstructures and granular composites to obtain the corresponding effective mechanical properties, respectively. Galvis et al. [START_REF] Galvis | Analysis of three-dimensional hexagonal and cubic polycrystals using the boundary element method[END_REF] evaluated the macroscopic elastic properties of hexagonal closed package (HCP) and face centered cubic (FCC) polycrystals using the BEM. Procházka et al. [START_REF] Procházka | A bem formulation for homogenization of composites with randomly distributed fibers[END_REF] applied a BEM-based homogenization method to study composites with randomly distributed fibres. Kamiński [START_REF] Kamiński | Boundary element method homogenization of the periodic linear elastic fiber composites[END_REF] applied a BEM-based homogenization method to study periodic linear elastic fibre composites.

Finite Difference Method

The Finite Difference Method (FDM) was one of the first method used in computer numerical simulation of PDE and is still widely used today [START_REF] Mitchell | The Finite Difference Method in Partial Differential Equations[END_REF]. The method divides the solution domain into a difference grid, using a finite number of grid nodes instead of a continuous solution domain. The Finite Difference Method uses methods such as Taylor series expansion to discretize the derivatives in the control equations by replacing the dif-ference quotient of the function values at the grid nodes, thus creating a system of algebraic equations with the values at the grid nodes as the unknowns.

The advantages of the Finite Difference Method are the direct transformation of differential problems into algebraic problems, the absence of the need to construct form functions, the absence of unitary and integral analysis, the simplicity of mathematical modeling, the ease of programming, and the ease of parallelism. The disadvantage is that irregular areas are cumbersome to deal with, e.g. using jagged shapes to represent curved shapes.

Abudull and E [START_REF] Abdulle | Finite difference heterogeneous multi-scale method for homogenization problems[END_REF] developed the Finite Difference heterogeneous multiscale method (FD-HMM) for solving multiscale parabolic problems. The method is based on the heterogeneous multiscale method (HMM) and heterogeneous discretization, specifically focused on fine scales representing the size of tiny regions in the spatial domain. FD-HMM consists of two parts: a macroscopic solution generated on a coarse grid using the known data extracted from the microscopic model solution, and a microscopic solution that is solved from the original equations over a sparse (heterogeneous) spatial domain. Chen et al. [START_REF] Chen | Application of the finite difference heterogeneous multiscale method to the richards' equation[END_REF][START_REF] Chen | New scheme of finite difference heterogeneous multiscale method to solve saturated flow in porous media[END_REF] used FD-HMM to handle the issues associated with unsaturated water flow in random porous material and later improved the FD-HMM scheme to simulate not only steady saturated flow problems in geostatistical stochastic porous media, but also transient saturated flow problems. Jaworska et al. [START_REF] Jaworska | On two-scale analysis of heterogeneous materials by means of the meshless finite difference method[END_REF] used meshless Finite Difference method for two-scale analysis of heterogeneous materials. Orive et al. [START_REF] Orive | Finite difference approximation of homogenization problems for elliptic equations[END_REF] presented a method of the homogenization of the Finite Difference schemes approximating a family of elliptic equations with rapidly oscillating coefficients.

Other computational homogenization methods

In addition to the computational homogenisation methods described above there are other methods such as phase field, Generalized Method of Cells, High-Fidelity Generalized Method of Cells, and Variational Asymptotical Method for Unit Cell Homogenization.

The Phase field method is a mathematical approach to interfaces based on Ginzburg-Landau theory [START_REF] Cyrot | Ginzburg-landau theory for superconductors[END_REF]. One of its key ideas is to use the gradient of the order parameter to trace the interface, which has the advantage of not displaying the traced interface and saves a lot of effort in handling the interface grid. A disadvantage is that the phase field method requires prior knowledge of the material properties of the system under study, such as density, heat capacity and latent heat, among other properties. However, some properties are difficult to obtain experimentally such as interfaces and kinetic growth coefficients. Furthermore in the field of homogenization, the Voigt/Taylor model and the Reuss/Sachs model are the two classical homogenization assumptions explored so far in phase field theory, and Ammar et al. [START_REF] Ammar | Combining phase field approach and homogenization methods for modelling phase transformation in elastoplastic media[END_REF] then first combined the phase field approach and the homogenization method to model the phase transformation in elastoplastic media.

The generalized method of cells (GMC) [START_REF] Aboudi | Micromechanical analysis of composites by the method of cells-update[END_REF] is a semi-analytical technique for determining the effective characteristics of composites. This approach discretizes a representative volume element (RVE) of the composite into a regular grid of subcells. Integrals over subcell boundaries are used to satisfy equilibrium and compatibility on an average basis across subcells. One advantage of GMC over other numerical techniques is that it can calculate the entire set of a composite's effective elastic properties in a single step rather than solving a series of boundary value problems with varying boundary conditions [START_REF] Banerjee | On the generalized method of cells and the prediction of effective elastic properties of polymer bonded explosives[END_REF]. In the field of homogenization, this method has also been used by several researchers, such as Greengard et al. [START_REF] Greengard | On the numerical evaluation of elastostatic fields in locally isotropic two-dimensional composites[END_REF] who used this method for the numerical evaluation of the elastic field of two-phase composites. Moghaddam [START_REF] Ghorbani Moghaddam | A multi-scale computational model using generalized method of cells (gmc) homogenization for multi-phase single crystal metals[END_REF] created a multi-scale computational model using Generalized Method of Cells (GMC) homogenization for multi-phase single crystal metals.

The High-Fidelity Generalized Method of Cells (HFGMC) is a micromechanics technique that can be used to simulate nonlinear composite materials [START_REF] Aboudi | Micromechanics of composite materials: a generalized multiscale analysis approach[END_REF]. In contrast to the Generalised Method of Cells, HFGMC provides the necessary shear-coupling between the local normal and shear deformation fields and the macroscopically applied average strains [START_REF] Aboudi | High-fidelity generalized method of cells for inelastic periodic multiphase materials[END_REF]. This method's core computing effort is spent solving sets of simultaneous linear algebraic equations repeatedly in order to determine local/global field quantities and effective qualities for heterogeneous materials with a periodic microstructure [START_REF] Ricks | Solution of the nonlinear high-fidelity generalized method of cells 245 micromechanics relations via order-reduction techniques[END_REF]. The advantage of HFGMC is that the microstructure of the composite can be selectively geometrically refined as few number of rectangular subcells is sufficient to generate the effective linear and nonlinear response for the overall composite with high accuracy [START_REF] Haj-Ali | Nonlinear micromechanical formulation of the high fidelity generalized method of cells[END_REF]. The disadvantage is that when dealing with non-linear materials, the computational running time becomes too long, especially when representing more detailed microstructures [START_REF] Ricks | Solution of the nonlinear high-fidelity generalized method of cells 245 micromechanics relations via order-reduction techniques[END_REF]. HFGMC has been performed to model different physical effects in periodic composites, such as viscoelasticviscoplastic micromechanics [START_REF] Aboudi | Micromechanically established constitutive equations for multiphase materials with viscoelastic-viscoplastic phases[END_REF], damage evolution effects [START_REF] Aboudi | The effect of anisotropic damage evolution on the behavior of ductile and brittle matrix composites[END_REF], etc.

Variational Asymptotical Method for Unit Cell Homogenization (VAMUCH) is a generalpurpose micromechanics algorithm based on Finite Elements. It is based on the variational asymptotic method [START_REF] Vl Berdichevskii | Variational-asymptotic method of constructing a theory of shells[END_REF] for homogenising heterogeneous materials. It can be used to calculate the material properties of arbitrary microstructures of arbitrary heterogeneous materials. The advantages of this method are that it has inherent variational properties, can be solved directly using numerical methods and, secondly, can calculate material properties in different directions simultaneously, making it more efficient than the Finite Element micromechanics method which requires repeated runs under different loading conditions. Variational Asymptotical Method for Unit Cell Homogenization approach was first proposed by Yu et al. [START_REF] Yu | Variational asymptotic method for unit cell homogenization of periodically heterogeneous materials[END_REF] for periodic heterogeneous materials and has since been developed for predicting the thermoelastic properties of heterogeneous materials [START_REF] Yu | A variational asymptotic micromechanics model for predicting thermoelastic properties of heterogeneous materials[END_REF] and for predicting the conductivity of composites [START_REF] Tang | A variational asymptotic micromechanics model for predicting conductivities of composite materials[END_REF]. It has since been developed by other researchers to predict the effective properties of heterogeneous materials, such as elastic [START_REF] Bouhala | An advanced numerical method for predicting effective elastic properties of heterogeneous composite materials[END_REF], elasto-viscoplastic [START_REF] Zhang | A micromechanics approach to homogenizing elasto-viscoplastic heterogeneous materials[END_REF], hyperelastic [START_REF] Zhang | A variational asymptotic method for unit cell homogenization of hyperelastic heterogeneous materials[END_REF], etc. For the GMC, HFGMC, VAMUCH and other methods mentioned above, YU et al. [START_REF] Yu | A critical evaluation of the predictive capabilities of various advanced micromechanics models[END_REF] has done a critical assessment on the predictive power of the micromechanical models. stage sintering. In addition to this, voxel has also shown great competence in the geometric modelling of numerical models. For example, Hello et al. [START_REF] Hello | Numerical simulations of woven composite materials with voxel-fe models[END_REF] used voxel-FE models for the numerical simulation of any kind of material architectures (polycrystalline aggregates, woven composites, particles reinforced composites). Therefore, the use of voxels as a representative model for heterogeneous materials is definitely the best choice for this work.

The commonly used Finite Element Method is not the best choice when using voxelized mesh as a premise. This is because the Finite Element method does not take advantage of voxel discretization. In terms of convergence speed, the use of voxel is not as fast as with conformal mesh. There is no advantage in terms of computation time for large problems, and the sparse matrices generated have more non-zero terms and are more memory intensive. For other methods, the FFT method presupposes periodicity, i.e. the microstructure is regarded as the fundamental unit of an infinitely periodic medium [START_REF] Willot | Numerical homogenization of random media: the fft method[END_REF], which strongly limits the use of the method since the internal structure of most real materials is random and not periodic. For the BEM method, it deals well with isotropic problems, but is too complicated for anisotropic problems. The solution matrix resulting from the Boundary Element formulation is unsymmetric and fully populated with non-zero coefficients. Although the Boundary Element method requires only a relatively small number of nodes and elements to achieve the same accuracy as the Finite Element solution, it still occupies a large amount of memory when faced with very large problems. So the Finite Difference method is definitely the more suitable method for this work, because of its natural suitability for voxelized mesh. Secondly, with the enhancement of computing power of computers and owing to the natural advantage of FDM in parallel computing, the proposed method is scalable in terms of computational efficiency. Moreover, depending on the different Finite Difference models, it can be made to deal with very large problems with far fewer non-zero terms in the sparse matrix than the FEM.

Heterogeneous materials are increasingly used for their superior overall properties, such as porous media, which are widely used in the electronics and biomedical industries, so determining the equivalent thermal conductivity (ETC) of heterogeneous materials is essential for the correct design of industrial equipment that may be subjected to severe thermal loads during use. The high cost of using the Finite Element Method in the face of some large problems has prompted the search for alternatives to the Finite Element Method.

Therefore, in this thesis, the thermal problem will be taken as the central problem and the Finite Difference Method will be used to deal with the anisotropic diffusion problem. The main objective of the paper is to present a numerical method and software for calculating the equivalent thermal conductivity of heterogeneous material is proposed, based on a voxelized mesh, using Finite Difference Method (FDM).

Outline of the thesis

This thesis is organized as follows:

In Chapter 2, the Günter scheme, a Finite Difference model applicable to anisotropic diffusion, is introduced and a three-dimensional extension of the model is given. An overview of homogenisation theory is also given and the approach of some researchers to loading periodic and mixed homogeneous boundary conditions is summarised.

In Chapter 3, a 5-point model and an integral model (2D) inspired by Gunther's model is developed. The method of loading periodic and mixed uniform boundary conditions for the Günter scheme is also presented.

In Chapter 4, results are shown comparing the Finite Difference method with the Finite Element method and Finite Element + pixel (voxel) for different RVEs where the shapes of inclusions in the RVEs include crosses, circles, ellipses, spheres and cylinders. Applications to sintered silver materials are also shown.

Chapter 5 focuses on some complementary studies, including the dynamic model of our method and the ADI method, and the thesis concludes with a summary and discussion in Chapter 6.

35

Chapter 2

Theoretical background 2.1 Introduction

Following on from the previous chapter where we introduced the different numerical methods for computational homogenization and the main purpose of the thesis, in this chapter, the focus will be on the theoretical background of the thesis. First of all, the thermal problem will be the main object of study of the thesis, as it is more suitable as a beginning of a new approach. Therefore, in this chapter, the main physical parameters in the thermal problem will be presented. The equivalent thermal conductivity of heterogeneous and complex materials is investigated by introducing homogenisation theory. To achieve the above, the Günter scheme will be used as the central finite difference model for this thesis, as it deals well with discontinuous thermal conductivity problems and shows great stability and sufficient accuracy in simulating anisotropic diffusion. Therefore, in this chapter, the two-dimensional formulation of the Günter scheme for isotropic diffusion as well as for anisotropic diffusion is re-demonstrated. The corresponding three-dimensional model extensions and homogenisation theories are also given. The different topics detailed are as follow:

Section 2.1 introduces the basic formulation of the thermal problem, the physical parameters used and the isotropic and anisotropic thermal conductivity tensor that we use. Section 2.2 re-demonstrates Günter's symmetric two-dimensional model and gives the extension of the three-dimensional model. Section 2.3 focus on the homogenization method we use and the choice of boundary conditions.

Thermal problems

In this work, the main objective is to study the equivalent properties of heterogeneous materials. Then the finite difference method is chosen as the core of the numerical approach and the anisotropic diffusion problem is taken as the main object of study. For diffusion problems, it has a wide range of applications, including thermal problems, plasma problems, image processing, etc. First, this work will consider the steady-state thermal problem.

Fields

Temperature field

Temperature field (T ) is a collection of temperatures at various points in the object at various times. In general the temperature field of an object is a function of coordinates and time, i.e.

T = f (x; y ; z; t) (2.1)
The temperature field can be divided into two categories depending on whether it varies with time: steady state temperature and non-steady state temperature. The unit is Kelvin (K).

Temperature Gradient

Temperature Gradient (∇T ) is the rate of change of temperature with position in a given @T @x @T @y @T @z 3 7 7 7 7 5

(2.2)

Heat flux

The heat flux (q) is the heat energy per unit of time that flows through a surface of unit area. It is a vector field, its unit in the international system of units is Watt per square metre (W=(m 2 )).

Thermal conductivity

Thermal conductivity is defined mathematically by Fourier's law. It is a parameter that indicates a substance's capacity to transmit heat through a medium.

Thermal conductivity tensor

In the study of heat conduction in solids, isotropic thermal conductivity tensor is often used, that is, the thermal conductivity does not depend on the direction. i.e. (2.3) However, there are numerous natural and synthetic materials whose thermal conductivity varies with direction, forming anisotropic materials. For example, thermal conductivity in wood varies along the grain, through the grain, and circumferentially, while thermal conductivity in laminates varies across and along the laminate.

D =
Crystals can be classified according to their symmetry into seven main systems: triclinic, monoclinic, orthorhombic, hexagonal, tetragonal, trigonal, and cubic systems [START_REF]Heat Conduction in Anisotropic Solids[END_REF]. The cubic system is same as isotropic thermal conductivity tensor. In triclinic system, conductivity coefficients are not constrained by symmetry considerations. There are nine components and they can be non-zero, and we have (2.4)

D = 2
In monoclinic system, due to symmetry considerations, some of the components become zero. Hence, we obtain (2.5)

D = 2
In orthorhombic system, the conductivity coefficients are given as follows:

D = 2 6 6 6 6 4
D 11 0 0 0 D 22 0 0 0 D 33 3 7 7 7 7 5
(2.6)

In trigonal, tetragonal, and hexagonal system, we have (2.8)

D =
Onsager's [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] principles of irreversible thermodynamics demonstrate that when the fluxes (i.e., q) and forces (i.e., ∇T ) are linearly linked to one another as indicated by Eq.2.13, the phenomenological coefficients follow the reciprocity relation. Casimir [START_REF] Casimir | On onsager's principle of microscopic reversibility[END_REF] discusses how Onsager's reciprocity relation is applied to the thermal conductivity coefficients related with heat conduction in anisotropic materials. As a result, it is reasonable to assume that the conductivity coefficients D ij obey the reciprocity relation

D ij = D ji i; j = 1; 2; 3 (2.9)
Additionally, as noted in [START_REF] Prigogine | Introduction to Thermodynamics of Irreversible Processes[END_REF], irreversible thermodynamics dictates that the coefficients D 11 , D 22 , and D 33 be positive, i.e.

D ii > 0 (2.10)
Hence the magnitude of the coefficients D ij is constrained by the requirement [START_REF] Onsager | Reciprocal relations in irreversible processes[END_REF] for

i = j. D ii D jj > D 2 ij F or i = 1; 2; 3 (2.11)

Examples of material's thermal conductivity

Below are the thermal conductivity of several common metallic isotropic materials and the woven composites.

T ( • C) (kg =m 3 ) D(W=(m • K)) D(W=(m • K)) C/C-Sic (2D) 25

Equations

With the steady diffusion equation, we have:

div ( - → q ) = Q (2.12)
where Q is a a heat source, q the heat flux. According to Fourier's law, in the case of isotropic diffusion, the heat flux can be expressed as:

- → q = -D • ∇T (2.13)
where T represents temperature, D the thermal conductivity independent of direction, and q the heat flux.

In the case of anisotropic solids, for a three-dimensional problem, the heat flux is composed of q 1 , q 2 , and q 3 along the directions in the rectangular coordinate system Ox, Oy, and Oz.

The general expression involving the heat flux, thermal gradient, and conductivity tensor is given by: [START_REF] Morel | A local support-operators diffusion discretization scheme for hexahedral meshes[END_REF][START_REF] Morel | A local support-operators diffusion discretization scheme for quadrilateralr-zmeshes[END_REF][START_REF] Morel | A cell-centered lagrangian-mesh diffusion differencing scheme[END_REF] used the MFDM to derive a cell-centered diffusion differencing scheme that gives a sparse matrix representation. In contrast, the traditional method of support operators gives a dense matrix representation. Günter et al. [START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF] offer two discretization techniques(Asymmetric and symmetric schemes) on rectangular grids that employ the finite-difference method and SOM conditions. They take a cautious approach by discretizing fluxes on the dual mesh. The Günter scheme, originally applied to magnetized plasma simulations, reduces the spurious vertical diffusion observed in standard finite difference simulations when the anisotropy becomes large due to its discretization of the parallel operator [START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF]. Thanks to Günter's model, several researchers have studied the diffusion problem in magnetized plasma. For example, Günter [START_REF] Günter | A mixed implicit-explicit finite difference scheme for heat transport in magnetised plasmas[END_REF] applies her scheme and explicit/implicit domain decomposition techniques to the problem of time-dependent heat conduction in a two-dimensional strongly anisotropic medium (magnetized plasma), using a spatial derivative formulation that avoids contamination of the vertical direction by parallel heat flow. Holzl et al. [START_REF] Hölzl | Numerical modeling of diffusive heat transport across magnetic islands and highly stochastic layers[END_REF] performed numerical simulations of diffusive heat transfer between magnetic islands and extremely random layers. Van es et al. [START_REF] Van Es | Finite-difference schemes for anisotropic diffusion[END_REF] developed several discretization schemes and applied them to the anisotropic heat diffusion equation on a non-aligned grid. Soler et al. [START_REF] Soler | A new conservative finite-difference scheme for anisotropic elliptic problems in bounded domain[END_REF] proposed a new conservative finite-difference scheme using non-aligned Cartesian grids and interpolations aligned along a parallel diffusion direction for anisotropic elliptic problems in bounded domain. In addition to the above-mentioned researchers, some new approaches have been proposed for anisotropic diffusion. For example, Yang et al. [START_REF] Yang | Two uniform tailored finite point method for strongly anisotropic and discontinuous diffusivity[END_REF] proposed two Tailored Finite Point Methods (TFPM) for the diffusion equation, which are valid for strongly anisotropic tensor diffusion and interfacial layers. The five-point scheme uses the values on the grid points and their derivatives. The second scheme constructs a four-point scheme for each cell based on the interface conditions. And later, Yang et al. [START_REF] Yang | A new tailored finite point method for strongly anisotropic diffusion equation on misaligned grids[END_REF] applied this method to on misaligned grids. Chamarthi et al. [START_REF] Chamarthi | First order hyperbolic approach for anisotropic diffusion equation[END_REF] presented a high order finite difference solver for anisotropic diffusion problems based on the first-order hyperbolic system method. Feng et al. [START_REF] Feng | Sixth order compact finite difference schemes for poisson interface problems with singular sources[END_REF] proposed a sixth order compact finite difference scheme on uniform Cartesian grids for the Poisson interface problem with singular sources.

2 6 6 6 6 4
q 1 q 2 q 3 3 7 7 7 7 5 = - 2
The Günter symmetry scheme will be used as the core method in this work due to its powerful ability to solve anisotropic diffusion problems. At the same time this method discretizes the fluxes on a dual grid located at the center of the cell, and this coincides with our method of determining the thermal conductivity of the voxelized grid.

2D scheme

There are two schemes that are presented in Günter's article [START_REF] Günter | Modelling of heat transport in magnetised plasmas using non-aligned coordinates[END_REF], the asymmetric scheme and the symmetric scheme. The former defines the solution point on the boundary and the latter defines the solution point at the centre of the cell. Since we have chosen voxel discretization, where a voxel represents only one material when defining the material, the asymmetric model is not applicable to this work due to the inability to identify the material attribution of the points on the boundary. So the symmetric model will be used. Here is an attempt to re-demonstrate the Günter scheme formulation using interpolation. The first finite-difference scheme is shown in Fig. 2.1 In an isotropic case, first of all, we use the interpolation method through the points (i, j), (i+1, j),(i+1, j+1), et (i, j+1).

Symmetric scheme of Günter

T (x; y ) = a 1 + a 2 x + a 3 y + a 4 xy ; x ∈ [0; ∆x]; y ∈ [0; ∆y ] (2.15) 0 B B B B B B B @ 1 0 0 0 1 ∆x 0 0 1 ∆x ∆y ∆x∆y 1 0 ∆y 0 1 C C C C C C C A • 0 B B B B B B B @ a1 a2 a3 a4 1 C C C C C C C A = 0 B B B B B B B @ T i;j T i+1;j T i+1;j+1 T i;j+1 1 C C C C C C C A (2.16)
Then, we can get the coefficient a 1 ,...a 4 : Finally we get the same gradient formula as Günter at the central point of cell:

a 1 = T i;j ; a 2 = T i+1;j -T i;j ∆x ; a 3 = T i;j+1
@T @x | i+ 1 2 ;j+ 1 2 = T i+1;j+1 + T i+1;j -T i;j+1 -T i;j 2∆x ; @T @y | i+ 1 2 ;j+ 1 2 = T i;j+1 + T i+1;j+1 -T i+1;j -T i;j 2∆y : (2.19) i,j+1 i+1,j+1 i+1,j D i+1/2,j+1/2 D i+1/2,j-1/2 q=-D• Δ T Δ •q Figure 2.2: Symmetric scheme point of cell. - → q = -D • ∇T; q i+ 1 2 ;j+ 1 2 = -D i+ 1 2 ;j+ 1 2 • " @T @x | i+ 1 2 ;j+ 1 2 ; @T @y | i+ 1 2 ;j+ 1 2 « T (2.20)
Finally, with the heat flux, the formulation of the thermal conduction operator is:

∇ • - → q = q x;i+ 1 2 ;j+ 1 2 + q x;i+ 1 2 ;j-1 2 -q x;i-1 2 ;j+ 1 2 -q x;i-1 2 ;j-1 2 2∆x +
q y ;i+ 1 2 ;j+ 1 2 + q y ;i-1 2 ;j+ 1 2 -q y ;i+ 1 2 ;j-1 2 -q y ;i-1 2 ;j-

(D i-1=2;j-1=2 + D i-1=2;j+1=2 ) (0:25=∆x 2 + 0:25=∆y 2 ) (-D i-1=2;j-1=2 -D i-1=2;j+1=2 -D i+1=2;j-1=2 -D i+1=2;j+1=2 ) (0:25=∆x 2 -0:25=∆y 2 ) (D i+1=2;j-1=2 + D i+1=2;j+1=2 )) j-1 (0:25=∆x 2 + 0:25=∆y 2 ) D i-1=2;j-1=2 (-0:25=∆x 2 + 0:25=∆y 2 ) (D i-1=2;j-1=2 + D i+1=2;j-1=2 ) (0:25=∆x 2 + 0:25=∆y 2 ) D i+1=2;j-1=2
In the case of anisotropic solids, first of all, the diffusion tensor is:

D(x) = 2 6 4 D 11 (x) D 12 (x) D 21 (x) D 22 (x) 3 7 5 (2.23)
Then, with the new diffusion tensor, the heat flow is:

- → q = -D • ∇T; 2 6 4 q x;i+ 1 2 ;j+ 1 2 q y ;i+ 1 2 ;j+ 1 2 3 7 5 = - 2 6 4 D 11 (x) D 12 (x) D 21 (x) D 22 (x) 3 7 5 • 2 6 4 @T @x | i+ 1 2 ;j+ 1 2 @T @y | i+ 1 2 ;j+ 1 2 3 7 5 (2.24)
Finally, the formulation of finite difference in thermal conduction is :

∇ • - → q = -( `∆x 2 T i+1;j+1 -∆x 2 T i+1;j +∆x 2 T i;j+1 -∆x 2 T i;j ´D22 i+1=2;j+1=2
+ (∆x∆y T i+1;j+1 +∆x∆y T i+1;j -∆x∆y T i;j+1 -∆x∆y T i;j ) D21 i+1=2;j+1=2 + (∆x∆y T i+1;j+1 -∆x∆y T i+1;j +∆x∆y T i;j+1 -∆x∆y T i;j ) D12 i+1=2;j+1=2 + `∆y 2 T i+1;j+1 +∆y 2 T i+1;j -∆y 2 T i;j+1 -∆y 2 T i;j ´D11 i+1=2;j+1=2

+ `-∆x 2 T i+1;j +∆x 2 T i+1;j-1 -∆x 2 T i;j +∆x 2 T i;j-1 ´D22 i+1=2;j-1=2
+ (-∆x∆y T i+1;j -∆x∆y T i+1;j-1 +∆x∆y T i;j +∆x∆y T i;j-1 ) D21 i+1=2;j-1=2 + (∆x∆y T i+1;j -∆x∆y T i+1;j-1 +∆x∆y T i;j -∆x∆y T i;j-1 ) D12 i+1=2;j-1=2 + `∆y 2 T i+1;j +∆y 2 T i+1;j-1 -∆y 2 T i;j -∆y 2 T i;j-1 ´D11 i+1=2;j-1=2

+ `∆x 2 T i;j+1 -∆x 2 T i;j +∆x 2 T i-1;j+1 -∆x 2 T i-1;j ´D22 i-1=2;j+1=2
+ (∆x∆y T i;j+1 +∆x∆y T i;j -∆x∆y T i-1;j+1 -∆x∆y T i-1;j ) D21 i-1=2;j+1=2 + (-∆x∆y T i;j+1 +∆x∆y T i;j -∆x∆y T i-1;j+1 +∆x∆y T i-1;j ) D12 i-1=2;j+1=2 + `-∆y 2 T i;j+1 -∆y 2 T i;j +∆y 2 T i-1;j+1 +∆y 2 T i-1;j ´D11 i-1=2;j+1=2 + `-∆x 2 T i;j +∆x 2 T i;j-1 -∆x 2 T i-1;j +∆x 2 T i-1;j-1 ´D22 i-1=2;j-1=2 + (-∆x∆y T i;j -∆x∆y T i;j-1 +∆x∆y T i-1;j +∆x∆y T i-1;j-1 ) D21 i-1=2;j-1=2 + (-∆x∆y T i;j +∆x∆y T i;j-1 -∆x∆y T i-1;j +∆x∆y T i-1;j-1 ) D12 i-1=2;j-1=2 + `-∆y 2 T i;j -∆y 2 T i;j-1 +∆y 2 T i-1;j +∆y 2 T i-1;j-1 ´D11 i-1=2;j-1=2 )=(4∆x 2 ∆y 2 )

(2.25)

By simplifying we get the final relation: Table 2.5: 2D coefficient for T i;j with anisotropic conductivity Due to the heterogeneous material, the thermal conduction does not vary linearly so the interpolation method is not applicable in this case. In this regard, the 9-points integral method is used to solve for the temperature gradient at the grid points.

j Coef i i -1 i i + 1 j+1 -D11 i-1=2;j+1=2 =(4∆x 2 ) -D22 i-1=2;j+1=2 =(4∆y 2 ) +(D12 i-1=2;j+1=2 + D21 i-1=2;j+1=2 )=(4∆x∆y ) (D11 i-1=2;j+1=2 + D11 i+1=2;j+1=2 ) =(4∆x 2 )- (D22 i-1=2;j+1=2 + D22 i+1=2;j+1=2 ) =(4∆y 2 )+ (D21 i+1=2;j+1=2 -D12 i+1=2;j+1=2 ) -D21 i-1=2;j+1=2 + D12 i-1=2;j+1=2 ) =(4∆x∆y ) -D11 i+1=2;j+1=2 =(4∆x 2 ) -D22 i+1=2;j+1=2 =(4∆y 2 ) +(D12 i+1=2;j+1=2 + D21 i+1=2;j+1=2 )=(4∆x∆y ) j -(D11 i-1=2;j-1=2 + D11 i-1=2;j+1=2 ) =(4∆x 2 )+ (D22 i-1=2;j-1=2 + D22 i-1=2;j+1=2 ) =(4∆y 2 )+ (D21 i-1=2;j+1=2 + D12 i-1=2;j-1=2 ) -D21 i-1=2;j-1=2 -D12 i-1=2;j+1=2 ) =(4∆x∆y ) (D11 i+1=2;j+1=2 + D11 i+1=2;j-1=2 + D11 i-1=2;j-1=2 + D11 i-1=2;j+1=2 ) =(4∆x 2 )+ (D22 i+1=2;j+1=2 + D22 i+1=2;j-1=2 + D22 i-1=2;j-1=2 + D22 i-1=2;j+1=2 ) =(4∆y 2 )+ (D21 i+1=2;j+1=2 + D12 i+1=2;j+1=2 -D21 i+1=2;j-1=2 -D12 i+1=2;j-1=2 -D21 i-1=2;j+1=2 -D12 i-1=2;j+1=2 + D21 i-1=2;j-1=2 + D12 i+1=2;j+1=2 =(4∆x∆y ) -(D11 i+1=2;j-1=2 + D11 i+1=2;j-1=2 ) =(4∆x 2 )+ (D22 i+1=2;j+1=2 + D22 i+1=2;j-1=2 ) =(4∆y 2 )+ (D21 i+1=2;j+1=2 + D12 i+1=2;j-1=2 ) -D21 i+1=2;j-1=2 -D12 i+1=2;j+1=2 ) =(4∆x∆y ) j-1 -D11 i-1=2;j-1=2 =(4∆x 2 ) -D22 i-1=2;j-1=2 =(4∆y 2 ) +(D12 i-1=2;j-1=2 + D21 i-1=2;j-1=2 )=(4∆x∆y ) (D11 i+1=2;j-1=2 + D11 i-1=2;j-1=2 ) =(4∆x 2 )- (D22 i+1=2;j-1=2 + D22 i-1=2;j-1=2 ) =(4∆y 2 )- (D21 i+1=2;j-1=2 -D12 i+1=2;j-1=2 ) -D21 i-1=2;j-1=2 + D12 i-1=2;j-1=2 ) =(4∆x∆y ) -D11 i+1=2;j-1=2 =(4∆x 2 ) -D22 i+1=2;j-1=2 =(4∆y 2 ) +(D12 i+1=2;j-1=2 + D21 i+1=2;j-1=2 )=(4∆x∆y ) Gradient operator i,j+1 i,j i+1,j i+1,j+1 i+1/2,j+1/2 Δx Δy i-1,j+1 i-1,j i-1,j-1 i,j-1 i+1,j-1 A1 A2 A3 A4 i-1/2,j+1/2 i-1/2,j-1/2 i+1/2,j-1/2
Although the temperature gradient is not continuous for the whole region, it is still continuous within each pixel (voxel). So we can still use interpolation to derive the equation for the temperature gradient within each of the four pixels (voxels) from the temperature values on the grid.

With Fig. 2.3, we can obtain the formula for the temperature gradient within each pixel.

@T @x ˛A1 (x; y ) = T i;j -T i-1;j ∆x + T i;j -T i;j-1 -T i-1;j + T i-1;j-1 ∆x∆y y (2.26) @T @y ˛A1 (x; y ) = T i;j -T i;j-1 ∆y + T i;j -T i;j-1 -T i-1;j + T i-1;j-1 ∆x∆y x (2.27) @T @x ˛A2 (x; y ) = T i;j -T i-1;j ∆x + T i;j+1 -T i;j -T i-1;j+1 + T i-1;j ∆x∆y y (2.28) @T @y ˛A2 (x; y ) = T i;j+1 -T i;j ∆y + T i;j+1 -T i;j -T i-1;j+1 + T i-1;j ∆x∆y x (2.29) @T @x ˛A3 (x; y ) = T i+1;j -T i;j ∆x + T i+1;j -T i+1;j-1 -T i;j + T i;j-1 ∆x∆y y (2.30) @T @y ˛A3 (x; y ) = T i;j -T i;j-1 ∆y + T i+1;j -T i+1;j-1 -T i;j + T i;j-1 ∆x∆y x (2.31)
@T @x ˛A4 (x; y ) = T i+1;j -T i;j ∆x + T i+1;j+1 -T i+1;j -T i;j+1 + T i;j ∆x∆y y (2.32) @T @y ˛A4 (x; y ) = T i;j+1 -T i;j ∆y + T i+1;j+1 -T i+1;j -T i;j+1 + T i;j ∆x∆y x (2.33)

Through integral, we can get:

1 |A total | Z A total @T @x dA total = 1 |A total | `ZA 1 @T @x ˛A1 (x; y )dA 1 + Z A 2 @T @x ˛A2 (x; y )dA 2 + Z A 3 @T @x ˛A3 (x; y )dA 3 + Z A 4 @T @x ˛A4 (x; y )dA 4 ´(2.34)
Finally, we can obtain: @T @x ˛i;j = @T @x ˛i-1=2;j-1=2 + @T @x ˛i-1=2;j+1=2 + @T @x ˛i+1=2;j-1=2 + @T @x ˛i+1=2;j+1=2 4

(2.35) @T @y ˛i;j = @T @y ˛i-1=2;j-1=2 + @T @y ˛i-1=2;j+1=2 + @T @y ˛i+1=2;j-1=2 + @T @y ˛i+1=2;j+1=2 4

(2.36) 51

In the same way, we can also get the descriptions of approximation for heat flux as:

q x | i;j = q x ˛i-1=2;j-1=2 + q x ˛i-1=2;j+1=2 + q x ˛i+1=2;j-1=2 + q x ˛i+1=2;j+1=2 4
(2.37)

q y | i;j =
q y ˛i-1=2;j-1=2 + q y ˛i-1=2;j+1=2 + q y ˛i+1=2;j-1=2 + q y ˛i+1=2;j+1=2

4

(2.38)

3D scheme

After re-demonstrating the 2D model, we can use the same method to implement the 3D extension. Take the 1/8 model as an example. First, a cartesian coordinate system { -→ e x ; -→ e y ; -→ e z } is defined in Fig. 2.4. The subscripts {i; j; k} denote the spatial discretized grid points along the directions {x; y ; z}, respectively. The gradient of the temperature ∆T at the center point (x; y ; z) = (i + 1=2; j -1=2; k + 1=2) can be expressed as: @T @x ˛i+ 1 2 ;j-1 2 ;k+

1 2 = 1 4∆x
• (T i+1;j-1;k+1 + T i+1;j-1;k + T i+1;j;k+1 + T i+1;j;k -T i;j-1;k+1 -T i;j-1;k -T i;j;k+1 -T i;j;k ) @T @y ˛i+ 1 2 ;j-1 2 ;k+

1 2 = 1 4∆y
• (T i+1;j;k+1 + T i+1;j;k + T i;j;k+1 + T i;j;k -T i+1;j-1;k+1 -T i+1;j-1;k -T i;j-1;k+1 -T i;j-1;k ) @T @z ˛i+ 1 2 ;j-1 2 ;k+

1 2 = 1 4∆z
• (T i+1;j;k+1 + T i+1;j-1;k+1 + T i;j;k+1 + T i;j-1;k+1

-T i+1;j;k -T i+1;j-1;k -T i;j;k -T i;j-1;k ) (2.39)
By inserting these terms into the first formula of Eq.2.13 and applying the diffusion tensor D, we obtain the heat flux -→ q at the center point:

-→ q i+ 1 2 ;j+ 1 2 ;k+ 1 2 = -D i+ 1 2 ;j+ 1 2 ;k+ 1 2

• " @T @x ˛i+ 1 2 ;j+ 1 2 ;k+ 1 2 ; @T @y ˛i+ 1 2 ;j+ 1 2 ;k+ 1 2 ; @T @z ˛i+ 1 2 ;j+ 1 2 ;k+

1 2 « T (2.40)
To take the divergence over the heat flux, we have the 3D formulation of the thermal conduction problem:

∇ • - → q = 1 4∆x
• (q x;i+ 1 2 ;j+ 1 2 ;k+ 1 2 + q x;i+ 1 2 ;j+ 1 2 ;k-1 2 + q x;i+ 1 2 ;j-1 2 ;k+ 1 2 + q x;i+ 1 2 ;j-1 2 ;k--q x;i-1 2 ;j+ 1 2 ;k+ 1 2 -q x;i-1 2 ;j-1 2 ;k+

1 2 -q x;i-1 2 ;j+ 1 2 ;k-1 2 -q x;i-1 2 ;j- ;k-1 2 )+ 1 4∆y
• (q y ;i+ 1 2 ;j+ 1 2 ;k+ 1 2 + q y ;i-1 2 ;j+ 1 2 ;k+ 1 2 + q y ;i+ 1 2 ;j+ 1 2 ;k-1 2 + q y ;i-1 2 ;j+ 1 2 ;k-1 2 -q y ;i+ 1 2 ;j-1 2 ;k+ 1 2 -q y ;i-1 2 ;j-1 2 ;k+ 1 2 -q y ;i+ 1 2 ;j-1 2 ;k-

1 2 -q y ;i-1 2 ;j- ;k-1 2 )+ 1 4∆z
• (q z;i+ 1 2 ;j+ 1 2 ;k+ 1 2 + q z;i-1 2 ;j+ 1 2 ;k+ 1 2 + q z;i+ 1 2 ;j-1 2 ;k+ 1 2 + q z;i-1 2 ;j-1 2 ;k+ 1 2 -q z;i+ 1 2 ;j+ 1 2 ;k-

1 2 -q z;i+ 1 2 ;j-1 2 ;k-1 2 -q z;i-1 2 ;j+ 1 2 ;k-1 2 -q z;i-1 2 ;j-1 2 ;k-1 2 ) (2.41)
The expansion of Eq. 2.41 is given in Appendix A.

By simplifying we get the final relation in Table .2.6 for isotropic case.

Table 2.6: 3D Günter scheme's coefficient for T i;j with isotropic conductivity

i-1 i i+1 k+1 j+1 -D i-1 2 ;j+ 1 2 ;k+ 1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) (-D i+ 1 2 ;j+ 1 2 ;k+ 1 2 -D i-1 2 ;j+ 1 2 ;k+ 1 2 ) (1/(16∆z 2 )+1/(16∆y 2 ) -1/(16∆x 2 )) -D i+ 1 2 ;j+ 1 2 ;k+ 1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) k+1 j (-D i-1 ) (-1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) k-1 j+1 -D i-1 2 ;j+ 1 2 ;k-1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) (-D i+ 1 2 ;j+ 1 2 ;k-1 2 -D i-1 2 ;j+ 1 2 ;k-1 2 ) (1/(16∆z 2 )+1/(16∆y 2 ) -1/(16∆x 2 )) -D i+ 1 2 ;j+ 1 2 ;k-1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) k-1 j (-D i-1 2 ;j+ 1 2 ;k-1 2 -D i-1 2 ;j-1 2 ;k-1 2 ) (1/(16∆z 2 )-1/(16∆y 2 ) +1/(16∆x 2 )) (-D i+ 1 2 ;j+ 1 2 ;k-1 2 -D i+ 1 2 ;j-1 2 ;k-1 2 -D i-1 2 ;j+ 1 2 ;k-1 2 -D i-1 2 ;j-1 2 ;k-1 2 ) (1/(16∆z 2 )-1/(16∆y 2 ) -1/(16∆x 2 )) (-D i+ 1 2 ;j+ 1 2 ;k-1 2 -D i+ 1 2 ;j-1 2 ;k-1 2 ) (1/(16∆z 2 )-1/(16∆y 2 ) +1/(16∆x 2 )) k-1 j-1 -D i-1 2 ;j-1 2 ;k-1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) (-D i+ 1 2 ;j-1 2 ;k-1 2 -D i-1 2 ;j-1 2 ;k-1 2 ) (1/(16∆z 2 )+1/(16∆y 2 ) -1/(16∆x 2 )) -D i+ 1 2 ;j-1 2 ;k-1 2 (1/(16∆z 2 )+1/(16∆y 2 ) +1/(16∆x 2 )) 54
In the case of anisotropic solids 3D, the diffusion tensor is: (2.42)

D(x) = 2
Then, with the new diffusion tensor, the heat flow is:

- → q = -D • ∇T; 2 6 6 6 6 4
q x;i+ 1 2 ;j+ 1 2 ;k+ 1 2 q y ;i+ 1 2 ;j+ 1 2 ;k+ 1 2 q z;i+ 1 2 ;j+ 1 2 ;k+

+ D11 i-1 2 ;j+ 1 2 ;k+ 1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j+ 1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j+ 1 2 ;k+ 1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) +(-D12 i+ 1 2 ;j+ 1 2 ;k+ 1 2 +D21 i+ 1 2 ;j+ 1 2 ;k+ 1 2 -D12 i+ 1 2 ;j-1 2 ;k+ 1 2 -D21 i+ 1 2 ;j-1 2 ;k+ 1 2 -D12 i-1 2 ;j+ 1 2 ;k+ 1 2 -D21 i-1 2 ;j+ 1 2 ;k+ 1 2 +D12 i-1 2 ;j-1 2 ;k+ 1 2 +D21 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆x∆y ) +(-D13 i+ 1 2 ;j+ 1 2 ;k+ 1 2 +D31 i+ 1 2 ;j+ 1 2 ;k+ 1 2 -D13 i+ 1 2 ;j-1 2 ;k+ 1 2 +D31 i+ 1 2 ;j-1 2 ;k+ 1 2 +D13 i-1 2 ;j+ 1 2 ;k+ 1 2 -D31 i-1 2 ;j+ 1 2 ;k+ 1 2 +D13 i-1 2 ;j-1 2 ;k+ 1 2 -D31 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆x∆z) +(-D23 i+ 1 2 ;j+ 1 2 ;k+ 1 2 -D32 i+ 1 2 ;j+ 1 2 ;k+ 1 2 +D23 i+ 1 2 ;j-1 2 ;k+ 1 2 -D32 i+ 1 2 ;j-1 2 ;k+ 1 2 -D23 i-1 2 ;j+ 1 2 ;k+ 1 2 +D32 i-1 2 ;j+ 1 2 ;k+ 1 2 +D23 i-1 2 ;j-1 2 ;k+ 1 2 -D32 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆y ∆z) -D11 i+ ;j+ 1 2 ;k+ 1 2 / (16 ∆x 2 ) + D22 i+ 1 2 ;j+ 1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i+ ;j+ 1 2 ;k+ 1 2 / (16 ∆z 2 ) -D11 i+ ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) + D22 i+ 1 2 ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i+ ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) +(D12 i+ 1 2 ;j+ 1 2 ;k+ 1 2 -D21 i+ ;j+ 1 2 ;k+ 1 2 -D12 i+ ;j-1 2 ;k+ 1 2 +D21 i+ ;j-1 2 ;k+ 1 2 )/(16 ∆x∆y ) -(D13 i+ ;j+ 1 2 ;k+ 1 2 +D31 i+ ;j+ 1 2 ;k+ 1 2 D13 i+ 1 2 ;j-1 2 ;k+ 1 2 +D31 i+ ;j-1 2 ;k+ 1 2
)/(16 ∆x∆z)

+(-D23 i+ 1 2 ;j+ 1 2 ;k+ 1 2 +D32 i+ ;j+ 1 2 ;k+ 1 2 +D23 i+ ;j-1 2 ;k+ 1 2 -D32 i+ ;j-1 2 ;k+ 1 2
)/(16 ∆y ∆z)

m+1 j-1 -D11 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) -D22 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) -(D12 i-1 2 ;j-1 2 ;k+ 1 2 +D21 i-1 2 ;j-1 2 ;k+ 1 2
)/(16 ∆x∆y )

+(D13 i-1 2 ;j-1 2 ;k+ 1 2 +D31 i-1 2 ;j-1 2 ;k+ 1 2
)/(16 ∆x∆z)

+(D23 i-1 2 ;j-1 2 ;k+ 1 2 +D32 i-1 2 ;j-1 2 ;k+ 1 2
)/(16 ∆y ∆z)

+ D11 i+ 1 2 ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) -D22 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) +(+D12 i+ 1 2 ;j-1 2 ;k+ 1 2 -D21 i+ 1 2 ;j-1 2 ;k+ 1 2 -D12 i-1 2 ;j-1 2 ;k+ 1 2 +D21 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆x∆y ) +(-D13 i+ 1 2 ;j-1 2 ;k+ 1 2 +D31 i+ 1 2 ;j-1 2 ;k+ 1 2 +D13 i-1 2 ;j-1 2 ;k+ 1 2 -D31 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆x∆z) +(+D23 i+ 1 2 ;j-1 2 ;k+ 1 2 +D32 i+ 1 2 ;j-1 2 ;k+ 1 2 +D23 i-1 2 ;j-1 2 ;k+ 1 2 +D32 i-1 2 ;j-1 2 ;k+ 1 2 ) /(16 ∆y ∆z) -D11 i+ ;j-1 2 ;k+ 1 2 / (16 ∆x 2 ) -D22 i+ ;j-1 2 ;k+ 1 2 / (16 ∆y 2 ) -D33 i+ ;j-1 2 ;k+ 1 2 / (16 ∆z 2 ) +(D12 i+ 1 2 ;j-1 2 ;k+ 1 2 +D21 i+ ;j-1 2 ;k+ 1 2 )/(16 ∆x∆y ) -(D13 i+ ;j-1 2 ;k+ 1 2 +D31 i+ ;j-1 2 ;k+ 1 2
)/(16 ∆x∆z)

+(D23 i+ 1 2 ;j-1 2 ;k+ 1 2 +D32 i+ ;j-1 2 ;k+ 1 2
)/(16 ∆y ∆z)

) / (16 ∆x ∆z ) +(D23 i + 1 2 ;j + 1 2 ;k + 1 2 -D23 i + 1 2 ;j + 1 2 ;k -1 2 -D23 i + 1 2 ;j -1 2 ;k + 1 2 +D23 i + 1 2 ;j -1 2 ;k -1 2 +D32 i + 1 2 ;j + 1 2 ;k + 1 2 -D32 i + 1 2 ;j + 1 2 ;k -1 2 -D32 i + 1 2 ;j -1 2 ;k + 1 2 +D32 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y ∆z ) m j-1 (-D11 i -1 2 ;j -1 2 ;k + 1 2 -D11 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x 2 ) +(-D22 i -1 2 ;j -1 2 ;k + 1 2 -D22 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y 2 ) +(D33 i -1 2 ;j -1 2 ;k + 1 2 +D33 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆z 2 ) +(-D12 i -1 2 ;j -1 2 ;k + 1 2 -D12 i -1 2 ;j -1 2 ;k -1 2 -D21 i + 1 2 ;j -1 2 ;k + 1 2 -D21 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆y ) +(-D13 i -1 2 ;j -1 2 ;k + 1 2 +D13 i -1 2 ;j -1 2 ;k -1 2 +D31 i -1 2 ;j -1 2 ;k + 1 2 -D31 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆z ) +(-D23 i -1 2 ;j -1 2 ;k + 1 2 +D23 i -1 2 ;j -1 2 ;k -1 2 +D32 i -1 2 ;j -1 2 ;k + 1 2 -D32 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y ∆z ) (+D11 i + 1 2 ;j -1 2 ;k + 1 2 +D11 i + 1 2 ;j -1 2 ;k -1 2 +D11 i -1 2 ;j -1 2 ;k + 1 2 +D11 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x 2 ) +(-D22 i + 1 2 ;j -1 2 ;k + 1 2 -D22 i + 1 2 ;j -1 2 ;k -1 2 -D22 i -1 2 ;j -1 2 ;k + 1 2 -D22 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y 2 ) +(D33 i + 1 2 ;j -1 2 ;k + 1 2 +D33 i + 1 2 ;j -1 2 ;k -1 2 +D33 i -1 2 ;j -1 2 ;k + 1 2 +D33 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆z 2 ) +(D12 i + 1 2 ;j -1 2 ;k + 1 2 +D12 i + 1 2 ;j -1 2 ;k -1 2 -D12 i -1 2 ;j -1 2 ;k + 1 2 -D12 i -1 2 ;j -1 2 ;k -1 2 -D21 i + 1 2 ;j -1 2 ;k + 1 2 -D21 i + 1 2 ;j -1 2 ;k -1 2 +D21 i + 1 2 ;j -1 2 ;k + 1 2 +D21 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆y ) +(D13 i + 1 2 ;j -1 2 ;k + 1 2 -D13 i + 1 2 ;j -1 2 ;k -1 2 -D13 i -1 2 ;j -1 2 ;k + 1 2 +D13 i -1 2 ;j -1 2 ;k -1 2 +D31 i + 1 2 ;j -1 2 ;k + 1 2 -D31 i + 1 2 ;j -1 2 ;k -1 2 -D31 i -1 2 ;j -1 2 ;k + 1 2 +D31 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆z ) +(-D23 i + 1 2 ;j -1 2 ;k + 1 2 +D23 i + 1 2 ;j -1 2 ;k -1 2 -D23 i -1 2 ;j -1 2 ;k + 1 2 +D23 i -1 2 ;j -1 2 ;k -1 2 +D32 i + 1 2 ;j -1 2 ;k + 1 2 -D32 i + 1 2 ;j -1 2 ;k -1 2 +D32 i -1 2 ;j -1 2 ;k + 1 2 -D32 i -1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y ∆z ) (-D11 i + 1 2 ;j -1 2 ;k + 1 2 -D11 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x 2 ) +(-D22 i + 1 2 ;j -1 2 ;k + 1 2 -D22 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y 2 ) +(D33 i + 1 2 ;j -1 2 ;k + 1 2 +D33 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆z 2 ) +(D12 i + 1 2 ;j -1 2 ;k + 1 2 +D12 i + 1 2 ;j -1 2 ;k -1 2 +D21 i + 1 2 ;j -1 2 ;k + 1 2 +D21 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆y ) +(+D13 i + 1 2 ;j -1 2 ;k + 1 2 -D13 i + 1 2 ;j -1 2 ;k -1 2 -D31 i + 1 2 ;j -1 2 ;k + 1 2 +D31 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆x ∆z ) +(-D23 i + 1 2 ;j -1 2 ;k + 1 2 +D23 i + 1 2 ;j -1 2 ;k -1 2 +D32 i + 1 2 ;j -1 2 ;k + 1 2 -D32 i + 1 2 ;j -1 2 ;k -1 2 ) / (16 ∆y ∆z ) i-1 i i+1 m-1 j+1 -D11 i- ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i- ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i- ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) +(D12 i-1 2 ;j+ 1 2 ;k-1 2 +D21 i- ;j+ 1 2 ;k-1 2 )/(16 ∆x∆y ) -(D13 i- ;j+ 1 2 ;k-1 2 +D31 i- ;j+ 1 2 ;k-1 2
)/(16 ∆x∆z)

+(D23 i-1 2 ;j+ 1 2 ;k-1 2 +D32 i- ;j+ 1 2 ;k-1 2
)/(16 ∆y ∆z)

+ D11 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) +(-D12 i+ 1 2 ;j+ 1 2 ;k-1 2 +D21 i+ 1 2 ;j+ 1 2 ;k-1 2 +D12 i-1 2 ;j+ 1 2 ;k-1 2 -D21 i-1 2 ;j+ 1 2 ;k-1 2 ) /(16 ∆x∆y ) -(-D13 i+ 1 2 ;j+ 1 2 ;k-1 2 +D31 i+ 1 2 ;j+ 1 2 ;k-1 2 +D13 i-1 2 ;j+ 1 2 ;k-1 2 -D31 i-1 2 ;j+ 1 2 ;k-1 2 ) /(16 ∆x∆z) -(-D23 i+ 1 2 ;j+ 1 2 ;k-1 2 -D32 i+ 1 2 ;j+ 1 2 ;k-1 2 -D23 i-1 2 ;j+ 1 2 ;k-1 2 -D32 i-1 2 ;j+ 1 2 ;k-1 2
) /(16 ∆y ∆z)

-D11 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) -(D12 i+ 1 2 ;j+ 1 2 ;k-1 2 +D21 i+ 1 2 ;j+ 1 2 ;k-1 2
)/(16 ∆x∆y )

+(D13 i+ 1 2 ;j+ 1 2 ;k-1 2 +D31 i+ 1 2 ;j+ 1 2 ;k-1 2
)/(16 ∆x∆z)

+(D23 i+ 1 2 ;j+ 1 2 ;k-1 2 +D32 i+ 1 2 ;j+ 1 2 ;k-1 2
)/(16 ∆y ∆z)

m-1 j -D11 i- ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i- ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) -D11 i- ;j-1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i- ;j-1 2 ;k-1 2 / (16 ∆z 2 ) +(-D12 i-1 2 ;j+ 1 2 ;k-1 2 +D21 i- ;j+ 1 2 ;k-1 2 +D12 i- ;j-1 2 ;k-1 2 -D21 i- ;j-1 2 ;k-1 2
)/(16 ∆x∆y )

-(D13 i- ;j+ 1 2 ;k-1 2 +D31 i- ;j+ 1 2 ;k-1 2 +D13 i- ;j-1 2 ;k-1 2 +D31 i- ;j-1 2 ;k-1 2
)/(16 ∆x∆z)

-(-D23 i-1 2 ;j+ 1 2 ;k-1 2 +D32 i- ;j+ 1 2 ;k-1 2 +D23 i- ;j-1 2 ;k-1 2 -D32 i- ;j-1 2 ;k-1 2
)/(16 ∆y ∆z)

+ D11 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) + D11 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) +(-D12 i+ 1 2 ;j+ 1 2 ;k-1 2 +D21 i+ 1 2 ;j+ 1 2 ;k-1 2 -D12 i+ 1 2 ;j-1 2 ;k-1 2 -D21 i+ 1 2 ;j-1 2 ;k-1 2 -D12 i-1 2 ;j+ 1 2 ;k-1 2 -D21 i-1 2 ;j+ 1 2 ;k-1 2 +D12 i-1 2 ;j-1 2 ;k-1 2 +D21 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆x∆y ) -(-D13 i+ 1 2 ;j+ 1 2 ;k-1 2 +D31 i+ 1 2 ;j+ 1 2 ;k-1 2 -D13 i+ 1 2 ;j-1 2 ;k-1 2 +D31 i+ 1 2 ;j-1 2 ;k-1 2 +D13 i-1 2 ;j+ 1 2 ;k-1 2 -D31 i-1 2 ;j+ 1 2 ;k-1 2 +D13 i-1 2 ;j-1 2 ;k-1 2 -D31 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆x∆z) -(-D23 i+ 1 2 ;j+ 1 2 ;k-1 2 -D32 i+ 1 2 ;j+ 1 2 ;k-1 2 +D23 i+ 1 2 ;j-1 2 ;k-1 2 -D32 i+ 1 2 ;j-1 2 ;k-1 2 -D23 i-1 2 ;j+ 1 2 ;k-1 2 +D32 i-1 2 ;j+ 1 2 ;k-1 2 +D23 i-1 2 ;j-1 2 ;k-1 2 -D32 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆y ∆z) -D11 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j+ 1 2 ;k-1 2 / (16 ∆z 2 ) -D11 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) + D22 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) +(D12 i+ 1 2 ;j+ 1 2 ;k-1 2 -D21 i+ 1 2 ;j+ 1 2 ;k-1 2 -D12 i+ 1 2 ;j-1 2 ;k-1 2 +D21 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆x∆y )

+(D13 i+ 1 2 ;j+ 1 2 ;k-1 2 +D31 i+ 1 2 ;j+ 1 2 ;k-1 2 D13 i+ 1 2 ;j-1 2 ;k-1 2 +D31 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆x∆z)

-(-D23 i+ 1 2 ;j+ 1 2 ;k-1 2 +D32 i+ 1 2 ;j+ 1 2 ;k-1 2 +D23 i+ 1 2 ;j-1 2 ;k-1 2 -D32 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆y ∆z)

m-1 j-1 -D11 i- ;j-1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i- ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i- ;j-1 2 ;k-1 2 / (16 ∆z 2 ) -(D12 i- ;j-1 2 ;k-1 2 +D21 i- ;j-1 2 ;k-1 2 )/(16 ∆x∆y ) -(D13 i- ;j-1 2 ;k-1 2 +D31 i- ;j-1 2 ;k-1 2
)/(16 ∆x∆z)

-(D23 i- ;j-1 2 ;k-1 2 +D32 i- ;j-1 2 ;k-1 2
)/(16 ∆y ∆z)

+ D11 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) + D11 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i-1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) +(+D12 i+ 1 2 ;j-1 2 ;k-1 2 -D21 i+ 1 2 ;j-1 2 ;k-1 2 -D12 i-1 2 ;j-1 2 ;k-1 2 +D21 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆x∆y ) -(-D13 i+ 1 2 ;j-1 2 ;k-1 2 +D31 i+ 1 2 ;j-1 2 ;k-1 2 +D13 i-1 2 ;j-1 2 ;k-1 2 -D31 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆x∆z) -(+D23 i+ 1 2 ;j-1 2 ;k-1 2 +D32 i+ 1 2 ;j-1 2 ;k-1 2 +D23 i-1 2 ;j-1 2 ;k-1 2 +D32 i-1 2 ;j-1 2 ;k-1 2 ) /(16 ∆y ∆z) -D11 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆x 2 ) -D22 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆y 2 ) -D33 i+ 1 2 ;j-1 2 ;k-1 2 / (16 ∆z 2 ) +(D12 i+ 1 2 ;j-1 2 ;k-1 2 +D21 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆x∆y )

+(D13 i+ 1 2 ;j-1 2 ;k-1 2 +D31 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆x∆z)

-(D23 i+ 1 2 ;j-1 2 ;k-1 2 +D32 i+ 1 2 ;j-1 2 ;k-1 2
)/(16 ∆y ∆z) 58

Homogenization

The main purpose of homogenization is to determine the <<equivalent>> properties of heterogeneous materials based on the knowledge of their components properties and spacial arrangement.

The classical homogenization technique aims at replacing the actual non-homogeneous composite with a fictitious homogeneous component for the purpose of calculating the effective properties of the composite component through the knowledge of its topological structure. To achieve homogenization we need to build two scales: macroscopic and microscopic scales with a "Representative Volume Element" (RVE) and an "equivalent homogeneous medium" (EHM), respectively [START_REF] Zaoui | Continuum micromechanics: survey[END_REF].

General principles Representative Volume Element

Since calculations for large volumes are often prohibitive the definition of the RVE is very important. Over the past decades many researchers have given definitions: Hill [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF] noted that the RVE must contain sufficient microstructural information and therefore include samples with heterogeneity such as voids, fibers, inclusions, grains, etc. According to Drugan and Willis [START_REF] Drugan | A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites[END_REF], the RVE is the smallest material volume element in the composite for which the conventional spatially constant "overall modulus" macroscopic constitutive representation is an accurate enough model to describe mean constitutive response. Kanit et al. [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] state that the size of the RVE must be considered in terms of five parameters: property contrast, physical properties, component volume fraction, relative accuracy required for effective property estimation, and the number of realizations of the microstructure relevant to the calculation. It is also pointed out that for all types of boundary conditions, when the RVE is too small, the estimation of the effective properties is biased.

Hill Lemma

Hill's lemma [START_REF] Hill | Elastic properties of reinforced solids: Some theoretical principles[END_REF], alternatively referred to as the macrohomogeneity or Hill-Mandel condition, says that the average of the product of the stress ff and strain " tensors in any

Representative Volume Element (RVE) (micro level) equals the product of their averages (macro level). It establishes the necessary and sufficient criteria for the equivalence of elastic materials' energy and mechanical properties.

< ff : " >=< ff >:< " > (2.44)
For thermal problem:

< q : ∇T >=< q >:< ∇T > (2.45)

Spatial average

The spatial averages of the local heat flux and temperature gradient are defined by:

< q >= 1 |Ω| Z Ω q(x)dV < ∇T >= 1 |Ω| Z Ω ∇T (x)dV (2.46)
Where |Ω| is the total volume of the RVE, q(x) and T (x) can be determined by applying a given temperature boundary condition; < q > and < ∇T > can be calculated by integrating the corresponding local fields over the domain Ω, then dividing it by the total volume.

With voxel(pixel) models, the integration operation is reduced to the multiplication of each q(x) and ∇T (x) by the volume fraction.

For a 2D case, if the point is in the corner, the volume fraction is ∆x∆y =(4

* (L x * L y )).
If the point is on edge, the volume fraction is ∆x∆y =(2 * (L x * L y )). If the point is inside, the volume fraction is ∆x∆y =(L x * L y ), where L x and L y represent the model length, respectively.

For a 3D case, if the point is in the corner, the volume fraction is ∆x∆y ∆z=(8 * (L x * L y * L z )).

If the point is on edge, the volume fraction is ∆x∆y ∆z=(4 Consequently, the macroscopic equivalent thermal conductivity D * is defined such that:

* (L x * L y * L z )). If
< q >= -D * • < ∇T > (2.47)
To numerically obtain each component of the tensor D * in Eq.2.47, the RVE needs to be simulated three times by applying boundary conditions in three orthogonal directions, respectively. Then, in the two-dimensional mode, we can obtain the following equation:

-

2 6 6 6 6 6 6 6 4 < ∇T 1 x > < ∇T 1 y > 0 0 0 0 < ∇T 1 x > < ∇T 1 y > < ∇T 2 x > < ∇T 2 y > 0 0 0 0 < ∇T 2 x > < ∇T 2 y > 3 7 7
< q 1 x > < q 1 y > < q 2 x > < q 2 y > (2.48)
and in the three-dimensional mode, we can obtain the following equation:

- (2.49)

2 6 6 6 6 6 6 6 6 4 <∇T 1 x > <∇T 1 y > <∇T 1 z > 0 0 0 0 0 0 0 0 0 <∇T 1 x > <∇T 1 y > <∇T 1 z > 0 0 0 0 0 0 0 0 0 <∇T 1 x > <∇T 1 y > <∇T 1 z > <∇T 2 x > <∇T 2 y > <∇T 2 z > 0 0 0 0 0 0 0 0 0 <∇T 2 x > <∇T 2 y > <∇T 2 z > 0 0 0 0 0 0 0 0 0 <∇T 2 x > <∇T 2 y > <∇T 2 z > <∇T 3 x > <∇T 3 y > <∇T 3 z > 0 0 0 0 0 0 0 0 0 <∇T 3 x > <∇T 3 y > <∇T 3 z > 0 0 0 0 0 0 0 0 0 <∇T 3 x > <∇T 3 y > <∇T
= < q 1 x > < q 1 y > < q 1 z > < q 2 x > < q 2 y > < q 2 z > < q 3 x > < q 3 y > < q 3 z > 3 7 7

Voigt-Reuss analytical estimates

Voigt and Reuss [START_REF] Saeb | Aspects of computational homogenization at finite deformations: A unifying review from reuss' to voigt's bound[END_REF][START_REF] Firooz | A systematic study on homogenization and the utility of circular simplified rve[END_REF] can be used to estimate the equivalent properties by supplying analytical upper and lower bound for thermal conductivity. For example, with a composite material which is consist of matrix and fiber inclusions.

Voigt can offer the upper bound.

D hom = D m V m + D f V f (2.50)
Reuss can provide the lower bound.

D hom -1 = D -1 m V m + D -1 f V f (2.51)
where D hom , D m , and D f represent the thermal conductivity of the homogenized material, matrix, and inclusion, respectively, V m and V f denotes the volume fraction of matrix, and inclusion, respectively.

Boundary conditions

In order to study the thermal behavior of inhomogeneous materials in great detail, we need to specify the appropriate boundary conditions for the solid elements. In the choice of the boundary conditions, standard boundary conditions include kinematically uniform (KUBC), statically uniform (SUBC), and periodicity (PBC). Uniform displacement-traction (orthogonal mixed) boundary conditions (MUBC) are less well-known. They were introduced by Hazanov and Amieur [START_REF] Hazanov | On overall properties of elastic heterogeneous bodies smaller than the representative volume[END_REF], Hazanov [START_REF] Hazanov | Hill condition and overall properties of composites[END_REF]. Nemat-Nasser and Hori [START_REF] Nemat-Nasser | Bounds and estimates of overall moduli of composites with periodic microstructure[END_REF]. discovered that SUBC results in lower estimations of their macroscopic stiffness, while KUBC results in greater values. Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] showed that PBC-based forecasts lie between these two extremes. Ostoja-Starzewski [START_REF] Ostoja-Starzewski | Material spatial randomness: from statistical to representative volume element[END_REF] showed that the predictions of MUBC boundary conditions are between those obtained with SUBC and KUBC. However, Hazanov [START_REF] Hazanov | Hill condition and overall properties of composites[END_REF] proposed that mixed homogeneous boundary conditions must be orthogonal and can be achieved only in materials having at least orthogonal elastic symmetry, with the orthogonal axes aligned with the hexahedral volume elements' symmetry axes. Pahr et Zysset [START_REF] Pahr | Influence of boundary conditions on computed apparent elastic properties of cancellous bone[END_REF] discovered that when a periodic and orthotropic microstructured material is used, the proposed mixed boundary conditions have the exact same effective elastic properties as periodic boundary conditions, and are thus referred to as "periodicity compatible" mixed uniform boundary conditions (PMUBCs). Thus, we use periodic boundary conditions for composites with anisotropic conductivity and periodic structure; for composites with nonperiodic structure and orthotropic conductivity, we use MUBC boundary conditions.

KUBC

Boundary conditions of uniform displacement in elasticity case (Dirichlet, kinematics, KUBC):

u i (x) = " 0 ij • x j (x); ∀x ∈ @Ω (2.52)
Boundary conditions of uniform gradient of temperature in thermal case.

T (x) = T g 0 i • x j (x); ∀x ∈ @Ω (2.53)

SUBC

Boundary conditions of uniform traction in elasticity case (Neumann, static, SUBC):

t i (x) = ff 0 ij • n j (x) ∀x ∈ @Ω (2.54)
Boundary conditions of uniform heat flux in in thermal case.

Q = q 0 i • n j (x) ∀x ∈ @Ω (2.55)

MUBC

Boundary conditions of uniform displacement-traction in elasticity case (mixed orthogonal, MUBC):

(u i (x) -" 0 ij • x j (x)) • (t i (x) -ff 0 ij • n j (x)) = 0; ∀x ∈ @Ω (2.56)
Boundary conditions of uniform gradient of temperature-heat flux in thermal case:

(T (x) -T g 0 i • x j (x)) • (Q i (x) -q 0 i • n j (x)) = 0; ∀x ∈ @Ω (2.57)
where " 0 ij , ff 0 ij , T g 0 i , and q 0 i are used to indicate constant tensors.

PBC

In the case of periodic microstructures, periodic boundary conditions (PBC) may be used.

For elasticity case:

u i (x) k+ -u i (x) k-= " 0 ij ∆x k ; ∀x ∈ @Ω (2.58)
For thermal case:

T (x) k+ -T (x) k-= T g 0 i ∆x k ; ∀x ∈ @Ω (2.59)
where k+ and kdenote a pair of parallel boundary surfaces, ∆x k denotes the constant distance between parallel planes.

Computational homogenization

PBC Periodic boundary conditions have been widely used in the field of homogenization, Tian [START_REF] Tian | Periodic boundary condition and its numerical implementation algorithm for the evaluation of effective mechanical properties of the composites with complicated micro-structures[END_REF][START_REF] Tian | New numerical algorithm for the periodic boundary condition for predicting the coefficients of thermal expansion of composites[END_REF][START_REF] Tian | Numerical evaluation on the effective thermal conductivity of the composites with discontinuous inclusions: Periodic boundary condition and its numerical algorithm[END_REF] details the method of adding periodic boundary conditions, summarized as follows:

When a material's microstructure is geometrically periodic, an RVE modeling approach exploits the periodicity by simulating the material as an endless repeat of the smallest periodic unit cell. For example in figure 2

u i (BCC B ) -u i (ADD A ) = ∆x" 0 ix (2.60) u i (ADCB) -u i (A D C B ) = ∆y " 0 iy (2.61) u i (ABB A ) -u i (DCC D ) = ∆z" 0 iz (2.62)
for edges:

u i (AA ) -u i (DD ) = ∆z" 0 iz (2.63) u i (CC ) -u i (DD ) = ∆x" 0 ix (2.64) u i (BB ) -u i (CC ) = ∆z" 0 iz (2.65) u i (AB) -u i (A B ) = ∆y " 0 iy (2.66) u i (A B ) -u i (D C ) = ∆z" 0 iz (2.67) u i (DC) -u i (D C ) = ∆y " 0 iy (2.68) u i (BC) -u i (AD) = ∆x" 0 ix (2.69) u i (AD) -u i (A D ) = ∆y " 0 iy (2.70) u i (B C ) -u i (A D ) = ∆x" 0 ix (2.71)
for vertex:

u i (A) -u i (D) = ∆z" 0 iz (2.72) u i (C) -u i (D) = ∆x" 0 ix (2.73) u i (B) -u i (C) = ∆z" 0 iz (2.74) u i (A ) -u i (D ) = ∆z" 0 iz (2.75) u i (C ) -u i (D ) = ∆x" 0 ix (2.76) u i (B ) -u i (C ) = ∆z" 0 iz (2.77) u i (D) -u i (D ) = ∆y " 0 iy (2.78)
In equations 2.61-2.78, ∆x, ∆y and ∆z denote the x, y, and z lengths of the RVE, respectively. The subscript i indicates a direction between the coordinates x, y, and z. " 0 ix , " 0 iy and " 0 iz are macroscopic strains used to specify particular boundary conditions, i.e., ∆x" 0 ix , ∆y " 0 iy and ∆z" 0 iz can be globally defined by the variation of the length or shear displacement of the RVE model. The PBC for temperature is similarly created. Equations 2.61-2.78 are modified by substituting the temperature change T for u i . The following are the precise implementation equations for the periodic boundary condition used to predict the equivalent thermal conductivity of composites:

Face-BCC'B' and Face-ADD'A':

T F ace-BCC B -T F ace-ADD A = ∆T x (2.79)
Face-ADCB and Face-A'D'C'B':

T F ace-ADCB -T F ace-A D C B = ∆T y (2.80)
Face-ABB'A' and Face-DCC'D':

T F ace-ABB A -T F ace-DCC D = ∆T z (2.81)
Edge-AA', Edge-BB', Edge-CC' and Edge-DD':

T

T Edg e-BC -T Edg e-AD = ∆T x T Edg e-BC -T Edg e-B C = ∆T y T Edg e-B C -T Edg e-A D = ∆T x (2.

84)

Vertice-A, Vertice-B, Vertice-C and Vertice-D:

T V er tice-A -T V er tice-D = ∆T z T V er tice-C -T V er tice-D = ∆T x T V er tice-B -T V er tice-C = ∆T z (2.

85)

Vertice-A', Vertice-B', Vertice-C' and Vertice-D':

T V er tice-A -T V er tice-D = ∆T z T V er tice-C -T V er tice-D = ∆T x T V er tice-B -T V er tice-C = ∆T z (2.86)
Vertice-B and Vertice-B':

T V er tice-B -T V er tice-B = ∆T y (2.87)
where ∆T x , ∆T y and ∆T z are the temperature difference between Face-BCC'B' and Face ADD'A', between Face-ADCB and Face-A'D'C'B' and between Face-ABB'A' and Face-DCC'D', respectively. In order to calculate the components (D ij ) of the composite's equivalent thermal conductivity tensor, the RVE must be numerically simulated three times using suitable magnitudes ∆T x , ∆T y and ∆T z .

In anisotropic case, the thermal conductivity tensor D * has 9 components. ∆T x = 0; ∆T y = 0; and ∆T z = 0 (2.90)

In isotropic case, the tensor D * is diagonal and can be reduced as D * = k * I, where the k * is the scalar thermal conductivity, and our research is based on this isotropic situation.

MUBC

When the material's microstructure is not periodic, we can use the MUBC boundary condition. However, due to the applicability of MUBC, this boundary condition is limited to isotropic and orthotropic thermal conductivity. For example in the same figue 2. Then we have MUBC equations for face:

u i (BCC B ) = " 0 ix • x x (BCC B ) u i (ADD A ) = " 0 ix • x x (ADD A ) u i (ADCB) = ff 0 iy • n y (ADCB) u i (A D C B ) = ff 0 iy • n y (A D C B ) u i (ABB A ) = ff 0 iz • n z (ABB A ) u i (DCC D ) = ff 0 iz • n z (DCC D ) (2.91)
for edges:

u i (AA ) = " 0 ix • x x (AA ) u i (DD ) = " 0 ix • x x (DD ) u i (AD) = " 0 ix • x x (AD) u i (A D ) = " 0 ix • x x (A D ) u i (BB ) = " 0 ix • x x (BB ) u i (CC ) = " 0 ix • x x (CC ) u i (BC) = " 0 ix • x x (BC) u i (B C ) = " 0 ix • x x (B C ) t i (AB) = ff 0 iy • n y (AB) + ff 0 iz • n z (AB) t i (A B ) = ff 0 iy • n y (A B ) + ff 0 iz • n z (A B ) t i (DC) = ff 0 iy • n y (DC) + ff 0 iz • n z (DC) t i (D C ) = ff 0 iy • n y (D C ) + ff 0 iz • n z (D C ) (2.92)
for vertices:

u i (A) = " 0 ix • x x (A) u i (B) = " 0 ix • x x (B) u i (C) = " 0 ix • x x (C) u i (D) = " 0 ix • x x (D) u i (A ) = " 0 ix • x x (A ) u i (B ) = " 0 ix • x x (B ) u i (C ) = " 0 ix • x x (C ) u i (D ) = " 0 ix • x x (D ) (2.93)
In equations 2.91-2.93, " 0 ix is macroscopic strain and ff 0 iy and ff 0 iz are macroscopic stress. In thermal, when we use the MUBC boundary condition to solve for the homogenized heat conductivity coefficient, we usually use the specified temperature as the thermal boundary condition in one direction and the adiabatic condition in the other two directions, where the heat flux is set to 0 (as if it were insulated). Therefore, the boundary conditions of MUBC are:

Face-ADD'A' and Face-BCC'B':

T F ace-ADD A = T lef t T F ace-BCC B = T r ig ht (2.94)
Face-ADCB and Face-A'D'C'B':

q F ace-ADCB • n y = 0 q F ace-A D C B • n y = 0 (2.95)
Face-ABB'A' and Face-DCC'D':

q F ace-ABB A • n z = 0 q F ace-DCC D • n z = 0 (2.
q Edg e-AB • n y = 0 q Edg e-AB • n z = 0 q Edg e-DC • n y = 0 q Edg e-DC • n z = 0 q Edg e-D C • n y = 0 q Edg e-D C • n z = 0 q Edg e-A B • n y = 0 q Edg e-A B • n z = 0 (2.

99)

Vertice-A, Vertice-D, Vertice-D' and Vertice-A':

T V er tice-A = T V er tice-D = T V er tice-D = T V er tice-A = T lef t (2.100)
Vertice-B, Vertice-C, Vertice-B' and Vertice-C':

T V er tice-B = T V er tice-C = T V er tice-C = T V er tice-B = T r ig ht (2.101)
In order to calculate the component of the equivalent thermal conductivity tensor of the composite (D ij ), the model requires three numerical simulations using suitable magnitudes ∆T x , ∆T y and ∆T z as in the case of periodic boundary conditions.

Numerical solvers

When we have completed the above section, we will have a complete system of linear equations and the final result is obtained by solving these equations. The methods for solving algebraic equations are generally divided into direct and iterative solutions.

The direct method is a method which, after a finite number of steps of arithmetic operations, leads to an exact solution of a system of linear equations (if there are no rounding errors in the calculation). It is commonly used to solve systems of low-order dense matrix equations and some large systems of sparse matrix equations (e.g. large systems of band equations).

The iterative method is a method of approximating the exact solution of a system of linear equations step by step using some kind of limit process. Advantages: few storage units, simple programming, the original coefficient matrix remains constant throughout the computation, etc.; Disadvantages: convergence and convergence speed problems exist. It is often used to solve large systems of sparse matrix equations [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. In this work, the Eigen library is used to develop software as it encompasses a variety of direct and iterative solvers.

Conlusion

In this chapter, we focus on the description of some common physical parameters in thermal problems and on the reproduction of isotropic and anisotropic expressions for the Günter scheme, as well as the extensions of the corresponding three-dimensional models.

A summary of homogenization theory is presented, focusing on Hill Lemma, boundary conditions, analytical methods, and loading methods for MUBC, and periodic boundary conditions. Finally, a brief introduction to common matrix solving methods is given.

Chapter 3

Extension of FD models and implementation of boundary condition

Introduction

In this chapter, we focus on our approach, where Section 3.2 introduces our own model developed after being inspired by the Günter scheme, including the 5-point model and the integral model. Secondly, in Section 3.3, we describe in detail our method of loading periodic boundary conditions when using the Günter scheme and the method of making it matrix symmetric and loading various MUBC boundary conditions. Finally, in Section 3.4, we introduce the high-performance computing methods we use, mainly focusing on the CPU method, the two solvers we choose, Eigen and Pardiso, and their comparison. Then, we can get the coefficient a1,...a3:

a 1 = T i;j ; a 2 = T i+1;j ∆x - T i;j ∆x ; a 3 = T i;j+1 ∆y - T i;j ∆y ; (3.3)
The gradients are determined at the centre:

@T @x | i+ 1 2 ;j+ 1 2 = T i+1;j -T i;j ∆x ; @T @y | i+ 1 2 ;j+ 1 2 = T i;j+1 -T i;j ∆y ; @T @x | i+ 1 2 ;j-1 2 = T i+1;j -T i;j ∆x ; @T @y | i+ 1 2 ;j-1 2 = T i;j -T i;j-1 ∆y ; @T @x | i-1 2 ;j+ 1 2 = T i;j -T i-1;j ∆x ; @T @y | i-1 2 ;j+ 1 2 = T i;j+1 -T i;j ∆y ; @T @x | i-1 2 ;j-1 2 = T i;j -T i-1;j ∆x ; @T @y | i-1 2 ;j-1 2 = T i;j -T i;j-1 ∆y (3.4)
We may obtain the heat flux by equation 2.20 and the final description for isotropic using the same procedure 2.21 as for the symmetric scheme : ∇ • -→ q = -(((∆y 2 + ∆x∆y ) T i+1;j +(∆x∆y + ∆x 2 ) T i;j+1 +(-∆y 2 -2∆x∆y -∆x 2 ) T i;j ) D i+1=2;j+1=2 +((∆y 2 -∆x∆y ) T i+1;j +(-∆y 2 + 2∆x∆y -∆x 2 ) T i;j +(∆x 2 -∆x∆y )

T i;j-1 ) D i+1=2;j-1=2 +((∆x 2 -∆x∆y ) T i;j+1 +(-∆y 2 + 2∆x∆y -∆x 2 ) T i;j +(∆y 2 -∆x∆y ) T i-1;j ) D i-1=2;j+1=2 +((-∆y 2 -2∆x∆y -∆x 2 ) T i;j +(∆x∆y + ∆x 2 ) T i;j-1 +(∆y 2 + ∆x∆y ) T i-1;j ) D i-1=2;j-1=2 )=(2∆x 2 ∆y 2 ) (3.5)
By simplifying we get the final relation: With equation 2.24 and 2.21, we can get the 77 Table 3.1: coefficient for T i;j of 5-point isotropic scheme

j Coef i i -1 i i + 1 j+1 (0:25=∆x 2 + 0:25=∆y 2 ) D i-1=2;j+1=2 (-0:25=∆x 2 + 0:25=∆y 2 ) (D i-1=2;j+1=2 + D i+1=2;j+1=2 ) (0:25=∆x 2 + 0:25=∆y 2 ) D i+1=2;j+1=2 j (0:25=∆x 2 -0:25=∆y 2 ) (D i-1=2;j-1=2 + D i-1=2;j+1=2 ) (0:25=∆x 2 + 0:25=∆y 2 ) (-D i-1=2;j-1=2 -D i-1=2;j+1=2 -D i+1=2;j-1=2 -D i+1=2;j+1=2 ) (0:25=∆x 2 -0:25=∆y 2 ) (D i+1=2;j-1=2 + D i+1=2;j+1=2 )) j-1 (0:25=∆x 2 + 0:25=∆y 2 ) D i-1=2;j-1=2 (-0:25=∆x 2 + 0:25=∆y 2 ) (D i-1=2;j-1=2 + D i+1=2;j-1=2 ) (0:25=∆x 2 + 0:25=∆y 2 ) D i+1=2;j-1=2
descriptions for anisotropic :

∇ • - → q = -( `∆x T i;j+1 -∆x 2 T i;j ´k22 i+1=2;j+1=2 + (∆x∆y T i+1;j -∆x∆y T i;j ) k21 i+1=2;j+1=2
+ (∆x∆y T i;j+1 -∆x∆y T i;j ) k12 i+1=2;j+1=2 + `∆y 2 T i+1;j -∆y 2 T i;j ´k11 i+1=2;j+1=2

+ `∆x T i;j-1 -∆x 2 T i;j ´k22 i+1=2;j-1=2 + (∆x∆y T i;j -∆x∆y T i+1;j ) k21 i+1=2;j-1=2 + (∆x∆y T i;j -∆x∆y T i;j-1 ) k12 i+1=2;j-1=2 + `∆y 2 T i+1;j -∆y 2 T i;j ´k11 i+1=2;j-1=2 + `∆x T i;j+1 -∆x 2 T i;j ´k22 i-1=2;j+1=2 + (∆x∆y T i;j -∆x∆y T i-1;j ) k21 i-1=2;j+1=2 + (∆x∆y T i;j -∆x∆y T i;j+1 ) k12 i-1=2;j+1=2 + `∆y 2 T i-1;j -∆y 2 T i;j ´k11 i-1=2;j+1=2
+ `∆x T i;j-1 -∆x 2 T i;j ´k22 i-1=2;j-1=2 + (∆x∆y T i-1;j -∆x∆y T i;j ) k21 i-1=2;j-1=2 + (∆x∆y T i;j-1 -∆x∆y T i;j ) k12 i-1=2;j-1=2 + `∆y 2 T i-1;j -∆y 2 T i;j ´k11 i-1=2;j-1=2 )=(2∆x 2 ∆y 2 )

(3.6)
By simplifying we get the final relation: Table 3.2: coefficient for T i;j of 5-point anisotropic scheme

j Coef i i -1 i i + 1 j+1 -(D22 i+ 1 2 ;j+ 1 2 + D22 i-1 2 ;j+ 1 2 )=(2∆y 2 ) +(D12 i-1 2 ;j+ 1 2 -D12 i+ 1 2 ;j+ 1 2 )=(2∆x∆y ) j (D21 i-1 2 ;j+ 1 2 -D21 i-1 2 ;j-1 2 )=(2∆x∆y ) -(D11 i- ;j+ 1 2 + D11 i-1 2 ;j-1 2 )=(2∆x 2 ) (D11 i+ 1 2 ;j+ 1 2 + D11 i+ 1 2 ;j-1 2 D11 i-1 2 ;j+ 1 2 + D11 i-1 2 ;j-1 2 )=(2∆y 2 ) +(D21 i+ 1 2 ;j+ 1 2 + D12 i+ 1 2 ;j+ 1 2 -D21 i+ 1 2 ;j-1 2 -D12 i+ 1 2 ;j-1 2 -D21 i+ 1 2 ;j-1 2 -D12 i+ 1 2 ;j-1 2 +D21 i-1 2 ;j-1 2 + D12 i-1 2 ;j-1 2 )=(2∆x∆y ) +(D22 i+ 1 2 ;j+ 1 2 + D22 i+ 1 2 ;j-1 2 D22 i-1 2 ;j+ 1 2 + D22 i-1 2 ;j-1 2 )=(2∆x 2 ) (D21 i+ 1 2 ;j-1 2 -D21 i+ 1 2 ;j+ 1 2 )=(2∆x∆y ) -(D11 i+ 1 2 ;j+ 1 2 + D11 i+ 1 2 ;j-1 2 )=(2∆x 2 ) j-1 -(D22 i+ 1 2 ;j-1 2 + D22 i-1 2 ;j-1 2 )=(2∆y 2 ) +(D12 i+ 1 2 ;j-1 2 -D12 i-1 2 ;j-1 2 )=(2∆x∆y )

Integral scheme

Since we cannot find the value of heat flow at the point, we choose the integral method to circumvent this problem by solving for ∇q through the integral property.

Z S ∇qdS = I l q • ndl (3.7)
As shown in Fig. 3.2, we have selected three different integration domains for the calculation, (a) for the square integration domain, (b) for the diamond integration domain, and (c) for the circular integration domain. In the three schemes, a is the parameter that controls the square, rhombus, and circle, a ∈ (0; min(dx; dy )]. We establish the coordinate system with the point (i, j) as the origin and obtain the following temperature expression (T(x,y)) by interpolation:

i,j+1 i,j i,j-1 i+1,j i+1,j-1 i+1,j+1 i-1,j+1 i-1,j i-1,j-1 a a q1 y •n q3 y •n q1 x •n q2 y •n q4 y •n q1 x •n q4 x •n q3 x •n (a) Square i,j+1 i,j i,j-1 i+1,j i+1,j-1 i+1,j+1 i-1,j+1 i-1,j i-1,j-1 a a q1•n q2•n q3•n q4•n (b) Rhombus i,j+1 i,j-1 i+1,j i+1,j-1 i+1,j+1 i-1,j+1 i-1,j i-1,j-1 i,j a a q1•n q2•n q3•n q4•n (c) Circle
" T i;j dx dy - T i;j-1 dx dy - T i-1 ;j dx dy + T i-1 ;j-1 dx dy « xy + " T i;j dy - T i;j-1 dy « y + " T i;j dx - T i-1 ;j dx « x + T i;j " T i;j+1 dx dy - T i;j dx dy - T i-1 ;j+1 dx dy + T i-1 ;j dx dy « xy + " T i;j+1 dy - T i;j dy « y + " T i;j dx - T i-1 ;j dx « x + T i;j " T i+1 ;j dx dy - T i+1 ;j-1 dx dy - T i;j dx dy + T i;j-1 dx dy « xy + " T i;j dy - T i;j-1 dy « y + " T i+1 ;j dx - T i;j dx « x + T i;j " T i+1 ;j+1 dx dy - T i+1 ;j dx dy - T i;j+1 dx dy + T i;j dx dy « xy + " T i;j+1 dy - T i;j dy « y + " T i+1 ;j dx - T i;j dx « x + T i;j (3.8)
The following heat flow expression can be obtained by Eq.3.8:

q1 = 2 6 4 -D i-1 2 ;j-1 2 "" T i;j dx dy - T i;j-1 dx dy - T i-1 ;j dx dy + T i-1 ;j-1 dx dy " y + T i;j dx - T i-1 ;j dx " -D i-1 2 ;j-1 2 ""
T i;j dx dy -

T i;j-1 dx dy - T i-1 ;j dx dy + T i-1 ;j-1 dx dy " x + T i;j dy - T i;j-1 dy " q2 = 2 6 4 -D i-1 2 ;j+ 1 2 "" T i;j+1
dx dy -

T i;j dx dy - T i-1 ;j+1 dx dy + T i-1 ;j dx dy " y + T i;j dx - T i-1 ;j dx " -D i-
With different integration domains, we end up with the same expression which includes coefficient a.

Dirichlet boundary condition

This condition gives the value that the unknown function needs to take on the boundary.

Using Fig.

T i = T lef t i ∈ [0; 3] T j = T r ig ht j ∈ [12; 15]
(3.10)

Neumann boundary condition

The Neumann boundary condition specifies the differentiation of the solution of the differential equation on the boundary, generally given the value of the heat flux at the boundary.

Using Fig. 3.3-(b) as an example, it is known that Φ q = 0 at the top boundary, when the thermal conductivity is isotropic and orthotropic, so for points within the top bound we have @T @x | x;j = 0. We can load Neumann boundary conditions using the following method.

Φ q T i,j T i,j-1 T i,j-2 T i,j+1 T i+1,j T i-1,j+1 T i+1,j+1 T i-1,j T i-1,j-1 T i+1,j-1 D i-1/2,j-1/2 D i+1/2,j-1/2 D i-1/2,j-3/2 D i+1/2,j-3/2 D i-1/2,j-5/2 D i+1/2,j-5/2
(b) Neumann boundary conditions

@T @x | i;j = T i;j+1 -T i;j-1 2∆y = 0 (3.11)
Loading the Neumann boundary conditions and the equation of equilibrium at point 7 in Fig. 3.3-(a) as an example, with the coefficients in Table 2.4, we can obtain the following results using the virtual point method.

coef i-1;j-1 T 2 + coef i-1;j T 3 + coef i-1;j+1 T 3 + coef i;j-1 T 6 + coef i;j T 7 +coef i;j+1 T 7 + coef i+1;j-1 T 10 + coef i+1;j T 11 + +coef i+1;j+1 T 11 = 0 (3.12) Φ q = 0 =⇒ T 3 -T 2 = 0 T 7 -T 6 = 0 T 11 -T 10 = 0 (3.13)
Bringing Eq.3.13 to Eq.3.12, we obtain the final equation:

(coef i-1;j-1 + coef i-1;j+1 )T 2 + coef i-1;j T 3 + (coef i;j-1 + coef i;j+1 )T 6 +coef i;j T 7 + (coef i+1;j-1 + coef i+1;j+1 )T 10 + coef i+1;j T 11 = 0 (3.14)

Neumann alone

In this method, we impose boundary conditions without using virtual points and without loading the equilibrium equations on the boundary. For example, for point (i,j) in Fig. 3.3-(b), we know the value of heat flow on the boundary, and according to Equation 2.20, we use Forward or Backward differences instead of central differences. Second, for the conductivity D i;j , we can use D i-1=2;j-1=2 and D i+1=2;j-1=2 to do the arithmetic average.

Forward difference:

@T @x | i;j = -3T i;j + 4T i+1;j -T i+2;j 2∆x (3.15) Backward difference: @T @x | i;j = 3T i;j -4T i-1;j + T i-2;j 2∆x (3.16)
That is, taking the point 7 in Fig. 3.3-(a) as an example, we have the following equation:

Φ q = 0 =⇒ 3T 7 -4T 6 + T 5 2∆y = 0 (3.17)

Neumann without virtual point 1

This method loads Neumann boundary conditions and equilibrium equations on the boundary without using the virtual points.

∇ • - → q = @q x @x + @q y @y (3.18)
For the part @qx @x , we can find it by the central difference, using q xi-1=2;j and q xi+1=2;j . For the part @qy @y , when the conductivity is isotropic or orthotropic, @qy @y = -D i;j @ 2 T i;j @y 2 . Similarly D i;j can be found by arithmetic averaging. @ 2 T i;j @y 2 can be obtained from Taylor's formula by the following equation

@ 2 T i;j @y 2 = 8T i;j-1 -T i;j-2 -7T i;j -(6∆y ) @T @y | i;j 2∆y 2 (3.19)
Bringing Eq.3.19 into Eq.3.18, combined with Φ q = 0, the final equation is as follows:

∇ • - → q = -D i+1=2;j-1=2 (T i+1;j -T i;j ) + D i-1=2;j-1=2 (T i;j -T i-1;j ) ∆x 2 - (D i+1=2;j-1=2 + D i-1=2;j-1=2 ) 2 8T i;j-1 -T i;j-2 -7T i;j 2∆y 2 (3.20)

Neumann without virtual point 2

The difference between this method and the previous one is the part of solving for @qy @y . For example in Fig. 3.3-(b), in this method, we use the Forward or Backward differences formula to solve for @qy @y using q y i;j , q y i;j-1 , qy i;j-2 . Where q i;j is obtained by Neumann boundary condition, q y i;j-1 and q y i;j-2 are obtained by q y i-1=2;j-1=2 , q y i+1=2;j-1=2 , q y i-1=2;j-3=2 , q y i+1=2;j-3=2 , q y i-1=2;j-5=2 and q y i+1=2;j-5=2 by arithmetic averaging, respectively.

Combined with Φ q = 0, the final equation is as follows:

∇ • - → q = -D i+1=2;j-1=2 (T i+1;j -T i;j ) + D i-1=2;j-1=2 (T i;j -T i-1;j ) ∆x 2
+(-(q y i-1=2;j-1=2 + q y i+1=2;j-1=2 + q y i+1=2;j-3=2 + q y i-1=2;j-3=2 ) +(q y i+1=2;j-3=2 + q y i-1=2;j-3=2 + q y i-1=2;j-5=2 + q y i+1=2;j-5=2 )=4)=(2∆y )

(3.21)

Neumann without virtual point 3

The difference between this method and method 2 is in solving the @qy @y part. For example in Fig. 3.3-(b), in this method, we will use the new formula obtained from Taylor's formula instead of the forward or backward difference formula.

@q @y

| i;j = 8q i;j -9q i;j-1=2 + q i;j-3=2 3∆y @q @y | i;j = -8q i;j + 9q i;j+1=2 -q i;j+3=2 3∆y

Where q i;j is obtained by Neumann boundary condition, q i;j-1=2 and q i;j-3=2 are obtained by q i-1=2;j-1=2 , q i+1=2;j-1=2 , q i-1=2;j-3=2 , q i+1=2;j-3=2 by arithmetic averaging, respectively.

Combined with Φ q = 0, the final equation is as follows:

∇ • - → q = -D i+1=2;j-1=2 (T i+1;j -T i;j ) + D i-1=2;j-1=2 (T i;j -T i-1;j ) ∆x 2 +(-9 (q y i-1=2;j-1=2 + q y i+1=2;j-1=2) 2 + (q y i+1=2;j-3=2 + q y i-1=2;j-3=2) 2 )=(3∆y ) (3.23)

MUBC on pixel models

As Fig. 2.4, and taking point 5 as an example, we can obtain the following:

coef i-1;j-1 T 0 + coef i-1;j T 1 + coef i-1;j+1 T 2 + coef i;j-1 T 4 + coef i;j T 5 +coef i;j+1 T 6 + coef i+1;j-1 T 8 + coef i+1;j T 9 + +coef i+1;j+1 T 10 = 0 (3.24)
The circle and star are loaded with the Dirichlet boundary condition and Neumann boundary condition, respectively. Using the virtual point method as an example, we end up with the following total matrix. It is worth noting that this virtual point method can be used for the 5-point scheme, but is not fully applicable to the Günter scheme. In the steady-state case, due to the specialization of the Günter model (inclusion of corner points), the virtual point method does not apply to the 2-D SUBC boundary condition and leads to matrix singularities that prevent the matrix from being solved using direct methods. In the 3D case, both the MUBC boundary condition and the SUBC boundary condition lead to matrix singularities. To solve this problem we can either choose the transient mode to solve, as detailed in Chapter 4, or use other methods to load the Neumann boundary conditions. model with 36 points, where 20 virtual points and 16 unknown points T i;j , i=0,1,2,3, j=0,1,2,3. Since the final matrix system formed Ax=b has only one index for the vector x consisting of unknown points, we number the unknown points in the way as follows:

PBC on pixel models

p = i * (Ny + 1) + j (3.26)
We take the Günter model as an example, given the temperature at point 0 (T lef t ), the temperature gradient in the x-direction is not 0 (Tgx = 0), and the temperature gradient in the y-direction is 0 (Tgy= 0), where 0-15 numbers are the model point numbers and 0'-15' are the virtual point numbers. The square is loaded with the equilibrium equation introduced by Table 2.4, and taking point 10 as an example, we can obtain the following:

coef i-1;j-1 T 5 + coef i-1;j T 6 + coef i-1;j+1 T 7 + coef i;j-1 T 9 + coef i;j T 10 +coef i;j+1 T 11 + coef i+1;j-1 T 13 + coef i+1;j T 14 + coef i+1;j+1 T 15 = 0 (3.27)
The circle is loaded with the periodic boundary condition, we can obtain the following:

T 3 -T 0 = 0 T 12 -T 0 = T g x * L x T 15 -T 0 = T g x * L x (3.28)
where T g x is the temperature gradient in the x-direction and L x is the length of model.

The star is loaded with the periodic boundary condition and the equilibrium equation and taking point 1 as an example, we can obtain the following using the imaginary point method.

coef i-1;j-1 T 8 + coef i-1;j T 9 + coef i-1;j+1 T 10 + coef i;j-1 T 0 + coef i;j T 1 +coef i;j+1 T 2 + coef i+1;j-1 T 4 + coef i+1;j T 5 + coef i+1;j+1 T 6 = 0 (3.29) T 8 -T 8 = T g x * L x T 9 -T 9 = T g x * L x T 10 -T 10 = T g x * L x (3.30)
Bringing Eq.3.30 to Eq.3.29, we obtain the final equation: where ∆T x = T g x * L x

coef i-1;j-1 T 8 + coef i-1;j T 9 + coef i-1;j+1 T 10 + coef i;j-1 T 0 + coef i;j T 1 +coef i;j+1 T 2 + coef i+1;j-1 T 4 + coef i+1;j T 5 + coef i+1;j+1 T 6 = T g x * L x * (coef i-1;j-1 + coef i-1;j + coef i-1;j+1 ) (3.31)

Symmetrization

Under the boundary conditions of PBC, the matrix is asymmetric due to the presence of virtual grid points. In this case, we cannot use iterative solutions to solve the matrix, which will greatly reduce our solving speed, especially when dealing with 3D problems. So we propose a method to make the matrix symmetric by transforming it. First, we found the reason for the matrix asymmetry, in Fig. 3.4-(a), taking point 4 as an example, point 4 is related to its surrounding points:

coef i-1;j-1 T 2 + coef i-1;j T 0 + coef i-1;j+1 T 1 + coef i;j-1 T 6 + coef i;j T 4 +coef i;j+1 T 5 + coef i+1;j-1 T 10 + coef i+1;j T 8 + coef i+1;j+1 T 9 = T g y * L y * (coef i-1;j-1 + coef i;j-1 + coef i+1;j-1 ) (3.33)
With point 6, we have the quation:

coef i-1;j-1 T 1 + coef i-1;j T 2 + coef i-1;j+1 T 3 + coef i;j-1 T 5 + coef i;j T 6 +coef i;j+1 T 7 + coef i+1;j-1 T 9 + coef i+1;j T 10 + coef i+1;j+1 T 11 = 0 (3.34)
We find that there is T 6 in Eq.3.33, but there is no T 4 in Eq.3.34, which leads to matrix asymmetry. But due to the periodic boundary condition we have T -T 4 = T g y * L y , so bringing this relation to Eq.3.34, we can obtain the relation that satisfies the requirement:

coef i-1;j-1 T 1 + coef i-1;j T 2 + coef i-1;j+1 T 3 + coef i;j-1 T 5 + coef i;j T 6 + coef i;j+1 2 T 7 + coef i;j+1 2 T 4 + coef i+1;j-1 T 9 + coef i+1;j T 10 + coef i+1;j+1 T 11 = -T g y * L y coef i;j+1 2 (3.35)
The next step is the adjustment of the coefficients to make them equal. It is worth noting that there are several special points in the matrix A where there is a multiplicative relationship that cannot reach symmetry requires special treatment, namely point 1, 2, 4, 7, 8, 11, 13, 14. Taking point 1 as an example, if we follow the previous treatment of point 4, there will be a 2-fold relationship that cannot reach symmetry, so we directly replace point 4 in the point 1 relationship with point 7 to obtain the following formula.

coef i-1;j-1 T 8 + coef i-1;j T 9 + coef i-1;j+1 T 10 + coef i;j-1 T 0 + coef i;j T 1 +coef i;j+1 T 2 + coef i+1;j-1 T 7 + coef i+1;j T 5 + coef i+1;j+1 T 6 = T g x * L x * (coef i-1;j-1 + coef i-1;j + coef i-1;j+1 ) + T g y * L y * coef i+1;j-1 (3.36)
To summarize, we take Fig3.4-(b) as an example and transform the matrix to make it symmetric by the following steps:

• Divide by 2 in the overall formula for the points in the circle and the star.

• Add relationship to the points in the square.

• Reset relationship to the points in the star.

Finally we can obtain a symmetric matrix with conductivity D i;j = 1, ∆x = ∆y = 1 as follows:

PBC on voxel models

The method of loading periodic boundary conditions under the voxel model is similar to that of the pixel, except that the 2D pixel model is a 9-point model with periodic boundary conditions loaded on 4 corner points and 4 edges, while the 3D voxel model is a 27-point model with periodic boundary conditions to be loaded on 8 vertices, angles, and 6 faces.

As shown in Fig3.5, it is a 4x4x4 voxel grid with 216 points, where virtual points and 64 unknown points T i;j;k , i=0,1,2,3, j=0,1,2,3, k=0,1,2,3. Since the final matrix system formed Ax=b has only one index for the vector x consisting of unknown points, we number the unknown points in the way as follows:

p = k * (Nx + 1) * (Ny + 1) + i * (Ny + 1) + j (3.38)
In Fig. 3.5, we are given a temperature at point 0 with a non-zero temperature gradient in the X direction (Tgx = 0), a temperature gradient in the Y direction (Tgy= 0), and a temperature gradient in the Z direction (Tgz= 0). The square is loaded with the equilibrium equation introduced by Table 2.6. We can obtain the following equation using point The circle is loaded with the periodic boundary condition, we can obtain the following:

coef i-1;j-1;k-1 T 0 + coef i-1;j;k-1 T 1 + coef i-1;j+1;k-1 T 2 + coef i;j-1;k-1 T 4 + coef i;j;k-1 T 5 +coef i;j+1;k-1 T 6 + coef i+1;j-1;k-1 T 8 + coef i+1;j;k-1 T 9 + coef i+1;j+1;k-1 T 10 + coef i-1;
T 3 -T 0 = 0 T 12 -T 0 = T g x * L x T 15 -T 0 = T g x * L x T 48 -T 0 = 0 T 51 -T 0 = 0 T 60 -T 0 = T g x * L x T 63 -T 0 = T g x * L x (3.40)
The star is loaded with the periodic boundary condition and the equilibrium equation and taking point 1 as an example, we can obtain the following using the virtual point method.

coef i-1;j-1;k-1 T 40 + coef i-1;j;k-1 T 41 + coef i-1;j+1;k-1 T 42 + coef i;j-1;k-1 T 32 + coef i;j;k-1 T 33 +coef i-1;j;k T 9 + coef i-1;j+1;k T 10 + coef i;j-1;k T 0 + coef i;j;k T 1 + coef i;j+1;k T +coef i+1;j-1;k T 4 + coef i+1;j;k T 5 + coef i+1;j+1;k T 6 + coef i-1;j-1;k+1
T 40 -T 40 = T g x * L x T 41 -T 41 = T g x * L x T 42 -T 42 = T g x * L x T 32 -T 32 = 0 T 33 -T 33 = 0 T 34 -T 34 = T 36 -T 36 = 0 T 37 -T 37 = 0 T 38 -T 38 = T 8 -T 8 = T g x * L x T 9 -T 9 = T g x * L x T 10 -T 10 = T g x * L x T 24 -T 24 = T g x * L x T 25 -T 25 = T g x * L x T 26 -T 26 = T g x * L x (3.42)

MUBC on voxel models

Fig. 3.6 shows a 4x4x4 voxel grid, where we are given temperatures on the left and right sides as T lef t and T r ig ht and insulation in the other directions, i.e. Φ q = 0, and we take the same numbering as for the periodic boundary conditions. Here we use the Neumann without virtual point 1 method as an example. The square is loaded with the equilibrium equation introduced by Table 2.6. We can also get Eq.3.39 on point 21. The circle is loaded with the periodic boundary condition, we can obtain the following: ∇ • -→ q = @q x @x + @q y @y + @q z @z @q x @x = q 1+ 1 2 ;2+ 1 2 ;0+

T (0; j; k) = T lef t T (3; j; k) = T r ig ht
1 2 + q 1+ 1 2 ;1+ 1 2 ;0+ 1 2 -q 0+ 1 2 ;2+ 1 2 ;0+ 1 2 -q 0+ 1 2 ;1+ 1 2 ;0+ 1 2 2∆x @q y @y = q 1+ 1 2 ;2+ 1 2 ;0+ 1 2 -q 1+ 1 2 ;1+ 1 2 ;0+ 1 2 + q 0+ 1 2 ;2+ 1 2 ;0+ 1 2 -q 0+ 1 2 ;1+ 1 2 ;0+ 1 2 2∆y @q z @z = - (D 0+ 1 2 ;2+ 1 2 ;0+ 1 2 + D 1+ 1 2 ;2+ 1 2 ;0+ 1 2 + D 0+ 1 2 ;1+ 1 2 ;0+ 1 2 + D 1+ 1 2 ;1+ 1 2 ;0+ 1 2 ) 4
8T 1;2;1 -T 1;2;2 -7T 1;2;0 2∆z 2 (3.44)

Code and HPC

HPC, shortened for high-performance computing, refers to the utilization of pooled computing capacity to execute data-intensive computations that cannot be performed on ordinary workstations, such as simulation, modeling, and rendering. We frequently encounter situations when dealing with various computing problems where a general-purpose computer is unable to complete the task in a reasonable amount of time due to a large amount of computation required, or where the amount of data required and the available resources make the computation impossible. HPC techniques efficiently circumvent these restrictions by utilizing specialized or high-end hardware or by combining the processing capacity of many units. The primary rationale for utilizing HPC is to make use of the rising amount of computing power accessible to individuals as the number of central processing units (CPUs)

and nodes continue to grow. With sufficient computing power, more operations may be performed per unit of time, improving the computational speed of a given model. There are two primary approaches to high-performance computing: conventional CPU clusters and GPGPU. This article focuses mostly on the CPU direction, testing, and comparing current solvers.

Code

First of all, the software we have created is in C++. The software process is illustrated in Fig. 3.7, starting with the reading of the input file, which can be divided into two types.

One type of input file uses the classical RVE included in the software, which contains the model size, the delineation of pixels (voxel meshes), the thermal conductivity of the matrix, and the inclusion, the choice of boundary conditions, and the method of solving the matrix.

The second type of RVE is a more complex calculation, where the input file differs from the first one in that the material number of each pixel (voxel) is given and the thermal conductivity of each material. After reading the input file, the software will create a sparse matrix based on the choice of the input file. And in the next step, the matrix is solved according to the chosen method of joining the matrices, and the temperature gradient, heat flow values, and other parameters are solved by the calculated temperature solution.

Finally, the different parameters of the solution are rooted and the VTK visualization file is created which can be used in ParaView for subsequent analysis.

Solvers

There are many libraries of linear operations available in C++, which is why we have chosen C++ as our development language. So, in the process of software development, we can use them to improve our software. Since our matrix is a sparse, large-scale matrix, we need the solver to have the following characteristics:

• The solution method has high performance. Because the dimension of our sparse matrix is relatively large, the solution's version is more elevated.

• An efficient matrix storage method for storing sparse and intermediate calculation results.

• The solution process is stable, and the result is correct.

In this part, we tested the speed of different solvers of Eigen and Pardiso. My computer processor is an Intel Core i7-8700 with 6 Cores and 12 threads.

Eigen

Eigen is a well-known library of linear operations models. Eigen supports dense matrix and sparse matrix operations and has built-in solvers for both matrices. In addition to its direct solver and iterative solver, Eigen can also connect third-party solvers such as SuperLU, PastixLU, SPQR.

The advantages of Eigen:

• Eigen is free.

• Eigen does not need to be installed. Just load the header file while using it

• Eigen is small.

• Eigen is cross-platform and is currently used on several different operating systems, hardware platforms, and compilers.

• Eigen, compared to some other C++ template libraries, is relatively easy on the compiler.

The different solvers of Eigen

Since our matrix is a large-scale sparse matrix, in Eigen, there are several methods available to solve linear systems when the coefficient matrix is sparse. So we can choose the appropriate solver depending on the type and size of our matrix to get good performance. To speed up the convergence of BiCGSTAB, use preconditioner IncompleteLUT.

Pardiso

The Pardiso [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with pardiso[END_REF][START_REF] Schenk | Pardiso: a high-performance serial and parallel sparse linear solver in semiconductor device simulation[END_REF] software package is a high-performance, robust, memory efficient, and easy-to-use software for solving large asymmetric and unsymmetrical linear systems of equations on shared memory and distributed memory architectures. The solver uses a combination of left and right level 3 BLAS supernode techniques and parallels LU, LDLT, or LLT factorization to improve the performance of sparse sequential and parallel numerical factorization. The Pardiso function interface supports several types of matrices, including real, complex, symmetric, and asymmetric matrices.

The advantages of Pardiso:

• Pardiso is free for academics.

• pardiso shows an excellent efficiency of computation and parallelism.

The different Pardiso solvers

Compared to Eigen, Pardiso has only two solvers. The first is the direct solver based on the LU factorization. The second is the classic iterative CG-based iterative solver. The difference is that Pardiso's direct solver uses the METIS command algorithm to reduce the LU factor filling and can realize multi-threaded calculations.

Conclusion

In this chapter, we first introduce two new models inspired by the Günter model: the In addition, a symmetry method is provided in which the matrix is changed to be symmetric under periodic boundary conditions, which can then be solved iteratively.

Finally, the software encoding process is given, along with the two decoders we use (Eigen and Pardiso). Different solution methods within the two decoders are also presented.

101 heterostructure models at different volume fractions is calculated and compared with the finite element results and the analytical solutions of Voigt and Reuss.

Numerical tests

In this section, various test cases of the diffusion equation are compared with the analytical solution for errors, including homogenized cases, heterogeneous cases with anisotropic diffusion without source term, and with source term.

Homogenized case with Scalar coefficient and source term

We consider a simple steady diffusion problem with scalar coefficient given in [START_REF] Chamarthi | First order hyperbolic approach for anisotropic diffusion equation[END_REF] with the following parameter: With the exact soulution is given by:

D = 2
T exact (x; y ) = x y ; x; y ∈ [1; 2]; (4.4)

Heterogeneous case with anisotropic tensor coefficient and source term

Next, we consider a simple stable diffusion problem with the following anisotropic tensor coefficients and source term given in [START_REF] Lou | Reconstructed discontinuous galerkin methods for hyperbolic diffusion equations on unstructured grids[END_REF]. With the exact solution is given by: T exact (x; y ) = 1 -tanh(100((x -0:5) 2 + (y -0:5) 2 ); x; y ∈ [0; 1];

D =
And source term is given by: for scheme 5-point. This proves that scheme Günter is order 2 on the heterogeneous case with anisotropic tensor coefficient and scheme 5-point is less than order 2.

S = 100(1 -A 2 )(200AB + C) ( 4

Heterogeneous case II with anisotropic tensor coefficient and source term

As in the previous case, we found that the 5-point model is smaller than second order when the anisotropic tensor is coefficient. So we choose the following case to verify it again. values varying from 10 0 to 10 9 . This indicates that the rate of convergence for scheme Günter is independent of anisotropy, but the convergence rate of the 5-point scheme is only related to anisotropy when ' < 2.

1D problem

First, we take a 1D model with PBC (T 0 = 0; T g x = 1) as an example to compare our method with FEM. As shown in Fig. 4 As shown in Table 4.1, by comparing the results it was found that the two methods obtained exactly the same results. This also proves the feasibility of our method in onedimensional problems.

Two-dimensional pixel models with PBC

In this section, we compare FDM results for periodic boundary conditions with the analytical solution, FEM and FEM+pixel, respectively. The FDM approach is validated by comparing the temperature field, temperature gradient field, heat flow field, and equivalent thermal conductivity and convergence rate. where Ω denotes a field consisting of a heterogeneous material(the matrix Ω \ Ω f ib and the fiber Ω f ib ). A temperature of 0 is imposed at point 0, and then a periodic boundary condition is imposed on the boundary with the following thermal boundary value problem:

Layered composite with analytical results

8 > > > > > > > > < > > > > > > > > : div - → q = 0 in Ω T = 0 on point 0 T Γ l -T Γr = 1 on Γ r , Γ l T Γt -T Γ b = 0 on Γ t , Γ b (4.14)
With isotope thermal conductivity: paring the results computed by our method exactly match the analytical method, which proves the accuracy of our method. It is worth noting that when dividing the pixels, it is necessary to make the boundaries of the pixels coincide with the fiber boundaries.

D = 8 > < > : 1 in Ω \ Ω f ib 100 in Ω f ib

FDM/FEM comparison

In this section, for a static heat conduction problem, we compare the distributions of temperature (T), temperature gradient (Tgx), and heat flow (Hfx) calculated by FDM and FEM in the corresponding meshes, and compare the thermal conductivity cases for both isotropic and anisotropic cases. The results of the FEM part are calculated by Comsol software. For periodic boundary conditions, Comsol software has an existing module, the principle is T (x 0 ) -T (x 1 ) = ∆T , for the two-dimensional model we can specify the corresponding two edges as periodic, in order to ensure the accuracy of the results, we choose to mesh one of the edges first and then copy the mesh to the corresponding edge to ensure that the mesh points correspond to each other after the mesh. Finally, we choose ParaView to visualize the temperature field, temperature gradient, and heat flow to compare the two methods and calculate the above values on the path middle (y=0.5)

and path diagonal (y=x) respectively.

In two-dimensional situations, we consider a domain defined by simple geometric shapes such as circular inclusion, cross, ellipse, and bcc. For the choice of thermal conductivity, we choose D f =D m = 100 instead of specific material parameters for the convenience of calculating isotropic conditions. anisotropic, the final equivalent thermal conductivity is anisotropic. Secondly, the equivalent thermal conductivity obtained for the FDM is also very close to the values of the FEM and FEM+pixel for the case of 2025 cells.

D m = 2

Anisotropic cross case

Isotropic disk case

Fig. 4.12 shows the convergence analysis of one disk inclusion of radius equal to 0.3. Although the slope difference between the two is not significant, the intercept of FDM is smaller than that of FEM+pixel, implying that FDM has a faster convergence rate than FEM+pixel. 4.5 shows the homogenised thermal conductivity obtained with the above mesh and by comparison we find that the FDM is still slightly better than the FEM+pixel.

secondly, the final equivalent thermal conductivity is obtained isotropic when the disc model's input thermal conductivity is isotropic. Secondly, when comparing the path data for the temperature gradients in (d) and (e), the changes occur mainly at the junction of the inclusion and matrix phases. However, among the differences between the three methods, the data for FDM are mainly between FEM and FEM+pixel, indicating that FDM is still better than FEM+pixel. In the comparison of the path data for heat flow in (e) and (f), we observe that the error components of the three methods drop as the grid density is increased relative to the isotropic model.

Anisotropic disk case

However, FEM+pixel retains a significant inaccuracy near the inflection point, whereas FDM does not. and by comparison we find that the FDM is still slightly better than the FEM+pixel.

secondly, the final equivalent thermal conductivity is obtained anisotropic when the disc model's input thermal conductivity is anisotropic. two methods, FDM and FEM+pixel, are not very different, the slope of FDM is slightly smaller than that of FEM+pixel, so FDM is superior to FEM+pixel. produces outcomes that are more similar to FEM than those of FEM+pixel. FEM+pixel has a more significant inaccuracy outside the inflection point than FDM. In the comparison of the path data for the heat flow, it can be seen that both FDM and FEM+pixel approaches are similar.

Table 4.7 shows the equivalent thermal conductivity of the above models, and according to the results, it can be found that for the elliptical model, the final equivalent thermal conductivity results are orthogonal when the input thermal conductivity is isotropic. Secondly, based on the analysis of the errors, it can be seen that the FDM is slightly better than the FEM. Table 4.8 shows the equivalent thermal conductivity of the above models. When the input thermal conductivity for the elliptical model is anisotropic, we discover that the final equivalent thermal conductivity derived is also anisotropic. Secondly, the error of the FDM is still smaller than that of the FEM+pixel. The comparable thermal conductivity of the models mentioned above is compared in Table 4.9. The comparison demonstrates that when the input thermal conductivity is isotropic, the resulting equivalent thermal conductivity is similarly isotropic, and the FDM value is higher than the FEM+pixel value. it can be seen that even though the intercept of FEM+pixel is smaller, its slope differs from that of FDM, so it is still FDM that has a better convergence rate than FEM+pixel.

The comparison of overall trends of the anisotropic case of one disk and four sectors 0 2 0 0 0 0 4 0 0 0 0 6 0 0 0 0 8 0 0 0 0 Num be r o f DOF After increasing the number of cells, the data fit well on all paths, but FEM+pixel still has large discrepancy values at the inflection points.

The following Table 4.10 compares the comparable thermal conductivity of the models mentioned above. The comparison demonstrates that while the input thermal conductivity is anisotropic, the resulting equivalent thermal conductivity is similarly anisotropic, and the FDM value is still greater than the FEM+pixel value.

Two-dimensional pixel model with MUBC

Due to the limitations of mixed uniform boundary conditions (MUBC), these boundary conditions are not fully applicable to anisotropic problems, so for the model shown in Fig. 4.7, we have only compared the FDM, FEM, and FEM+pixel methods for the isotropic case and the anisotropic case, and select the following data as input thermal conductivity.

For the isotropic thermal conductivity we choose the following: can be seen that the findings for the isotropic case are the same for the mixed uniform and periodic boundary conditions. Thirdly, the overall trend is the same for all three ways, regardless of whether it is the temperature field, the temperature gradient field, or the heat flow field, with the exception of the cross model's four inflection points.

D m = 2
The orthotropic instance of one cross inclusion with a side length of 0.3 is illustrated in According to this table, we find that in the isotropic case the results for the mixed uniform boundary condition are identical to those for the periodic boundary condition. We can see that the overall trend is comparable in the orthotropic and isotropic cases, and the results derived are similar, i.e. the three methods produce the same variance in the temperature field, and the temperatures on the paths are fairly similar. For the variation of the temperature gradient field and the heat flow field, the differences also exist only at the four inflection points with the cross model. The temperature field variation (see (a), (b), and (c)), the temperature gradient field variation (see (d), (e), and (f)), and the heat flow field variation (see (g), (h), and (i)) trends reveals that while all three methods exhibit a similar overall trend, the temperature gradient field and heat flow field in the region between the inclusion phase and matrix phase with the FEM method exhibit differences, primarily due to thermal conductivity discontinuities.

Orthotropic cross case

We used data from the paths with y = 0. The equivalent thermal conductivity of an isotropic disk case with MUBC is shown in Table 4.14; we obtained a similar conclusion with PBC, namely that when the input thermal conductivity is isotropic, the final equivalent thermal conductivity is also isotropic. FDM produces superior results to FEM+pixel.

Orthotropic disk case

The convergence of one disk inclusion of radius 0. 4. [START_REF] Reuss | Berechnung der fließgrenze von mischkristallen auf grund der plastizitätsbedingung für einkristalle[END_REF], we can see that the convergence rates for the various loading Neumann's methods are quite similar and that for the orthotropic disc model, FDM's method is always superior to FEM+pixel. Whereas for the disc model, the ultimate convergence rate is slightly variable for varied input thermal conductivities (isotropic versus orthotropic). The temperature field, temperature gradient field, and heat flow field for an orthotropic disk model are shown in Fig. 4.38. The FDM and FEM+pixel have 3600 cells each, while the FEM has 4042 cells. The comparison shows that the overall trend is similar to the isotropic case, i.e. the temperature field follows exactly the same trend, the overall trend in the temperature gradient field and the heat flow field is approximately the same but there are differences in the region between the inclusion phase and the matrix phase. 4.38, demonstrating that the three approaches agree well when comparing the temperature path. In the comparisons of temperature gradient and heat flow, it can be observed that as the number of cells rises, the FEM+pixel values converge to the FDM values. However, the FEM+pixel still exhibits a considerable value difference at the inflection point.

The equivalent thermal conductivity of the above models is compared in Table 4.16 using the MUBC boundary conditions. The comparison demonstrates that when the input thermal conductivity is orthotropic, the disc model also has an orthotropic equivalent thermal conductivity. Secondly, FDM is superior to FEM+pixel in terms of inaccuracy.

Isotropic ellipse case

The convergence analysis with the long axis of the ellipse equal to 0.6 with isotropic thermal conductivity is shown in Fig. (c)-(f) illustrate the temperature gradient and heat flow on the paths, respectively, where the differences between the three methods are most pronounced at the intersections of the inclusion and matrix phases, and where the FEM+pixel can exhibit significant differences, with the overall comparison indicating that the FDM is superior to the FEM+pixel. Table 4.18 shows the equivalent thermal conductivity of an isotropic ellipse case with MUBC. We achieved the same conclusions as to the PBC, namely that when the elliptical filler model's input thermal conductivity is isotropic, the final equivalent thermal conductivity is orthotropic. It can be seen from the error comparison that the FDM value is 157 The convergence analysis with the long axis of the ellipse equal to 0. The overall comparison does not differ significantly from the isotropic model described above. The overall trends in the results for the temperature, temperature gradient, and heat flow fields obtained by the three methods are generally the same, but there are still differences between the inclusion and matrix phases. as the number of cells increases, the FEM+pixel value eventually converges to the FDM value. However, there are still significant differences at the inflection points. The overall comparison is still better for FDM than for FEM+pixel.

Table 4.20 shows the equivalent thermal conductivity of an orthotropic ellipse case with MUBC. The findings indicate that while the input thermal conductivity is orthotropic, the end equivalent thermal conductivity is similarly orthotropic. The FDM is still superior to the FEM+pixel when the cell count is increased. we observe that the overall trend changes are the same for all three methods, whether in the temperature field, the temperature gradient field, or the heat flow field, but there are differences in the region between the inclusion phase and the matrix phase for the temperature gradient field and the heat flow field, due to discontinuities in thermal conductivity.

For the path data comparison (see Fig. Table 4.22 shows the equivalent thermal conductivity of an isotropic One disk and four sectors case with MUBC. We obtain the same results as the periodic boundary condition, in that when the input thermal conductivity is isotropic, the final equivalent thermal conductivity is also isotropic. Secondly, the FDM results are significantly better than the FEM+pixel results for the one disk and four sectors model.

Orthotropic one disk and four sectors case

The convergence analysis of one disk and four sectors inclusions of radius r= 0. 4. [START_REF] Kröner | Berechnung der elastischen konstanten des vielkristalls aus den konstanten des einkristalls[END_REF]. By comparing the coefficients, we discover that the slope of the trend line is somewhat smaller than in the isotropic case and that the intercept is also slightly more significant, indicating that the rate of convergence is slightly slower in the orthotropic case than in the isotropic case.

In comparison to FEM+pixel, while the slope of the trend line of FDM is slightly less than that of FEM+pixel, its intercept is much less, implying that FDM has a higher convergence rate than FEM+pixel. We opted to compare data on the y = 0:5 and y = x pathways (see Fig.

Evolution of the equivalent thermal conductivity

In this section we will compare the variation of the equivalent thermal conductivity with volume fraction for various 2D models with different boundary conditions. For the different The comparison shows that the equivalent thermal conductivity of all models becomes more significant as the volume fraction of fiber increases, which is logical. Secondly, the equivalent thermal conductivity results obtained for different boundary conditions are between the estimates given by Voigt and Reuss. Additionally, it can be seen that the solutions obtained for the KUBC boundary conditions are more significant than the results for PBC and MUBC (more pronounced for the cross and Disk and 4 sectors model and less so for Disk and Ellipse model), which is also consistent with the results from Suquet [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] and Ostoja-Starzewski [START_REF] Ostoja-Starzewski | Material spatial randomness: from statistical to representative volume element[END_REF] mentioned in their paper. Finally, in the isotropic case, the PBC results are almost indistinguishable from the MUBC results.

Three-dimensional model

In this section, we will compare FDM, FEM, and FEM+voxel under the 3D model. We have chosen mainly spherical and cylindrical fillings as representatives for the convergence analysis of FDM and FEM+voxel, the comparison of temperature fields, temperature gradient fields, heat flow fields, and equivalent thermal conductivity for the three methods FDM, FEM, and FEM+voxel. Since there is no difference between the results for PBC and MUBC in an isotropic medium, the data can be seen as either of these two boundary conditions in the next isotropic comparison. The data prove that our method and the analytical solution give identical results. Table 4.25 shows the equivalent thermal conductivity of 3D layered composite. The pairwise results show that for different volume fractions of the inclusion, our method is exactly the same as the analytical solution found by Voigt and Reuss. 4. [START_REF] Sharma | Size-dependent eshelby's tensor for embedded nanoinclusions incorporating surface/interface energies[END_REF] The comparison shows that the rate of convergence varies slightly with the radius of the sphere, but in most cases, the slope of the trend line is smaller for FDM than for FEM+voxel, which means that FDM converges faster. In a small number of cases, although the slope of the trend line is similar to that of FEM+pixel, the intercept is much smaller than that of FEM+voxel, which means that the convergence speed of FDM is ultimately faster. shows that the overall images of the temperature field, temperature gradient field, and heat flow field of the three methods match well and do not differ significantly. On the cross-section xOz, the temperature results agree well. For the temperature gradient and heat flow values in section xOz, it can be seen that the overall trend is the same, but due to the voxel grid, there are differences between the FDM and FEM+voxel and FEM values in the region between the inclusion and matrix phases. However, comparing FDM and FEM+voxel, the differences are smaller.

3D Layered composite model

(d) (e) (f) (g) (h) (i) (j) (k) (l) (m) (n) (o) (p) (q) (r)

Cylinder

Fig. 4.59 shows the convergence analysis of the cylindrical inclusion model, which includes the variation of the ETC with increasing DOF, and the error analysis at different R/L using the FEM values as reference. (a), (c), (e), (g), (i), (k), and (m) show the variation of ETC at different R/L. It can be seen that the ETC of the FRM+voxel method gradually converges to the value of FDM as the DOF increases, and the value of FDM is always between that of FEM and FEM+voxel. The error analysis in (b), (d), (f), (h), (j), (l) and

(n) shows that the FDM values are almost always below the FEM+voxel values, which means that the errors in the FDM are all smaller than those in the FEM+voxel. In Table 4.27, which shows the trend line parameters for both FDM and FEM+voxel methods at different R/L, it can be seen that the slope of FDM is smaller than that of FEM+voxel in most cases. Even though the slope of FDM is similar to that of FEM+voxel in the remaining cases, the intercept of FDM is smaller than that of FEM+voxel, which means that the convergence rate of FDM is better than that of FEM+voxel for different R/L cases. Secondly, we can also see that the convergence speed of the same method is slightly different at different R/L. The comparison of the Dirichlet boundary conditions imposed in the x-direction in Fig4.60

shows that both the overall trend and the cross-sectional temperatures match well in the comparison of the temperature fields of the three methods, with no significant differences.

(d) (e) (f) (j) (k) (l) (m) (n) (o) (p) (q) (r) (g) Tgz(°C/m) (h) (i)
For the comparison of temperature gradients and heat flow values, it can be seen that the overall trend is consistent, but due to the voxel grid, the two methods, FDM and FEM+voxel, still differ from the FEM values in the region between the inclusion and matrix phases. However, comparing FDM and FEM+voxel, the differences are smaller. The Dirichlet boundary condition in Fig4.61 is applied in the z-direction and the comparison shows that in this direction the three methods agree very well in both the temperature and temperature gradient fields and the heat flow field, the only remaining difference being that in the heat flow field comparison the surface of the cylinder is serrated rather than circular as in the finite element method due to the voxel. The overall comparison between FDM and FEM+voxel is identical. As low-temperature sintering of silver nanoparticles is becoming a reliable technology for solder die attachment of electronic and optical components [START_REF] Wang | Numerical modeling of low-temperature and low-pressure sintering of silver microparticles based on surface and grain boundary diffusion mechanisms[END_REF], it becomes necessary to develop numerical approaches with the ability to predict correctly the thermal performance of such materials. Therefore, in this section, our method is applied to the sintered silver model and compared with the equivalent thermal conductivity obtained from the FEM.

Geometrical model

We consider for the sake of simplicity three typical types of lattice structures: simple cubic crystal form (SC), body-centered cubic crystal form (BCC), and face-centered cubic crystal form (FCC), and their conjugated structures. The theoretical geometric models are shown in Fig. 4.63, where the spherical particles in Fig. 4.63a-c can represent either atomic/molecular clusters at mesoscopic scale, or stacked granules at the macroscopic scale.

The spheres keep in touch or overlap with each other (marked in color) to meet practical situations: for example, in a mechanical system with large number of loose particles, they are deformed due to squeezing actions and therefore form contact surfaces; in a heat transfer system, the clusters of sintered silver will fuse together and penetrate into each other with the increase of sintering time, and the porosity of the overall structure reduces as the overlap area increases. Consequently, the use of these simple basic models can provide references for complex structures. In our case, these ideal arrangements of particles offer good approximations for real powder compacts prepared by sintering and serve to illustrate different packing factors of the particles with more or less porous silver particle networks [START_REF] Wang | Numerical modeling of low-temperature and low-pressure sintering of silver microparticles based on surface and grain boundary diffusion mechanisms[END_REF].

It is noted that the ideal geometrical models that we considered possess the property of spatial periodicity. The representative cell can be obtained by cutting out along the orthogonal edges illustrated in Fig. 4.63d-f. In the subsequent study, the side length of the cell is taken to be the unit length L SC = L BCC = L F CC = 1. They are composed of two phases: a solid part consisting of silver contacting balls and a remaining void part filled with air. To calculate the equivalent thermal conductivity under normal conditions, it is necessary to ensure that the solid phases are connected as a single entity. However, considering only the solid phases will lead to the porosity of the studied structures varying only within a small range. Hence, in order to verify the robustness of our developed algorithm, their conjugated structure are also considered (see Fig. 4.63g-i), i.e., the gas phase and solid phase are swapped to obtain a broader range of porosity.

Results and comparisons

Convergence analysis The convergence analysis is applied to a case of the unit BCC structure, where the sphere radius are all equal to 0.46, thus resulting in a volume fraction of 79.89%. The boundary conditions we used are PBC such that ∆T x = 1, ∆T y = 0, and ∆T z = 0. We fix the temperature of one vertex on the left side to 19 • C so that the right side is 20 • C (room temperature). The thermal conductivity of silver is 429 W

• m -1 •K -1
and that of air is 0.0257 W• m -1 •K -1 . To study the convergence of thermal conductivity, different refinements of mesh-grid are used, i.e. from 10 thousand to 2 million elements.

The results are shown in Table 4.28 and in Fig. 4.64. According to Fig. 4.64, we find that firstly, the equivalent thermal conductivity exhibits some fluctuations but they are reasonable. These small fluctuations can be explained by the nature of voxel meshing which results in slight variations of the volume fraction of silver when the mesh density is increased. Secondly, the convergence rate of the finite difference method is very close to that of FEM, which shows the stability and robustness of the FD algorithm. Finally, when we compare the equivalent thermal conductivity obtained for a number of elements greater than 500 thousand elements, the results remain nearly unchanged (226 W• m -1 •K -1 by FDM, 224 W• m -1 •K -1 by FEM) and the error between the two methods is less than 1%.

Thereafter, we will use the FD system that we developed to calculate and compare the temperature field and the heat flow field for different geometrical models at different porosities. Consequently, the volume fraction of silver can take values in a wider range (i.e. from 0 to 1). We can see that the results obtained from FDM for different mesh densities and from FEM are superimposed very well at both low and high fraction values. We also introduce Maxwell's analytic formula [START_REF] Maxwell | A treatise on electricity and magnetism[END_REF]:

Simulation results of local fields

D eq D m = 1 + 3ffi ( D f +2Dm D f -Dm ) -ffi (4.22)
where D m represents the thermal conductivity of matrix, D f the thermal conductivity of fillers, D eq the equivalent conductivity of the whole structure, and ffi is the volume fraction of fillers. This formula is valid only in the case of low ffi according to Pietrak et al. [START_REF] Pietrak | A review of models for effective thermal conductivity of composite materials[END_REF] (under about 25%). It corresponds to the situation when the volume fraction of silver is above 75%. When comparing all the three results, it turns out that the FDM and the FEM curves can fit well with the Maxwell's formula not only in the expected fraction, but also between 50% and 75%. This proves the accuracy and effectiveness of the finite difference algorithm, even in a sparse mesh density of 45 3 . one can find that the value of this property for the SC structure is larger than the ones for the BCC/FCC structures for an identical fraction of silver (or porosity). It indicates that changing the inner organization pattern of a structure can significantly impact its thermal property.

Fig. 4.73 compares the equivalent thermal conductivities of the non-conjugated models with those of the corresponding conjugated models. It can be found that the results for the conjugated structures are larger than the ones for the non-conjugated structures for the three models, when making the comparison at the same porosity. This confirms that the conjugated structures are more conductive. By altering the topology or morphology of a structure, the thermal property can thus be changed greatly. conclusion in the previous section that the voxelized mesh does not simulate the interfaces well when the material with a larger thermal conductivity has a narrow connection area. By comparing the equivalent thermal conductivity in Case 2 and Case 3, it can be found that the FD simulation of the conjugated structure is better than that of the normal structure when the silver fraction remains the same. The relative errors of the three cases are all below 6%. , it is shown that the overall trend is the same, but in some junctions between different silver spheres, the finite element method is slightly better than the finite difference method, which is caused by the voxelized mesh, as the silver sphere connections inside the random system are more complex and require a more dense mesh to approximate finely the curved interfaces at these locations. diagonal preconditioner and initial guess, respectively. We can find that although Eigen can use multi-threads, it does not parallelize very well, and with 4 threads the parallel speed is the fastest. In Eigen, the direct solver cannot use multi-threads. Then, we will compare the speed of the direct solver and the different iterative solvers under 4 threads. With Fig. 4.78, we can see that when the DOF is less than 10 5 , the speed of the iterative solver is higher than that of the direct solver. As the DOF increases, the speed of the direct solver will get closer to the iterative solver, but at the same time, the direct solver will hit the upper limit of the calculation due to the excessive DOF. At the same time, we found that the preprocessor had little effect on the iterative solver. The initial assumption had a significant influence on the iterative solver; we only divide the temperature evenly at each point, increasing the speed by about 30%. It also shows that if the initial value closer to the result can be found, the rate of the iterative solver will increase even more.

Conclusion

In this chapter, we present our numerical results, comparing the results obtained by FDM, FEM, and FEM with pixel (voxel) for different models with periodic boundary conditions and mixed homogeneous boundary conditions. Through extensive comparisons, it is demonstrated that our FDM method outperforms FEM with pixel (voxel) in terms of convergence speed and that the results are more closely matched to those of FEM.

Applying our method to solving the equivalent thermal conductivity of sintered silver, we use a spherical filler model to model the sintered silver structure, consisting of two phases: silver and air. Three typical geometries (SC, BCC, FCC models) and stochastic models are analyzed. Comparisons are carried out with the finite element method. The following conclusions can be summarized. First, our proposed algorithm yielded comparable results with the FEM, whether dealing with the classical SC, BCC, FCC models, or complex stochastic models. Second, the calculation of the equivalent thermal conductivity of the conjugated 206 structures is more accurate than the typical structures for all geometrical models that we used. Finally, the FDM is more effective in dealing with conjugated structures because it is less sensitive to the density of mesh. This point indicates that the developed FDM is suitable for handling this type of material.

Finally, in the comparison between Eigen and Pardiso, it can be seen that Eigen is generally parallelized, in contrast, Pardiso can be well parallelized and Pardiso is much faster than Eigen as the DOF increases.

207

Chapter 5

Complementary research 5.1 Transient model

Due to the peculiarities of our finite-difference model, the matrix will be singular if the virtual point approach is used to load the Neumann boundary conditions. We can either utilize another method provided in the preceding section to load the Neumann boundary conditions or calculate using a transient model to fix this concern.

In comparison to the transient state, the steady-state has the advantage that it can be solved directly, and there is no time term truncation mistake, but the drawback is that the solution's uniqueness is in doubt. On the contrary, the benefit of the transient state is that it does not need consideration of the solution's existence. Generally, the matrix properties are superior to those of the steady-state matrix. The downside is that iterative calculation requires a large number of time steps, and there is a time term truncation error.

The two spatial dimensional thermal conduction equation with isotropic conductivity (2D) can be written as: @T (x; y ; z; t) @t = D[ @ 2 T (x; y ; z; t) @x 2 + @ 2 T (x; y ; z; t) @y 2 ] (5.1)

where D is conductivity. where n is the index of the time step,fl x = D∆t=∆x 2 ; fl y = D∆t=∆y 2 , ‹x 2 T n i;j = T i+1;j+1 + 2T i+1;j + T i+1;j-1 -2T i;j+1 -4T i;j -2T i;j-1 + T i-1;j+1 + 2T i-1;j + T i-1;j-1 ‹y 2 T n i;j = T i+1;j+1 + 2T i;j+1 + T i-1;j+1 -2T i+1;j -4T i;j -2T i-1;j + T i+1;j-1 + 2T i;j-1 + T i-1;j-1

Euler explicit scheme

Note that there is only one unknown in each equation concerning point (i,j). Therefore, there is no matrix inversion needed. The truncation error is O(∆t) + O(∆x 2 ).

Stability analysis

Substitute in the equations the solution of the form T n j =n exp(ix i k). The advantage of the Euler explicit scheme is that it does not require matrix inversion, which significantly reduces calculation and calculation time. But the disadvantage is that the limitation of stability conditions makes it too many iterations when solving significant problems, resulting in too long solution time. The advantage of this method is that due to unconditional stability, there is no limit in the step size. However, the disadvantage is that when solving large-scale problems, the space and time occupied by matrix inversion are too large.

The temperature variation results of the disk model obtained using the Euler implicit Although the Crank-Nicolson (CN) scheme is unconditionally stable and it is considered to have very high numerical accuracy and has improved the convergence speed compared to the implicit algorithm. However, this method also needs to solve a vast spare irreducible matrix. Consider a two-dimensional grid of size (nx+1)by (ny+1), the number of variables required by the system is N=(nx+1)(ny+1), which involves a matrix of N*N size to solve,so that the CN scheme can hardly be used for practical problems.

The temperature variation results of the disk model obtained using the Crank-Nicolson scheme are shown in Fig. 5.3. The structure has 900 pixel cells. When the step size is set to 0.002, 426 iterations are required to reach a tolerance of 1e-6.

ADI method

As the heat conduction problem changes from one-dimensional to multi-dimensional, the order of the formed algebraic equations has increased by one or two orders of magnitude.

How to solve these algebraic equations economically and effectively has become an important issue. Although the explicit format does not have the problem of solving algebraic equations, its stability conditions greatly limit the time step. When solving large-scale problems, the limited time step causes too many iterations, which seriously affects the solution time. When adopting the implicit format and semi-implicit format of the secondorder intercept, the value of each node must be solved at the same time in each time layer.

For three-dimensional problems, this consumes a lot of computer memory and calculation time. To alleviate this difficulty, Peaceman, Rachford [START_REF] Peaceman | The numerical solution of parabolic and elliptic differential equations[END_REF], and Douglas [START_REF] Douglas | A general formulation of alternating direction methods[END_REF] proposed the ADI method. Basiclly, when using this kind of method to solve the problem of unsteady heat conduction, a three-dimensional problem can be transformed into three concatenated one-dimensional implicit format problems.

Peaceman-Rachford method

First, the factorized form of 5.6 can be written to be:

( (5.9)

After simplification, 5.9 can be written as:

( Gunn method. The structure contains 900 pixel cells, and when the step size is set to 0.002, 412 iterations are required to achieve a tolerance of 1e-6.

Although the transient state can solve the matrix dissatisfaction, the final result can be obtained. However, we spent more time on iteration. Compared with our ultimate goal of solving homogenization, it is evident that steady-state calculation is more suitable for us and more effective. But ADI's idea is also used for steady-state. This idea can be a more clever parallel algorithm to save running time.

Alternating direction iteration

The ADI method to solve unsteady heat conduction is discussed from the perspective of a discrete format. ADI constructs a way of solving algebraic equations of multi-dimensional problems that combine direct and iterative methods. In other words, we can use the alternate direction implicit form, which solves the non-steady-state problem, to solve the steady-state problem, which is essentially the same. In the non-steady-state situation, the one-time level is advanced, which is equivalent to completing a level of iteration in the steady-state problem. In conventional iterative methods, such as Jacobi iteration or Gauss-Seidel iteration, the scanning direction in each iteration remains unchanged, that is, row by row or column by column. After scanning the entire field, one iteration is completed.

Repeat this method until the iteration is completed after convergence. If scanning in alternate directions is used, the convergence speed can often be accelerated. One scan is performed row by row (or column by column), and then a scan is performed column by column (or row by row), and two full-field scans form one iteration. In this way, we can use the alternating direction iteration method to solve the steady-state problem and effectively improve the solution speed and reduce the computer memory consumed by the solution. (16∆x∆y ) + k33 i+ 1 2 ;j+ 1 2 ;k+ (16∆z∆z) + k33 i+ 1 2 ;j+ 1 2 ;k-1 2 (16∆z∆z) + k33 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆z∆z) + k33 i+ 1 2 ;j-1 2 ;k-1 2 (16∆z∆z) + k32 i+ 1 2 ;j+ 1 2 ;k+ 1 2 (16∆y ∆z) -k32 i+ 1 2 ;j+ 1 2 ;k-(16∆y ∆z) -k32 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆y ∆z) + k32 i+ 1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k31 i+ 1 2 ;j+ 1 2 ;k+ 1 2 (16∆x∆z) + k31 i+ 1 2 ;j+ 1 2 ;k-1 2 (16∆x∆z) -k31 i+ 1 2 ;j-1 2 ;k+ (16∆x∆z) + k31 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆z) )T i+1;j;k +( k13 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆z) -k13 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆z) + k12 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆y ) + k12 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k11 i+ 1 2 ;j-1 2 ;k+ (16∆x∆x) -k11 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆x) -k23 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆y ∆z) + k23 i+ 1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k22 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆y ∆y ) -k22 i+ 1 2 ;j-1 2 ;k-(16∆y ∆y ) + k21 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆y ) + k21 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) + k33 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆z∆z) + k33 i+ 1 2 ;j-1 2 ;k-1 2 (16∆z∆z) + k32 i+ 1 2 ;j-1 2 ;k+ (16∆y ∆z) -k32 i+ 1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k31 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆z) + k31 i+ 1 2 ;j-1 2 ;k-

-k31 i-1 2 ;j-1 2 ;k+ 1 2 (16∆x∆z) + k31 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆z) )T i;j;k +( k13 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆z) -k13 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆z) -k13 i-1 2 ;j-1 2 ;k+ 1 2 (16∆x∆z) + k13 i-1 2 ;j-1 2 ;k-(16∆x∆z) + k12 i+ 1 2 ;j-1 2 ;k+ 1 2 (16∆x∆y ) + k12 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k12 i-1 2 ;j-1 2 ;k+ 1 2 (16∆x∆y ) -k12 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) + k11 i+ 1 2 ;j-1 2 ;k+ (16∆x∆x) + k11 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆x) + k11 i-1 2 ;j-1 2 ;k+ 1 2 (16∆x∆x) + k11 i-1 2 ;j-1 2 ;k- +( k13 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆z) + k12 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k11 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆x) -k23 i+ 1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k22 i+ 1 2 ;j-1 2 ;k-(16∆y ∆y ) + k21 i+ 1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k33 i+ 1 2 ;j-1 2 ;k-1 2 (16∆z∆z) -k32 i+ 1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) + k31 i+ 1 2 ;j-1 2 ;k- Abstract: Nowadays, heterogeneous materials are increasingly used for their superior overall properties, such as porous media, which are widely used in the electronics and biomedical industries, so determining the equivalent thermal conductivity (ETC) of heterogeneous materials is essential for the correct design of industrial equipment that may be subjected to severe thermal loads during use.

The main objective of this thesis is to calculate the homogenization of the thermal conductivity of heterogeneous materials using the finite difference method. Voxel was chosen for modeling heterogeneous materials and the Günter scheme will be employed as the primary technique for anisotropic thermal diffusion problems. The two-dimensional Günter system is re-demonstrated in this thesis, along with an extension to the three-dimensional model, as well as methods for loading periodic and mixed uniform boundary conditions. The three methods (FDM, FEM, and FEM+pixel(voxel)) are compared for 2D RVEs such as crosses, circles, and ellipses and for 3D RVEs such as spheres and cylinders. It is discovered that the developed FDM produces results that are consistent with those of FEM and FEM+pixel(voxel) and that the FDM outperforms FEM+pixel(voxel) in terms of convergence speed. This method has also been applied to sintered silver materials for the study of equivalent thermal conductivity.

Comparisons between the two methods (FDM and FEM) are carried out for the classical unit cells such as simple cubic, body-centered cubic, and face-centered cubic, as well as the silver-based stochastic model. The developed finite difference algorithm is valid, and consistent results are obtained. In addition to the Günter scheme, a 5-point model and an integral model have also been developed inspired by the Günter scheme.

For high-performance computing, the Eigen library and the Pardiso library are also detailed in the thesis. Both libraries contain both direct and iterative solutions for solving linear equations. However, while Eigen allows for parallel computation of only the iterative solution, Pardiso allows for parallel computation of both approaches, and the parallelism is significantly superior than that of Eigen. While Eigen is more straightforward to construct and more powerful, Pardiso is faster at tackling complex problems.

Maison du doctorat de l'Université Paris-Saclay 2 ème étage aile ouest, Ecole normale supérieure Paris-Saclay 4 avenue des Sciences, 91190 Gif sur Yvette, France

Finally, I would 3 . 7

 37 scheme Günter and 5-point for homogenized case . . . . . . . . . . . . . .

4. 11

 11 Comparison of path data of the anisotropic case of one cross inclusion of side length equal to 0.3 for thermal periodic boundary value problem. . . . 4.12 Convergence analysis of one disk inclusion of radius r= 0.3 . . . . . . . . . 4.13 Comparison of overall trends of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. . . . . . . 4.14 Comparison of path data of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. . . . . . . . . . 4.15 Convergence analysis with disk's radius equal to 0.3 with anisotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Comparison of overall trends of the anisotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. . . . . . . 4.17 Comparison of path data of the anisotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. . . . . . . 4.18 Convergence analysis with the long axis of the ellipse equal to 0.6 . . . . . 4.19 Comparison of overall trends of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. . . 4.20 Comparison of path data of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. . . . . 4.21 Convergence analysis with the long axis of the ellipse equal to 0.6 with anisotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Comparison of overall trends of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem.

4. 23

 23 Comparison of path data of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. . . 133 4.24 Convergence analysis of one disk and four sectors inclusions of radius r= 0.3 135 4.25 Comparison of overall trends of the isotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem.137 4.26 Comparison of path data of the isotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem.138 4.27 Convergence analysis of one disk and four sectors inclusions of radius r= 0.3 with anisotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . 139 4.28 Comparison of overall trends of the anisotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 4.29 Comparison of path data of the anisotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem.141 4.30 Comparison of overall trends of the isotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. . . . 144 4.31 Comparison of path data of the isotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. . . . . . . 145 4.32 Comparison of overall trends of the orthotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. . . 147 4.33 Comparison of path data of the orthotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. . . . 148 4.34 Convergence analysis with disk's radius equal to 0.3 with isotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.35 Comparison of overall trends of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . 4.36 Comparison of path data of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . . . . . 4.37 Convergence analysis with disk's radius equal to 0.3 with orthotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.38 Comparison of overall trends of the orthotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . 4.39 Comparison of path data of the orthotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . 4.40 Convergence analysis with the long axis of the ellipse equal to 0.6 with isotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . . . 4.41 Comparison of overall trends of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. . . 4.42 Comparison of path data of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. . . . . . 4.43 Convergence analysis with the long axis of the ellipse equal to 0.6 with orthotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . . . 4.44 Comparison of overall trends of the orthotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. .

4. 45

 45 Comparison of path data of the orthotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. . . 4.46 Convergence analysis of One disk and four sectors inclusions of radius r= 0.3 with isotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . 4.47 Comparison of overall trends of the isotropic case of One disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.48 Comparison of path data of the isotropic case of One disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. 4.49 Convergence analysis of one disk and four sectors inclusions of radius r= 0.3 with orthotropic thermal conductivity . . . . . . . . . . . . . . . . . . . . . 4.50 Comparison of overall trends of the orthotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.51 Comparison of path data of the orthotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. 4.52 Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with periodic boundary conditions. . . . 4.53 Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with MUBC boundary conditions. . . . 4.54 Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with different boundary conditions. . . . 4.55 3D layered composite model . . . . . . . . . . . . . . . . . . . . . . . . . . 4.72 Comparison of the equivalent thermal conductivity between the SC, BCC and FCC models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.73 Comparison of equivalent thermal conductivity between the non-conjugated and the conjugated models. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Equivalent thermal conductivity of an isotropic One disk and four sectors case with MUBC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.23 Coefficient of linear trend line for the orthotropic one disk and four sectors model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 Equivalent thermal conductivity of an orthotropic one disk and four sectors case with MUBC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.25 Equivalent thermal conductivity of 3D layered composite . . . . . . . . . . 4.26 Comparison of trendline parameters for Ball inclusion model with different R/L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.27 Comparison of trendline parameters for Cylindre inclusion model with different R/L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.28 Convergence results of different mesh levels by FED and FEM. . . . . . . . 4.29 Comparison of the equivalent thermal conductivity for random models with FDM and FEM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  1.1-(b)) (only cell edges exist in the material microstructure) and closed-cell structures (see Fig.1.1-(c)) (contain cell faces that separate each cell from its neighbours in the material microstructure). Fabrics are generally divided into woven bodies made by processes such as weaving, knitting, and sewing and non-woven bodies in which the fibers are mechanically, chemically or physically bonded into a possibly random web structure. Fabrics are generally divided into woven structures and non-woven structures. Woven structures are generally made of woven, knitted, stitched or knitted. Woven structures are inherently made of fibre and have a layered character, which can be divided into 2D woven (see Fig.1.1-(d)) and 3D woven structures (see Fig.1.1-(e))

  1.1-(f)) are defined by non-overlapping particles embedded in an interconnected matrix.

  (a) Olivine (b) Open-cell foam (c) Closed-cell foam (d) 2D woven fabric (From Texgen) (e) 3D woven fabric (From Texgen) (f) YAG(Y 3 AL 5 O 12

Figure 1 . 1 :

 11 Figure 1.1: Example of heterogeneous materials.

Fig. 1 .

 1 Fig.1.2-(b) and (c) show the model meshed by the Quadtree and Octree methods [48, 49].

Figure 1 . 2 :

 12 Figure 1.2: Mesh generation techniques

Figure 2 . 1 :

 21 Figure 2.1: Interpolation scheme

Figure 2 . 3 :

 23 Figure 2.3: Integral scheme for gradient operator

Figure 2 . 4 :

 24 Figure 2.4: Illustration of the 3D finite difference scheme based on voxel grid

  the point is on face, the volume fraction is ∆x∆y ∆z=(2 * (L x * L y * L z )). If the point is inside, the volume fraction is ∆x∆y ∆z=(L x * L y * L z ), where L x , L y and L z represent the model length, respectively.

. 5 ,

 5 PBCs are imposed between face ABCD and face A'B'C'D', between face ABB'A' and face DCC'D' and face ADD'A' and BCC'B'.

Figure 2 . 5 :

 25 Figure 2.5: RVE with periodic boundary conditions

T

  Edg e-AA -T Edg e-DD = ∆T z T Edg e-CC -T Edg e-DD = ∆T x T Edg e-BB -T Edg e-CC = ∆T z (2.82) Edg e-AB -T Edg e-A B = ∆T y T Edg e-AB -T Edg e-DC = ∆T z T Edg e-DC -T Edg e-D C = ∆T y (2.83) Edge-BC, Edge-B'C', Edge-A'D' and Edge-AD:

5 ,

 5 Dirichlet boundary conditions are imposed between face ADD'A' and face B CC'BB'. Neumann boundary conditions are imposed between face ADCB and face A'D'C'B' and between face ABB'A' and face DCC'D'.

  96) Edge-AD, Edge-DD', Edge-D'A', Edge-AA': T Edg e-AD = T Edg e-DD = T Edg e-D A = T Edg e-AA = T lef t (2.97) Edge-BC, Edge-CC', Edge-B'C', Edge-BB': T Edg e-BC = T Edg e-CC = T Edg e-B C = T Edg e-BB = T r ig ht (2.98) 71 Edge-AB, Edge-DC, Edge-D'C' and Edge-A'B':

3. 2 2 Figure 3

 223 Figure 3.1: 5-point scheme

Figure 3 .

 3 Figure 3.2: integral model

  3.3-(a) as an example, we load Dirichlet boundary conditions on the left and right boundaries, i.e. the left and right boundary values are T lef t and T r ig ht . We can obtain the following:

Figure 3 . 3 :

 33 Figure 3.3: MUBC boundary conditions

  3.3-(a) shows the 4*4 pixel schematic diagram of MUBC boundary conditions. The numbering of 0-15 in the figure is the point numbering, and 0'-15' are the virtual points outside the boundary. We take the Günter isotropic scheme as an example of loading fixed temperatures (T r ig ht ; T lef t ) on the left and right sides and insulation (Φ q = 0) on the top and bottom sides. When loading the MUBC boundary condition, we need to load both Dirichlet and Neumann boundary condition on the boundary. As shown in Fig.3.3-(a), the square is loaded with the equilibrium equation introduced by Table

Fig. 3 .

 3 Fig.3.4-(a) shows the schematic diagram of periodic boundary conditions for a 4*4 pixel

Figure 3 . 4 :

 34 Figure 3.4: Periodic boundary conditions

Figure 3 . 5 :

 35 Figure 3.5: PBC on 4x4x4 voxel

(3. 43 )Figure 3 . 6 : on 4x4x4 voxel 95 The

 433695 Figure 3.6: on 4x4x4 voxel

Figure 3 . 7 :

 37 Figure 3.7: Coding process

5 -

 5 point model and the integral model. The 5-point model does not have the problem of loading boundary conditions and is simple to implement. On the other hand, the integral model can be obtained by changing the parameters to obtain various models (including the Günter model). Secondly, we propose methods for loading MUBC and periodic boundary conditions, where we present several methods for loading Neumann boundary conditions, and also solve the problem of matrix singularities caused by the Günter model when loading Neumann boundary conditions using the virtual point method.

  + D y y ) 2 sin(ıx)sin(ıy ) (4.1) where D is diffusion tensor, S is source terme and ' is varied from 0 to 9. T exact (x; y ) = 1 2ı 2 sin(ıx)sin(ıy ); x; y ∈ [0; 1]; (4.2) m e Günte r s c he m e 5 -po int Line a r (Erro r s c he m e Günte r) Line a r (Erro r s c he m e 5 -po int ) e m e Gü n te r s c h e m e 5 -p o in t (c) Accuracy for different degree of anisotropy

Figure 4 . 1 :

 41 Figure 4.1: Solution contours by scheme Günter and order of accuracy obtained by scheme Günter and 5-point for homogenized case

Fig4. 2 -

 2 Fig4.2-(a) shows the solution contours by scheme Günter with 50*50 pixels. The order of

Figure 4 . 2 :

 42 Figure 4.2: Solution contours by scheme Günter and order of accuracy obtained by scheme Günter and 5-point for heterogeneous case with orthotropic tensor coefficient

. 7 )

 7 Fig4.3-(a) shows the solution contours by scheme Günter with 80*80 pixels. The order of

Figure 4 . 3 :

 43 Figure 4.3: Solution contours by scheme Günter and order of accuracy obtained by scheme Günter and 5-point for heterogeneous case with anisotropic tensor coefficient

11 )

 11 Fig4.4-(a) shows the solution contours by scheme Günter with 80*80 pixels. The order

Figure 4 . 4 :

 44 Figure 4.4: Solution contours by scheme Günter and order of accuracy obtained by scheme Günter and 5-point for heterogeneous case II with anisotropic tensor coefficient

Fig4. 4 -

 4 Fig4.4-(c) shows that the error remains invariant as the D xx value increases, with D xx

  Figure 4.5: 1D model

Fig. 4 . 6 -

 46 Fig.4.6-(a) shows a slat consisting of a heterogeneous material in the format of 3*1 pixels,

(4. 15 )Fig. 4 . 6 -

 1546 Fig.4.6-(b) and (c) show the obtained temperature fields and temperature gradients, respectively. The final equivalent thermal conductivity can be obtained with the 3*1 pixel as follows:

Figure 4 .

 4 Figure 4.6: 3*1 pixel model of Slat

Figure 4 . 7 :Fig. 4 .Figure 4 . 8 :

 47448 Figure 4.7: Simple geometric shapes

Figure 4 . 9 :

 49 Figure 4.9: Comparison of path data of the isotropic case of one cross inclusion of side length equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux

Fig. 4 . 116 Temperature

 4116 Fig.4.10 shows the comparison of overall trends of the anisotropic case of one cross inclusion of side length equal to 0.3 for the thermal periodic boundary value problem. The structure contains 2025 pixel cells for the FDM and 2025 quadrilateral finite elements for the FEM and FEM+pixels. By comparing the temperature field and the temperature gradient field (see (a)-(f)), we find that the distributions obtained by the three methods are quite consistent. If we compare the distribution of the heat flux field (see (g)-(i)), we can see that the general trend is consistent, except that the values of the intersections obtained by FEM and FEM+pixel are much smaller than those of FDM in the model, but the difference in the inflection points is changed from four to two compared to the isotropic model. The comparison of path data of the anisotropic case of one cross inclusion of side length equal to 0.3 for thermal periodic boundary value problem is shown in Fig.4.11. The comparison shows that the data on all paths match well, including the heat flow comparison, which

Figure 4 . 10 :

 410 Figure 4.10: Comparison of overall trends of the anisotropic case of one cross inclusion of side length equal to 0.3 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 11 :

 411 Figure 4.11: Comparison of path data of the anisotropic case of one cross inclusion of side length equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux

Fig. 4 .

 4 Fig.4.12-(a) shows the trend of the equivalent thermal conductivity as the DOF increases (denser grid), and Fig.4.12-(b) shows the error variation of both methods, FDM andFEM+pixel, with the FEM value as the reference, as well as the geometric error. The comparison shows that superior to the pixel grid characteristics, the geometric error is more significant for the FDM and FEM+pixel methods when the DOF is less. Secondly, in Fig.4.12-(b), the trend line expression for FDM is y = -0:54366665 * x -0:0731229, whereas, for FEM+pixel, the trend line expression is y = -0:54472008 * x + 0:04362687.

Fig. 4 .

 4 Fig.4.13 shows the comparison of overall trends of the isotropic case of one disk inclusion of radius equal to 0.3. The FDM and FEM+pixel contain 2025 cells, and 2257 cells

Figure 4 . 13 :

 413 Figure 4.13: Comparison of overall trends of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 14 :

 414 Figure 4.14: Comparison of path data of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

1 .

 1 nt the rm a l c o nduc tivity (W.m -Erro r FEM with Pixe l Line a r (Erro r FDM) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) (b) Error comparison

Figure 4 . 15 :Fig. 4 .

 4154 Figure 4.15: Convergence analysis with disk's radius equal to 0.3 with anisotropic thermal conductivity

Fig. 4 .

 4 Fig.4.15-(a) shows the trend of the equivalent thermal conductivity as the DOF increases (denser grid), and Fig.4.15-(b) shows the variation in error for both methods, FDM and 1283 using the FEM values as a reference and also showing the geometric error, where the expression for the error trend line for FDM is y = -0:55026384 * x -0:04117075 and the

Fig. 4 .

 4 Fig.4.16 shows the comparison of the temperature field, temperature gradient field, and heat flux field for the thermal periodic boundary value problem. The FDM and FEM+pixel contain 4225 cells, and 4724 cells are wrapped in the FEM. In the comparison of the temperature field (see (a), (b), and (c)), the images obtained for the three methods all show the same trend and are clearly different from the isotropic model. (d)-(i)illustrate the change in temperature gradient and heat flux. We find that the trends change in the same way for all three methods, with the differences mainly concentrated in the region between the inclusion and matrix phases. The comparison of path data of the anisotropic case of one disk inclusion of radius r equal to 0.3 for thermal periodic boundary value problem is shown in Fig.4.17. We choose paths with y = 0:5 and y = x respectively to compare, where (a) and (b) indicate the temperature change over the path, and it can be seen that the results of the three methods agree well.

Figure 4 . 16 :

 416 Figure 4.16: Comparison of overall trends of the anisotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 17 :

 417 Figure 4.17: Comparison of path data of the anisotropic case of one disk inclusion of radius equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

1 .

 1 Eq uiva le nt the rm a l c o nduc tivity (W.m -Erro r FEM with Pixe l Line a r (Erro r FDM) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) (b) Error comparison

Figure 4 . 18 : 6 Fig. 4 . 127 Temperature

 41864127 Figure 4.18: Convergence analysis with the long axis of the ellipse equal to 0.6

Figure 4 . 19 :

 419 Figure 4.19: Comparison of overall trends of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 20 :

 420 Figure 4.20: Comparison of path data of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Fig. 4 .

 4 Fig.4.19 shows the comparison of overall trends of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for the thermal periodic boundary value problem. The FDM and FEM+pixel contain 2025 cells, and 2015 cells are wrapped in the FEM. The temperature changes are shown in (a)-(c), all three methods produce the same overall full trend change. The variation of the temperature gradient and heat flux is shown in (d)-(i), respectively. They show that the overall trend variation is the same for all three methods, but the discontinuity in thermal conductivity causes errors at the junction between the inclusion and matrix phases. The comparison of path data of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem is shown in Fig.4.20. We have selected paths with y = 0:5 and y = x respectively for comparison. The comparison of temperatures on the paths shown in (a) and (b) illustrates the consistency of the three methods in terms of temperature results. In the comparison of the paths of the temperature gradients (see (c) and (d)), it can be seen that the values of FDM at the interface between the inclusion and matrix phases are between those of FEM and FEM+pixel. FDM

1 .

 1 nt the rm a l c o nduc tivity (W.m -Erro r) Erro r Ge o m e try Erro r FDM Erro r FEM with Pixe l Line a r (Erro r FDM) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) (b) Error comparison

Figure 4 . 21 :Fig. 4 .

 4214 Figure 4.21: Convergence analysis with the long axis of the ellipse equal to 0.6 with anisotropic thermal conductivity

Figure 4 . 22 :

 422 Figure 4.22: Comparison of overall trends of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 23 :

 423 Figure 4.23: Comparison of path data of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Fig. 4 .Fig. 4 .

 44 Fig.4.22 shows the comparison of overall trends of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for the thermal periodic boundary value problem. The FDM and FEM+pixel contain 9025 cells, and 9265 cells are wrapped in the FEM.Fig.4.22-(a)-(c) show the temperature variation, the overall trend is clearly different from the isotopic variation, and secondly, there is good agreement between the results of the three methods. The change in temperature gradient and heat flux is shown in (d)-(i), comparing the previous isotropic model. The results of the three methods agree better when the mesh density is increased. However, some errors still occur at the junction between the inclusion and matrix phases. The comparison of path data of the anisotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal periodic boundary value problem is shown in Fig.4.23. Comparingthe temperatures along the paths demonstrates that the three approaches accord pretty well. When the temperature gradients along the paths are compared, it is clear that all locations match precisely, except those between the inclusion and matrix phases. When examining the heat flow throughout the path, it can be seen that as the mesh density increases, the spots within the inclusion phase of the FDM and FEM+pixel model more closely resemble the FEM. However, the FEM+pixel generates far more significant inaccuracies at the inflection points.

Fig. 4 .

 4 Fig.4.25-(a)-(c) illustrate the changes in the temperature field, and it is clear that the overall trend is the same for all three methods. (d)-(i) illustrate the variation of temperature gradient and heat flux, respectively. The comparison demonstrates that the overall trend isthe same for all three methods, but there are differences between the inclusion and matrix phases.

Fig. 4 .

 4 Fig.4.26 illustrates the comparison of path data of the isotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem, and the comparison is made between the paths y = 0:5 and y = x. The temperatures on the paths shown in (a) and (b) demonstrate the consistency of the three methods in solving for temperature. In the comparison of temperature gradients on the path (see (c) and (d)), wefind that FDM is slightly better than FEM+pixel. Even though there are some differences between the three at the point between the inclusion and matrix phases, FEM+pixel has a significant error at the inflection point. (e) and (f) illustrate the comparison of heat flow on the paths, and it can be seen that the value of FDM remains between FEM and FEM+pixel in the different sections of the three methods, indicating that FDM is superior to FEM+pixel and that FEM+pixel has a more significant difference value at the inflection point.

Figure 4 . 25 :

 425 Figure 4.25: Comparison of overall trends of the isotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 26 :Fig. 4 .

 4264 Figure 4.26: Comparison of path data of the isotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

1 .

 1 nt the rm a l c o nduc tivity (W.m -Ge o m e try Erro r FDM Erro r FEM with Pixe l Line a r (Erro r FDM) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) (b) Error comparison

Figure 4 . 27 :

 427 Figure 4.27: Convergence analysis of one disk and four sectors inclusions of radius r= 0.3 with anisotropic thermal conductivity

Figure 4 . 28 :

 428 Figure 4.28: Comparison of overall trends of the anisotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 29 :Fig. 4 .

 4294 Figure 4.29: Comparison of path data of the anisotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal periodic boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

  The comparison of the overall trend of a cross-inclusion with side length equal to 0.3 for the thermal MUBC boundary value problem in the isotropic case is shown in Fig.4.30, where (a)-(c) represents the temperature field variation, (d)-(f) represents the temperature gradient field variation, and (g)-(i) represents the heat flow field variation (all three methods contain 2025 cells). To begin, in the scenario where the grid border corresponds with the inclusive phase boundary, the findings for FEM and FEM+pixel are identical. Second, it

Fig. 4 .Figure 4 . 30 :

 4430 Fig.4.312 for the thermal MUBC boundary value problem. We have compared primarily two pathways at y = 0.5 and y = x. It can be seen that the path comparison yields the same conclusions as Fig.4.30, where the temperature results on the paths of the three methods perfectly overlap, and where the temperature gradient and heat flow results on the paths perfectly overlap at all points, except for the values of the cross model's inflection points.

Figure 4 . 31 :

 431 Figure 4.31: Comparison of path data of the isotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Fig. 4 .

 4 Fig.4.32 shows the comparison of overall trends of the orthotropic case of one cross inclusion of side length equal to 0.3 for the thermal MUBC boundary value problem. Fig.4.33 shows the comparison of path data of the orthotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. All three of these methods contain 5625 cells.

Figure 4 . 32 :

 432 Figure 4.32: Comparison of overall trends of the orthotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

  ra ture FDM pa th dia g o na l FEM+pixe l pa th dia g o na l FEM pa th dia g o na l pa th dia g o na lFEM+pixe l pa th dia g o na l FEM pa th dia g o na l dia g o na l FEM+pixe l pa th dia g o na l FEM pa th dia g o na l (c) (d)

Figure 4 . 33 : 1 .

 4331 Figure 4.33: Comparison of path data of the orthotropic case of one cross inclusion of side length equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Figure 4 . 34 :

 434 Figure 4.34: Convergence analysis with disk's radius equal to 0.3 with isotropic thermal conductivity

Figure 4 . 35 :

 435 Figure 4.35: Comparison of overall trends of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 36 :

 436 Figure 4.36: Comparison of path data of the isotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

  5 and y = x to perform the comparison displayed in Fig.4.36. The comparison of temperatures on the paths (a) and (b) indicates the three techniques' agreement in terms of temperature values. The comparison of temperature gradient (see (c) and (d)) and heat flow path data (see (e) and (f)) reveals that the primary differences with FEM are concentrated in the region between the inclusion and matrix phases. However, the overall comparison is better for FDM than FEM+pixel, primarily because FEM+pixel has a tremendous difference value at the inflection point, whereas FDM has a much smaller difference value.

  3 with orthotropic thermal conductivity input for the mixed uniform boundary conditions (MUBC) is shown in Fig.4.37. As DOF increases, both FDM and FEM+pixel converge gradually towards the FEM value, and FDM is always between FEM and FEM+pixel. Fig 4.37-b illustrates a comparison of the errors using FEM values as a reference, and by comparing the trend line coefficients in Table

1 .

 1 nt the rm a l c o nduc tivity (W.m -Erro r) Erro r Ge o m e try Erro r FDM NNPG3 Erro r FEM with Pixe l Erro r FDM NNPG2 Erro r FDM NNPG1 Erro r FDM NA Erro r FDM NWPG Line a r (Erro r FDM NNPG3 ) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) Line a r (Erro r FDM NWPG) Line a r (Erro r FDM NA) Line a r (Erro r FDM NNPG1 ) Line a r (Erro r FDM NNPG2 ) (b) Error comparison

Figure 4 . 37 :

 437 Figure 4.37: Convergence analysis with disk's radius equal to 0.3 with orthotropic thermal conductivity

Figure 4 . 38 :

 438 Figure 4.38: Comparison of overall trends of the orthotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 39 :

 439 Figure 4.39: Comparison of path data of the orthotropic case of one disk inclusion of radius equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Fig. 4 .

 4 Fig.4.39 illustrates a comparison of the path data (y = 0:5 and y = x) for the pictures in Fig.4.38, demonstrating that the three approaches agree well when comparing the temper-

4 . 40 .

 440 By increasing the DOF, the FDM and FEM+pixel gradually converge towards the FEM (see Fig.4.40-a). The comparison of errors using the FEM values as a reference is shown in Fig.4.40-b. According to the trend line coefficients in Table4.40, the FDM converges much quicker than the FEM+pixel under the isotropic elliptical model and at a slightly different rate than the disc model.

Fig. 4 . 1 .

 41 Fig.4.41 shows the comparison of the temperature field, temperature gradient field, and heat flux field for the thermal periodic boundary value problem. The FDM and FEM+pixel

Figure 4 . 40 :

 440 Figure 4.40: Convergence analysis with the long axis of the ellipse equal to 0.6 with isotropic thermal conductivity

Fig. 4 .

 4 Fig.4.42 compares the path data for the temperature field, temperature gradient field, and heat flow field in Fig.4.41. comparison of temperatures along the lines illustrated in (a) and (b) demonstrates the consistency of the three approaches for calculating temperatures.

Figure 4 . 41 :

 441 Figure 4.41: Comparison of overall trends of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 42 :

 442 Figure 4.42: Comparison of path data of the isotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

1 .

 1 nt the rm a l c o nduc tivity (W.mo m e try Erro r FDM NNPG3 Erro r FEM with Pixe l Erro r FDM NNPG2 Erro r FDM NNPG1 Erro r FDM NA Erro r FDM NWPG Line a r (Erro r FDM NNPG3 ) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) Line a r (Erro r FDM NWPG) Line a r (Erro r FDM NA) Line a r (Erro r FDM NNPG1 ) Line a r (Erro r FDM NNPG2 ) (b) Error comparison

Figure 4 . 43 :

 443 Figure 4.43: Convergence analysis with the long axis of the ellipse equal to 0.6 with orthotropic thermal conductivity

  6 with orthotropic thermal conductivity is shown in Fig.4.43. As the DOF rises, both FDM and FEM+pixel gradually converge towards FEM, but the value of FDM is always between FEM and FEM+pixel (see Fig.4.43-a). Fig.4.43-b compares the FDM and FEM+pixel errors using the FEM values as a reference. Comparing the trend line parameters in Table 4.19 reveals that for the elliptical model, alternative loading Newman techniques and input thermal conductivities (isotropic or anisotropic) have no effect on the rate of convergence. However, it is evident that the FDM has a faster convergence rate than the FEM+pixel and that the FDM advantage is more significant for the elliptical model than for the disc model. The variations in the temperature field, temperature gradient field, and heat flow field for an orthotropic ellipse model are shown in Fig.4.44. The FDM and FEM+pixel have 2500 cells each, whilst the FEM has 2428 cells.

Fig. 4 .

 4 Fig.4.45 illustrates the comparison of the path data in Fig.4.44, demonstrating that the three approaches agree well when comparing the temperatures on the paths in (a) and (b). The comparison of temperature gradients and heat flow in (c)-(f) demonstrates that

Figure 4 . 44 :

 444 Figure 4.44: Comparison of overall trends of the orthotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

  ra ture FDM pa th dia g o na l FEM+pixe l pa th dia g o na l FEM pa th dia g o na l dia g o na lFEM+pixe l pa th dia g o na l FEM pa th dia g o na l dia g o na l FEM+pixe l pa th dia g o na l FEM pa th dia g o na l

Figure 4 . 45 :

 445 Figure 4.45: Comparison of path data of the orthotropic case of one ellipse inclusion of the long axis equal to 0.6 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

Figure 4 . 47 :

 447 Figure 4.47: Comparison of overall trends of the isotropic case of One disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 48 :

 448 Figure 4.48: Comparison of path data of the isotropic case of One disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

  4.48), we used data from paths with y = 0:5 and y = x. Whereas (a) and (b) illustrate temperature comparisons along the pathways, it is clear that the three methodologies agree pretty well. The comparison of temperature gradients in (c) and (d) demonstrates a distinction between the inclusion phase and the matrix phase. The heat flow comparisons in (e) and (f) demonstrate that the FDM values are substantially closer to the FEM values than the FEM+pixel values. The FEM+pixel values appear to have a more significant difference value at some inflection points.

  3 with orthotropic thermal conductivity is shown in Fig.4.49. where (a) depicts the change in equivalent thermal conductivity, i.e., the FDM and FEM+ pixel values gradually converge to the FEM value as the DOF increases, the FDM value is always between the FEM and FEM+ pixel values, and the curves derived using the various loading Neumann methods perfectly overlap. (b) compares the FDM and FEM+ pixel errors using the FEM values as a reference, and the trend line's coefficients are presented in Table

1 .

 1 Eq uiva le nt the rm a l c o nduc tivity (W.mo m e try Erro r FDM NNPG3 Erro r FEM with Pixe l Erro r FDM NNPG2 Erro r FDM NNPG1 Erro r FDM NA Erro r FDM NWPG Line a r (Erro r FDM NNPG3 ) Line a r (Erro r FEM with Pixe l) Line a r (Erro r Ge o m e try) Line a r (Erro r FDM NWPG) Line a r (Erro r FDM NA) Line a r (Erro r FDM NNPG1 ) Line a r (Erro r FDM NNPG2 ) (b) Error comparison

Figure 4 . 49 :

 449 Figure 4.49: Convergence analysis of one disk and four sectors inclusions of radius r= 0.3 with orthotropic thermal conductivity

Fig. 4 . 168 Temperature

 4168 Fig.4.50 show the variation of the temperature field, temperature gradient field, and heat flow field in the orthotropic case of one disk and four sectors inclusions of radius r

Figure 4 . 50 :

 450 Figure 4.50: Comparison of overall trends of the orthotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature gradient field computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Heat flux field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel.

Figure 4 . 51 :

 451 Figure 4.51: Comparison of path data of the orthotropic case of one disk and four sectors inclusions of radius equal to 0.3 for thermal MUBC boundary value problem. path middle(y = 0:5): (a) Temperature, (c) Temperature gradient, (e) Heat flux. path diagonal(y = x): (b) Temperature, (d) Temperature gradient, (f) Heat flux.

  4.51). The comparison of (a) and (b) demonstrates the uniformity of the three approaches for determining the temperatures. Second, when the path data for the temperature gradient and heat flow are compared in (c)-(e), the value of FEM+pixel, while gradually approaching the FEM and FDM values as the number of cells increases, nevertheless results in a high difference value at the inflection point.

  Anisotropic Disk and 4 secteurs

Figure 4 . 52 :

 452 Figure 4.52: Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with periodic boundary conditions.

  Orthotropic Disk and 4 secteurs

Figure 4 . 53 :

 453 Figure 4.53: Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with MUBC boundary conditions.

Fig. 4 . 1 . 1 . 1 .

 4111 Fig.4.52 and Fig.4.53 show the variation of equivalent thermal conductivity with volume fraction for various 2D models under periodic boundary conditions and MUBC boundary conditions, respectively. The periodic boundary condition includes isotropic and anisotropic input thermal conductivity cases, while the MUBC boundary condition includes isotropic and orthotropic input thermal conductivity cases. It can be seen that the equivalent thermal conductivity obtained by all three methods after convergence agrees well for either boundary condition under these displayed models.

Figure 4 . 54 :

 454 Figure 4.54: Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various 2D models with different boundary conditions.

Fig. 4 .

 4 Fig.4.55-(a) shows the 3D Layered composite model, with the dimensions of inclusion equal to L/3. We can use Voigt's and Reuss's values to find the analytical solution of the ETC for

Figure 4 . 55 :

 455 Figure 4.55: 3D layered composite model

Fig. 4 .

 4 Fig.4.56 shows the convergence analysis of the ball inclusion model, which includes the variation of ETC with increasing DOF, and the error analysis at different R/L using the FEM values as reference. As can be seen in the variation of the equivalent thermal conductivity shown in (a), (c), (e), (g), (i), (k), and (m), the FDM and FEM+voxel gradually converge towards the FEM value as the DOF increases. However, the FDM value remains closer to the FEM and FEM+voxel values.(b), (d), (f), (h), (j), (l) and (n) show the error analysis for different radii with FEM as a reference, respectively, and it can be seen through the images that the errors for FDM are all smaller than those for FEM+voxel. The trend line parameters for FDM and FEM+voxel are shown in Table4.[START_REF] Sharma | Size-dependent eshelby's tensor for embedded nanoinclusions incorporating surface/interface energies[END_REF] The comparison shows that the rate of convergence

Fig. 4 .

 4 Fig.4.57 shows an overall comparison of the temperature field, temperature gradient field, and heat flow field for the ball inclusion model as well as a comparison of the crosssectional data. The structure contains 8000 cubic voxel cells by FDM and FEM+voxel, and 60374 tetrahedral finite elements by FEM. The data structure for FEM and FEM+voxel

Figure 4 . 57 :

 457 Figure 4.57: Ball inclusion model with radius/L equal to 0.25 for thermal MUBC boundary value problem. Temperature field computed by: (a) FDM, (b) FEM, (c) FEM+Pixel. Temperature field in cross section xOz computed by: (d) FDM, (e) FEM, (f) FEM+Pixel. Temperature gradient field computed by:(g) FDM, (h) FEM, (i) FEM+Pixel. Temperature gradient field in cross section xOz computed by:(j) FDM, (k) FEM, (l) FEM+Pixel. Heat flux field computed by:(m) FDM, (n) FEM, (o) FEM+Pixel. Heat flux field in cross section xOz computed by:(m) FDM, (n) FEM, (o) FEM+Pixel.

Fig. 4 . 1 .

 41 Fig.4.58 shows the variation of ETC with a different fraction of the ball inclusion model.

Figure 4 .

 4 Figure 4.58: ETC with different fraction of ball inclusion model.

Fig. 4 .

 4 Fig.4.60 and Fig.4.61 show an overall comparison of the temperature field, temperature gradient field, and heat flow field for the spherical inclusion model, as well as a comparison of the cross-sectional data. The structure contains 8000 cubic voxel cells by FDM and FEM+voxel, and 59889 tetrahedral finite elements by FEM. The data structure for FEM and FEM+voxel is taken from Ansys, where the temperature values are taken from the nodes and the temperature gradient and heat flow values are taken from the cell centroids.

Figure 4 . 60 :

 460 Figure 4.60: Cylinder model with radius/L equal to 0.25 for thermal MUBC boundary value problem with Dirichlet en X-direction. Temperature field computed by: (a) FDM, (b) FEM, (c) FEM+Pixel. Temperature field in cross section xOz computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Temperature gradient field computed by: (g) FDM, (h) FEM, (i) FEM+Pixel. Temperature gradient field in cross section xOz computed by:(j) FDM, (k) FEM, (l) FEM+Pixel. Heat flux field computed by:(m) FDM, (n) FEM, (o) FEM+Pixel. Heat flux field in cross section xOz computed by:(m) FDM, (n) FEM, (o) FEM+Pixel.

Figure 4 . 61 :

 461 Figure 4.61: Cylinder model with radius/L equal to 0.25 for thermal MUBC boundary value problem with Dirichlet en Z-direction. Temperature field computed by:(a) FDM, (b) FEM, (c) FEM+Pixel. Temperature field in cross section xOz computed by:(d) FDM, (e) FEM, (f) FEM+Pixel. Temperature gradient field computed by: (g) FDM, (h) FEM, (i) FEM+Pixel. Temperature gradient field in cross section xOz computed by: (j) FDM, (k) FEM, (l) FEM+Pixel. Heat flux field computed by:(m) FDM, (n) FEM, (o) FEM+Pixel. Heat flux field in cross section xOz computed by:(m) FDM, (n) FEM, (o) FEM+Pixel.

Fig. 4 . 1 .

 41 62 shows the variation of ETC with a differentEq uiva le nt the rm a l c o nduc tivity (W.m -

Figure 4 .

 4 Figure 4.62: ETC with different fraction of cylinder inclusion model.

4. 8

 8 Application to REV of complex heterogeneous materials 4.8.1 sintered silver

Figure 4 . 63 : 1 Figure 4 . 64 :

 4631464 Figure 4.63: Illustration of geometrical models considered: (a)-(c) structures of simple cubic (SC), body-centered cubic (BCC) and face-centered cubic (FCC), the colored parts indicate the overlaps; (d)-(f) the corresponding periodical unit cells; (g)-(i) the corresponding conjugated structure models.

  To proceed with the comparison of local fields, several geometrical models are selected as examples in this section. The first one is the SC model with the volume fraction of silver equal to 67:18%, which is shown in Fig.4.65 (corresponding to a sphere radius equal to 0:55). The overall temperature field (see Fig.4.65a for FDM, Fig.4.65b for FEM) and heat flux field (see Fig.4.65c for FDM, Fig.4.65d for FEM) as well as their middle cross sections (see Fig.4.65e-h) are displayed. The unit structure is divided into 45 × 45 × 45 = 91125 cubic voxel cells by FDM, while it is divided into 31460 tetrahedral finite elements by FEM. We can see from the figure that the distributions of the temperature and the heat flux fields by the two methods are comparable and consistent, which demonstrates again the validity of the FD system we developed in this study. The results for the SC conjugated model are shown in Fig.4.66. The fraction of silver is 32:82%. The numbers of voxel cells and that of finite elements are the same as in the non-conjugated model. It can be seen that the temperature and the heat flux fields in this situation by the two methods are also matched and coordinated, which shows the effectiveness of the FD algorithm in the conjugated model.Temperature (°C) Heat flux (W•m -2 )

Figure 4 . 65 :

 465 Figure 4.65: SC model with the volume fraction of silver equal to 67:18%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOy by FDM, (f) temperature field in cross section xOy by FEM, (g) heat flux field in cross section xOy by FDM, (h) heat flux field in cross section xOy by FEM.

Figure 4 .

 4 Figure 4.66: SC conjugated model with the volume fraction of silver equal to 32:82%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOy by FDM, (f) temperature field in cross section xOy by FEM, (g) heat flux field in cross section xOy by FDM, (h) heat flux field in cross section xOy by FEM.

Figure 4 . 68 :Figure 4 . 67 :

 468467 Figure 4.68: BCC conjugated model with the volume fraction of silver equal to 20:11%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOy by FDM, (f) temperature field in cross section xOy by FEM, (g) heat flux field in cross section xOz by FDM, (h) heat flux field in cross section xOz by FEM.

Figure 4 . 69 :Figure 4 . 70 :

 469470 Figure 4.69: FCC model with the volume fraction of silver equal to 88:07%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOy by FDM, (f) temperature field in cross section xOy by FEM, (g) heat flux field in cross section xOy by FDM, (h) heat flux field in cross section xOy by FEM.

Figure 4 . 71 :

 471 Figure 4.71: Evolution of the equivalent thermal conductivity as a function of the volume fraction for the various models. The star symbols represent the values obtained by FEM. The cross symbols, the circle markers, the triangle-down markers correspond to the values of FDM with 45 × 45 × 45, 85 × 85 × 85 and 125 × 125 × 125 voxels, respectively. The tree markers designate the values calculated by theoretical formula of Maxwell.

Fig. 4 .

 4 Fig.4.72 compares the equivalent thermal conductivity of SC, BCC and FCC models.

Figure 4 . 72 :

 472 Figure 4.72: Comparison of the equivalent thermal conductivity between the SC, BCC and FCC models.

  nt the rm a l c o nduc tivity(w/(m .k)) nt the rm a l c o nduc tivity(w/(m .k)) nt the rm a l c o nduc tivity(w/(m .k)) FCC and FCC conjugated models

Figure 4 . 73 :

 473 Figure 4.73: Comparison of equivalent thermal conductivity between the non-conjugated and the conjugated models.

Figure 4 . 74 :

 474 Figure 4.74: Geometrical illustration of three random models. (a) case 1: volume fraction of silver = 66.3%, (b) case 2: volume fraction of silver = 73.4% , (c) case 3: conjugated structure with volume fraction of silver = 73.4% (i.e. the silver is for the matrix, and the air is for the sphere inclusion).

Fig. 4 .

 4 Fig.4.75 shows the random model with the fraction of silver equal to 73:4% (corresponding to an average radius of 0:146 and a standard deviation of 0:028). The structure contains 614125 cubic voxel cells by FDM, and 473920 tetrahedral finite elements by FEM. By comparing the temperature fields (see Fig.4.75a,b,e,f), we find that the distributions obtained by the two methods are quite consistent. If comparing the distribution of the heat flux field (see Fig.4.75c,d,g,h), it is shown that the overall trend is the same, but in some

Fig. 4 .

 4 Fig.4.[START_REF] Dirrenberger | Effective elastic properties of auxetic microstructures: Anisotropy and structural applications[END_REF] shows the conjugated random model with the fraction of silver equal to 73:4%

Figure 4 . 75 :

 475 Figure 4.75: Stochastic model with the volume fraction of silver equal to 73:4%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOz by FDM, (f) temperature field in cross section xOz by FEM, (g) heat flux field in cross section y Oz by FDM, (h) heat flux field in cross section y Oz by FEM.

  flux (W•m -2 ) -610-500 -400 -300 -200 -100 -0.027

Figure 4 . 76 :

 476 Figure 4.76: Stochastic conjugated model with the volume fraction of silver equal to 73:4%: (a) temperature field by FDM, (b) temperature field by FEM, (c) heat flux field by FDM, (d) heat flux field by FEM, (e) temperature field in cross section xOz by FDM, (f) temperature field in cross section xOz by FEM, (g) heat flux field in cross section y Oz by FDM, (h) heat flux field in cross section y Oz by FEM.

Figure 4 . 77 :

 477 Figure 4.77: Comparison between Eigen solvers

Figure 4 .

 4 Figure 4.78: comparison different solvers

  Dc solver with different threads

Figure 4 . 79 :

 479 Figure 4.79: Comparison between Pardiso solvers

Figure 4 . 80 :

 480 Figure 4.80: Dc and cg solver comparison under 4 and 6 threads

Figure 4 . 81 :

 481 Figure 4.81: Comparison Eigen and Pardiso under 4threads

(

  -n+1 --n )exp(ix j k) = fl x -n exp(ix j k)[exp(-i∆xk) -2 + exp(i∆xk)] --1 = fl x [exp(-i∆xk) -2 + exp(i∆xk)] -= 1 -4fl x sin 2 (∆xk=2) -| < 1 → fl x < 1=2 (5.3)In the same way, stability condition 2D and 3D: fl x + fl y < 1=2, fl x + fl y + fl z < 1=2.

Fig. 5 .Figure 5 . 1 :-- 1 =

 5511 Fig.5.1 shows the temperature variation results of the disk model obtained by the Euler explicit scheme. The structure contains 900 pixel cells, and when the step size is set to

Figure 5 . 2 :

 52 Figure 5.2: Disk model with Euler implicit scheme

Figure 5 . 3 :

 53 Figure 5.3: Disk model with Crank-Nicolson scheme

(5. 8 )

 8 This algorithm divides the time step from n to n+1 into two sub-time stages: n to n+1/2 and n+1/2 to n+1. In the first step, use the implicit format to solve in the xdirection, and make the y-direction use the explicit format to solve. In the second step, similarly use the x directions to solve the problem with the implicit form.

Fig. 5 .Figure 5 . 4 :

 554 Fig.5.4 shows the temperature variation results of the disk model obtained by the Peaceman-Rachford method. The structure contains 900 pixel cells, and when the step size is set to 0.002, 412 iterations are required to achieve a tolerance of 1e-6.

(5. 10 )

 10 Compared with the Crank-Nicolson method, the ADI method significantly saves computing time and computing space. Consider a 2-D mesh with sizeN=(nx+1)(ny+1). In the first step, there are ny+1 values, and each value takes O(nx+1) times to solve the system. So for a complete ADI iteration, it takes O(2*(nx+1)(ny+1)) time. And for the general problem, it takes O(I2N)=O(N), where I is the number of iterations to achieve a steady-state.

Fig. 5 .

 5 Fig.5.5 shows the temperature variation results of the disk model obtained by the Douglas-

Figure 5 . 5 :

 55 Figure 5.5: Disk model with Douglas-Gunn Method

Fig. 5 .Figure 5 . 6 : 6

 5566 Fig.5.[START_REF] Lee | Periodic bicontinuous composites for high specific energy absorption[END_REF] shows the temperature variation results of the disk model obtained by ADI method. The structure contains 900 pixel cells, and there is no need to set a step size here,

Titre:

  Calcul de l'homogénéisation de la conductivité thermique des matériaux hétérogènes par la méthode des différences finiesMots clés : Matériaux hétérogènes Méthode des différences finies conductivité thermiqueRésumé : De nos jours, les matériaux hétérogènes sont de plus en plus utilisés pour leurs propriétés générales supérieures, comme les milieux poreux, qui sont largement utilisés dans les industries électroniques et biomédicales. La détermination de la conductivité thermique équivalente (CTE) des matériaux hétérogènes est donc essentielle pour la conception correcte des équipements industriels qui peuvent être soumis à des charges thermiques sévères pendant leur utilisation. L'objectif principal de cette thèse est de calculer l'homogénéisation de la conductivité thermique de matériaux hétérogènes en utilisant la méthode des différences finies (FDM). Voxel a été choisi pour la modélisation des matériaux hétérogènes et le schéma de Günter sera employé comme technique principale pour les problèmes de diffusion thermique anisotrope. Le système de Günter bidimensionnel est redémontré dans cette thèse, ainsi qu'une extension au modèle tridimensionnel, de même que des méthodes pour charger des conditions aux limites périodiques et uniformes mixtes. Les trois méthodes (FDM, FEM, et FEM+pixel(voxel)) sont comparées pour des RVEs 2D tels que des croix, des cercles et des ellipses et pour des RVEs 3D tels que des sphères et des cylindres. On découvre que la méthode FDM développée produit des résultats qui sont cohérents avec ceux de FEM et FEM+pixel(voxel) et que la méthode FDM surpasse FEM+pixel(voxel) en termes de vitesse de convergence. Cette méthode a également été appliquée à des matériaux en argent fritté pour l'étude de la conductivité thermique équivalente. Des comparaisons entre les deux méthodes (FDM et FEM) sont effectuées pour les cellules unitaires classiques telles que le cubique simple, le cubique centré sur le corps et le cubique centré sur la face, ainsi que pour le modèle stochastique à base d'argent. L'algorithme de différences finies développé est valide, et des résultats cohérents sont obtenus. En plus du schéma de Günter, un modèle à 5 points et un modèle intégral ont également été développés en s'inspirant du schéma de Günter. Pour le calcul haute performance, la bibliothèque Eigen et la bibliothèque Pardiso sont également détaillées dans la thèse. Ces deux bibliothèques contiennent à la fois des solutions directes et itératives pour la résolution d'équations linéaires. Cependant, alors qu'Eigen ne permet le calcul parallèle que de la solution itérative, Pardiso permet le calcul parallèle des deux approches, et le parallélisme est nettement supérieur à celui d'Eigen. Alors que Eigen est plus simple à construire et plus puissant, Pardiso est plus rapide pour traiter les problèmes complexes. Title: Computation homogenization of heterogeneous materials' thermal conductivity by the Finite Difference Method Keywords: Heterogeneous materials Finite Difference Method thermal conductivity

  include fixed discretization topology and size, large memory requirements, and lack of geometric representation. Due to the advantages of voxel, numerous researchers have chosen voxel models for heterogeneous materials. For example, Lebensohn et al.[START_REF] Lebensohn | N-site modeling of a 3d viscoplastic polycrystal using fast fourier transform[END_REF] first extended voxel-based FFT-based approach to crystal viscoplasticity. Kim et al.[START_REF] Kim | Voxel-based meshing and unit-cell analysis of textile composites[END_REF] pre-

sented a novel voxel-based meshing and analysis techniques based on bilinear and trilinear elements to construct textile composites. Lebensohn et al.

[START_REF] Lebensohn | Fast fourier transform-based modeling for the determination of micromechanical fields in polycrystals[END_REF]

proposed a voxel-based fast fourier transform-based modeling for the determination of micromechanical fields in polycrystals. Segurado et al.

[START_REF] Segurado | Simulation of the deformation of polycrystalline nanostructured ti by computational homogenization[END_REF]

simulated the deformation of polycrystalline nanostructured Ti with a voxel-FE model. Kabel et al.

[START_REF] Kabel | Use of composite voxels in fft-based homogenization[END_REF]

proposed a method which use of composite voxels in FFT-based homogenization. Mareau et al.

[START_REF] Mareau | Different composite voxel methods for the numerical homogenization of heterogeneous inelastic materials with fft-based techniques[END_REF]

presented different composite voxel approaches to inelastic heterogeneous materials with multiple length scales. Sossou et al.

[START_REF] Sossou | Design for 4d printing: A voxel-based modeling and simulation of smart materials[END_REF]

developed a voxel-based modeling framework for modeling and rapidly simulating the behavior of smart materials. Phung et al.

[START_REF] Phung | A voxel-based remeshing framework for the simulation of arbitrary three-dimensional crack growth in heterogeneous materials[END_REF]

developed a voxel-based meshing framework for modelling crack extension in heterogeneous materials. Berbenni et al.

[START_REF] Berbenni | A fast fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals[END_REF]

proposed a voxel-based fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals. Marano et al.

[START_REF] Marano | Fft-based simulations of slip and kink bands formation in 3d polycrystals: influence of strain gradient crystal plasticity[END_REF]

used a voxel model and an FFT-based approach to investigate the effect of strain gradients on crystal plasticity.

Table 2 .

 2

		1: Thermal conductivity of metallic materials
						20 • C
	Materials		Density	Specific heat Thermal conductivity
			(kg =m 3 ) C p (J=(kg • K))	D(W=(m • K))
	Pure aluminum		2710			902	236
	aluminum alloy (92Al-8Mg)	2610			904	107
	aluminum alloy (87Al-13Si)	2660			871	162
	Silver		10500		234	427
	Gold		19300		127	315
	Fine copper		8930			386	398
	Copper alloy (60Cu-40Ni)	8920			410	22.2
	Pure iron		7870			455	81.1
	Table 2.2: Thermal conductivity of woven composites
	Materials	Temperature	Density	In-surface thermal conductivity	Out-surface thermal conductivity

Table 2 .

 2

	3: Thermal conductivity for some Crystal at 30 • C [1]
	Crystal	system	D 11 , D 22 (W=(m • K) D 33 (W=(m • K)
	Quartz	Trigonal	6.5	11.3
	Calcite	Trigonal	4.2	5.0
	Bismuth	Trigonal	9.2	6.7
	Graphite Hexagonal	355	89

Table 3 .

 3 Imposition of BC on FD pixel/voxel models The three boundary conditions KUBC, SUBC, and MUBC, use fundamentally the Dirichlet boundary condition and the Neumann boundary condition, so in this Section, we will describe how the Dirichlet and Neumann boundary conditions are loaded.

		i-1	i	i+1
	j+1	-(a 2 *D i-1

3: coefficient for T i;j of integral isotropic scheme

  j-1;k T 16 +coef i-1;j;k T 17 + coef i-1;j+1;k T 18 + coef i;j-1;k T 20 + coef i;j;k T 21 + coef i;j+1;k T 22 +coef i+1;j-1;k T 24 + coef i+1;j;k T 25 + coef i+1;j+1;k T 26 + coef i-1;j-1;k+1 T 32 + coef i-1;j;k+1 T 33 +coef i-1;j+1;k+1 T 34 + coef i;j-1;k+1 T 36 + coef i;j;k+1 T 37 + coef i;j+1;k+1 T 38 + coef i+1;j-1;k+1 T 40

	+coef i+1;j;k+1 T 41 + coef i+1;j+1;k+1 T 42 = 0
	(3.39)

  T 24 + coef i-1;j;k+1 T 25 +coef i-1;j+1;k+1 T 26 + coef i;j-1;k+1 T 16 + coef i;j;k+1 T 17 + coef i;j+1;k+1 T 18 + coef i+1;j-1;k+1 T

	+coef i+1;j;k+1 T 21 + coef i+1;j+1;k+1 T 22 =
	(3.41)

Table 4 .

 4

			2: Equivalent thermal conductivity of layered composite
	volume
	fraction	0.04	0.12	0.2	0.28	0.36	0.44	0.52	0.6	0.68	0.76	0.84
	of fiber
	Dxx	1.04123	1.13482	1.24688	1.38351	1.55376	1.77179	2.06101	2.46305	3.05998	4.03877	5.93824
	Dyy	4.96	12.88	20.8	28.72	36.64	44.56	52.48	60.4	68.32	76.24	84.16
	Voigt	4.96	12.88	20.8	28.72	36.64	44.56	52.48	60.4	68.32	76.24	84.16
	Reuss	1.04123	1.13482	1.24688	1.38351	1.55376	1.77179	2.06101	2.46305	3.05998	4.03877	5.93824

Table 4 .

 4

		3: Equivalent thermal conductivity of isotropic cross
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2025	37.6957 6.8193e-9	5.4789e-9 37.6957	0
	FDM	2025	37.6830 3.1910e-14 1.9564e-14 37.6830 0.0003
	FEM+pixel	2025	37.6957 6.8193e-9	5.4789e-9 37.6957	0

Table 4 .

 4

[START_REF] Bika | Mechanical properties of agglomerates[END_REF]

shows the equivalent thermal conductivity of isotropic cross. According to the table, it can be found that for the cross model, when the input thermal conductivity is isotropic, the final equivalent thermal conductivity is isotropic. Secondly, the equivalent thermal conductivity obtained for the FDM is also very close to the values of the FEM and FEM+pixel for the case of 2025 cells.

Table 4 .

 4

		4: Equivalent thermal conductivity of anisotropic cross
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2025	82.3292 24.1298 24.1298 126.9628	0
	FDM	2025	82.2766 24.0740 24.0740 126.8848 0.0006
	FEM+pixel	2025	82.3292 24.1298 24.1298 126.9628	0

Table 4 .

 4 4 shows the equivalent thermal conductivity of anisotropic cross. According to the table, it can be found that for the cross model, when the input thermal conductivity is

Table 4 .

 4

		5: Equivalent thermal conductivity of isotropic disk case
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2257	1.7688 5.1364e-9	4.2535e-9 1.7688	0
	FDM	2025	1.7955 3.9454e-14 8.1789e-14 1.7955 0.0151
	FEM+pixel	2025	1.8009 5.8125e-10 3.5436e-10 1.8009 0.0182
	Table

Table 4 .

 4

[START_REF] Lee | Periodic bicontinuous composites for high specific energy absorption[END_REF]

shows the homogenised thermal conductivity obtained with the above mesh 124

Table 4 .

 4

	6: Equivalent thermal conductivity of anisotropic disk case
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	4724	8.8025 3.7809 3.7809 13.2840	0
	FDM	4225	8.8844 3.7786 3.7786 13.3964 0.0093
	FEM+pixel	4225	8.9046 3.7736 3.7736 13.4310 0.0116
	Isotropic ellipse case

Table 4 .

 4 7: Equivalent thermal conductivity of isotropic ellipse case

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2015	1.5330 4.7753e-8	4.7976e-8 1.2305	0
	FDM	2025	1.5520 2.3884e-15 1.5306e-14 1.2377 0.0124
	FEM+pixel	2025	1.5550 1.6371e-10 2.4327e-10 1.2391 0.0144
	Anisotropic ellipse case

Table 4 .

 4

	8: Equivalent thermal conductivity of anisotropic ellipse case
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	9265	7.6708 3.4101 3.4101 9.6658	0
	FDM	9025	7.6919 3.4121 3.4121 9.6975 0.0028
	FEM+pixel	9025	7.7065 3.4106 3.4106 9.7041 0.0047

Table 4 .

 4 9: Equivalent thermal conductivity of isotropic one disk and four sectors case

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	6784	3.6708 4.7092e-8 3.6482e-9 3.6708	0
	FDM	5625	3.7122 1.6851e-13 5.241e-14 3.7122 0.0113
	FEM+pixel	5625	3.7532 5.6434e-8 2.4563e-9 3.7532 0.0225

Table 4 .

 4 10: Equivalent thermal conductivity of anisotropic one disk and four sectors case

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	28702	20.2356 9.5355 9.5355 24.8120	0
	FDM	24025	20.3124 9.5112 9.5112 24.8886 0.0038
	FEM+pixel	24025	20.4105 9.5332 9.5332 24.9964 0.0086

Table 4 .

 4 11: Equivalent thermal conductivity of isotropic cross with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2025	37.6957 6.8193e-9	5.4789e-9 37.6957	0
	FDM	2025	37.6830 3.1910e-14 1.9564e-14 37.6830 0.0003
	FEM+pixel	2025	37.6957 6.8193e-9	5.4789e-9 37.6957	0

Table 4 .

 4 [START_REF] Cverna | ASM ready reference: thermal properties of metals[END_REF] shows the equivalent thermal conductivity of isotropic cross with MUBC.

Table 4 .

 4

	12: Equivalent thermal conductivity of orthotropic cross with MUBC
		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	5625	98.2706 -1.3997e-8 -4.0531e-8 151.9848	0
	FDM	5625	98.2630 -1.6352e-8 -1.1531e-8 151.9726 7.7337e-5
	FEM+pixel	5625	98.2706 -1.3997e-8 -4.0531e-8 151.9848	0
				146

Table 4 .

 4 13: Coefficient of linear trend line for the isotropic disk model

	y=ax+b	FDM 0	FDM 1	FDM 2	FDM 3	FDM 4	FEM+pixel	Geometry
	a	-0.60725365	-0.60337053	-0.60490769	-0.60534945	-0.60669458	-0.57979198	-0.88877376
	b	0.10134371	0.08756823	0.09272454	0.09427982	0.09977462	0.13237423	0.66552513

Table 4 .

 4 14: Equivalent thermal conductivity of an isotropic disk case with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	1764	1.7688 9.0216e-9	3.6979e-9 1.7688	0
	FDM	1600	1.7759 3.1629e-12 1.3452e-12 1.7759 0.00403
	FEM+pixel	1600	1.7853 3.7623e-10 4.1777e-10 1.7853 0.00935
	for an isotropic disk inclusion with a radius of 0.3 are depicted in Fig.4.35, where FDM
	and FEM+pixel each have 1600 cells and FEM has 1764 cells. Comparing the overall

Table 4 .

 4 15: Coefficient of linear trend line for the orthotropic disk model

	y=ax+b	FDM 0	FDM 1	FDM 2	FDM 3	FDM 4	FEM+pixel	Geometry
	a	-0.60191387	-0.58978822	-0.59453232	-0.59771438	-0.59906097	-0.57286603	-0.88877376
	b	0.07773614	0.03096675	0.04862944	0.06104253	0.06665823	0.09260261	0.66552513

Table 4 .

 4 [START_REF] Hashin | A variational approach to the theory of the effective magnetic permeability of multiphase materials[END_REF]: Equivalent thermal conductivity of an orthotropic disk case with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	4042	8.9974 6.0216e-9	5.2314e-9 13.6238	0
	FDM	3600	9.0668 1.2366e-10 2.5367e-10 13.7222 0.00771
	FEM+pixel	3600	9.0894 5.3792e-10 2.1014e-9 13.7573 0.01022

Table 4 .

 4 17: Coefficient of linear trend line for the isotropic ellipse model

	y=ax+b	FDM 0	FDM 1	FDM 2	FDM 3	FDM 4	FEM+pixel	Geometry
	a	-0.55525695	-0.56054928	-0.55770636	-0.55518292	-0.55702431	-0.47442006	-0.66190383
	b	-0.15879611	-0.13766961	-0.14912848	-0.15901608	-0.1517911	-0.25934387	0.27896002

Table 4 .

 4 18: Equivalent thermal conductivity of an isotropic ellipse case with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	1257	1.5330 2.8993e-7	8.4788e-8 1.2305	0
	FDM	1225	1.5436 3.1092e-13 1.2073e-13 1.2444 0.0069
	FEM+pixel	1225	1.5518 3.2252e-9 9.8487e-10 1.2467 0.0123
	Orthotropic ellipse case

Table 4 .

 4 19: Coefficient of linear trend line for the orthotropic ellipse model

	y=ax+b	FDM 0	FDM 1	FDM 2	FDM 3	FDM 4	FEM+pixel	Geometry
	a	-0.55214852	-0.55865893	-0.55326616	-0.54663055	-0.55308127	-0.46282871	-0.66190383
	b	-0.19898609	-0.17299922	-0.194714	-0.22070744	-0.19532906	-0.3348816	0.27896002

Table 4 .

 4 20: Equivalent thermal conductivity of an orthotropic ellipse case with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	2428	7.8935 1.5548e-8	7.8910e-8	9.7269	0
	FDM	2500	7.9321 1.7834e-11 9.3702e-12 9.6855 0.0049
	FEM+pixel	2500	7.9622 5.95586e-9 3.38902e-10 9.6972 0.0087
				161

Table 4 .

 4 22: Equivalent thermal conductivity of an isotropic One disk and four sectors case with MUBC discs of radius equal to 0.3 for the inclusions is shown in Fig.4.47. The FDM and FEM+pixel have 1600 cells each, while the FEM has 1883 cells. By comparison,

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	1883	3.6708 4.6623e-8	8.3775e-8 3.6708	0
	FDM	1600	3.6768 1.0934e-09 2.2057e-10 3.6768 0.0016
	FEM+pixel	1600	3.7551 2.4192e-10 3.3040e-10 3.7551 0.0230
	and four quarter

Table 4 .

 4 23: Coefficient of linear trend line for the orthotropic one disk and four sectors model

	y=ax+b	FDM 0	FDM 1	FDM 2	FDM 3	FDM 4	FEM+pixel	Geometry
	a	-0.45436224	-0.45559613	-0.45868457	-0.45747664	-0.45765636	-0.50867393	-0.7652408
	b	-0.13628578	-0.12835056	-0.1140718	-0.12042247	-0.11912843	0.37046881	0.18567703

Table 4 .

 4 24: Equivalent thermal conductivity of an orthotropic one disk and four sectors case with MUBC

		Number of cells	D xx	D xy	D y x	D y y	Error of D xx
	FEM	7484	20.5418 1.5279e-8	1.2570e-8 25.1330	0
	FDM	6400	20.6068 1.4852e-08 7.5961e-9 25.1867 0.0032
	FEM+pixel	6400	20.8319 3.2252e-9 7.3518e-10 25.4296 0.0141

Table 4 .

 4 [START_REF] Hill | A self-consistent mechanics of composite materials[END_REF] shows the equivalent thermal conductivity of an orthotropic one disk and four sectors case with MUBC. After comparison, we can deduce that when the input thermal conductivity of a disk and four sector model is orthotropic, the final equivalent thermal conductivity is orthotropic. When the errors are compared to the FEM values, it is clear that the FDM values are far smaller than the FEM+pixel values.

Table 4 .

 4 25: Equivalent thermal conductivity of 3D layered composite

	volume
	fraction	0.04	0.12	0.2	0.28	0.36	0.44	0.52	0.6	0.68	0.76	0.84
	of fiber
	Dxx	1.04123	1.13482	1.24688	1.38351	1.55376	1.77179	2.06101	2.46305	3.05998	4.03877	5.93824
	Dyy	4.96	12.88	20.8	28.72	36.64	44.56	52.48	60.4	68.32	76.24	84.16
	Dzz	4.96	12.88	20.8	28.72	36.64	44.56	52.48	60.4	68.32	76.24	84.16
	Voigt	4.96	12.88	20.8	28.72	36.64	44.56	52.48	60.4	68.32	76.24	84.16
	Reuss	1.04123	1.13482	1.24688	1.38351	1.55376	1.77179	2.06101	2.46305	3.05998	4.03877	5.93824

Table 4 .

 4

	26: Comparison of trendline parameters for Ball inclusion model with different
	R/L
	R/L	0.1	0.15	0.2	0.25	0.3	0.35	0.4
	Slope of FDM	-0.3287 -0.36850 -0.3469 -0.3241 -0.35286 -0.3034 -0.3567
	Slope of FEM +voxel	-0.2866 -0.2137	-0.3348 -0.3214 -0.3712	-0.3265 -0.3670
	Intercept of FDM -1.1906 -0.7613	-0.6071 -0.5393 -0.1906	-0.3751 0.0294
	Intercept of FEM +voxel	-1.1227 -1.2342	-0.3505 -0.2599 0.1484	0.0033 0.3600

Table 4 .

 4 27: Comparison of trendline parameters for Cylindre inclusion model with different

	R/L
	R/L	0.1	0.15	0.2	0.25	0.3	0.35	0.4
	Slope of FDM	-0.3525 -0.6052 -0.2941 -0.5224 -0.4359 -0.5807 -0.4396
	Slope of FEM +voxel	-0.3660 -0.5477 -0.3308 -0.3511 -0.4078 -0.4744 -0.4111
	Intercept of FDM -0.6632 0.5845 -0.7150 0.2399 0.20779 0.8697	0.5125
	Intercept of FEM +voxel	-0.3977 0.3994 -0.3003 -0.2908 0.2768	0.52634 0.5333

Table 4 .

 4 28: Convergence results of different mesh levels by FED and FEM.

	Number of voxel cells by FDM	15625	91125	274625 614125 1953125
	Equivalent thermal conductivity 233.092 232.035 229.258 226.226 224.789
	Number of finite elements	17677	36961	132968 446553 2348718
	Equivalent thermal conductivity 227.54	226.1	224.38	223.86	223.35

  By comparison, we select SparseLU as our direct solver, ConjugateGradient as our iterative solver, and BiCGSTAB if our matrix is asymmetric. By default, iterations start with x = 0 as the initial estimate of the solution. We can control the startup using the solveWith-Guess method. In addition, in the iterative solver, if we compile user code with OpenMP enabled, we can take advantage of multithreading.

	Time (s)	10	CG thread1 CG thread2 CG thread4 CG thread6 CG thread8 CG thread12		Time (s)	10 1 10 0 10 1 10 2 10 3	CG_is thread1 CG_is thread2 CG_is thread4 CG_is thread6 CG_is thread8 CG_is thread12
		10					10 2
		10					10 3
		10 2	10 3	10 4 Number of DOF 10 5	10 6	10 7	10 2	10 3	10 4 Number of DOF 10 5	10 6	10 7
			(a) CG with different threads		(b) CG initial guess with different threads
	4.9 Comparison of solvers 4.9.1 The different solvers of Eigen 10 1 10 10 10 10 Time (s) CGD thread1 CGD thread2 CGD thread4 CGD thread6 CGD thread8 CGD thread12 10 1 10 0 10 1 10 2 10 3 Time (s)	CGD_is thread1 CGD_is thread2 CGD_is thread4 CGD_is thread6 CGD_is thread8 CGD_is thread12
	10 2					10 2
	10 3					10 3
		10 2	10 3	10 4 Number of DOF 10 5	10 6	10 7	10 2	10 3	10 4 Number of DOF 10 5	10 6	10 7
			(c) CGD with different threads		(d) CGD initial guess with different
	Then we compare SparseLU, ConjugateGradient, ConjugateGradient with initial assump-
	tion, ConjugateGradient with DiagonalPreconditioner, and ConjugateGradient with Diag-
	onalPreconditioner with initial assumption with different threads.
	Fig.4.77 shows the solution speed for different iterative solutions with different threads.
						201

  Douglas-Gunn methodDouglas and Gunn developed another ADI method based on scheme Crank-Nicolson; this method is solved by two sub-steps:

	(1 -	fl x 2	‹ 2 x )T i;j n+1=2	= (1 +	fl x 2	‹ 2 x + fl y ‹ 2 y )T n i;j
	(1 -	fl y 2	‹ 2 y )T n+1 i;j	= (1 +	fl x 2	‹ 2 x +	fl y 2	‹ 2 y )T n i;j +	fl x 2	‹ 2 x T i;j n+1=2

;j+[START_REF]Heat Conduction in Anisotropic Solids[END_REF]

;j+ 1 2

;j-1

;j+ 1 2 ;k+ (16∆y ∆z)

Acknowledgements

This work was funded by the China Scholarship Council (CSC). Firstly, we would like to thank my supervisor, Zhiqiang Feng, who gave me this opportunity to pursue my PhD three years ago and with his help succeeded in getting the CSC funding and helped me plan my direction as well as offered valuable advice and direction during my numerous

*T i;j+1 *T i;j+1 *T i;j+1 )/(dx*dy*T i+1;j+1 *T i;

Chapter 4

Comparison of numerical results and the application of REV to complex heterogeneous materials

Introduction

In chapter 2 we introduce the two-dimensional Günter model and the extension to the three-dimensional model. In chapter 3 we present methods for loading periodic boundary conditions and for mixing homogeneous boundary conditions. Our main objective is to deal with homogenization, which means solving for the equivalent properties of the composite, and for the moment we only consider the thermal case.

In this chapter, we present our numerical results, which include different classical heterogeneous structural models in 1D-2D-3D. We also present the results of the finite element method and the analytical solution for comparison. We verify the accuracy of our model calculations by comparing the calculated temperature fields, temperature gradients and heat flows under different models. The equivalent thermal conductivity of these classical and (c), it is clear that the overall trend of variation is identical for all three techniques.

(d)-(i) illustrate the temperature gradient and heat flow fields, respectively, and it can be seen that the overall trend is similar for all three methods, with the errors concentrated primarily in the region between the inclusion and matrix phases, which is primarily caused by the thermal conductivity discontinuity.

The comparison of path data of the isotropic case of one disk inclusion of radius r equal to 0.3 for thermal periodic boundary value problem is shown in Fig. 4.14. We have mainly chosen the two paths with y = 0:5 and y = x for comparison. When the temperatures on the paths in (a) and (b) are compared, it is clear that the three approaches agree well.

Comparing temperature gradients along the paths depicted in (c) and (d) reveals that the differences are mostly concentrated at the intersections of the inclusion and matrix phases.

In the comparison of heat flux depicted in (e) and (f), (d) shows similar results for FDM and FEM+pixel on the points between the inclusion and matrix phases, but with differences from the FEM results. This disparity also causes a discrepancy between the FEM results and the other two methods for points within the inclusion phase. However, the FDM is slightly better than the FEM+pixel. Whereas (e) shows more FEM+pixel problems at corner points, showing significant differences, FDM does not show such problems. Among them, 5625 pixel cells are included in FDM and FEM+pixel, and 6784 quadrilateral cells are included in FEM.

Isotropic one disk and four sectors case

Isotropic one disk and four sectors case

The convergence analysis of One disk and four sectors inclusions of radius r= 0.

Evolution of the equivalent thermal conductivity

In this section, we calculate the equivalent thermal conductivity of the aforementioned models by using Eq.2.47 for various volume fraction of silver. To make the volume fraction of silver vary in a given unit cell (SC, BCC or FCC), the positions of the corresponding silver particles are kept unchanged while their radii are uniformly increased to reduce the voids. The results of the simulations are shown in Fig. 4.71. We draw the following observations:

(1) The equivalent thermal conductivities of SC, BCC, FCC, and their conjugated models increase with the increase of the silver fraction. This is consistent with the fact that the conductivity of silver (429 W• m -1 •K -1 ) is larger than the air's one (0.0257 W• m -1 •K -1 ).

(2) From Fig. 4.71a,c,e, one can find that as the silver volume fraction increases, a better agreement between the FDM and the FEM results is found. The largest difference between the two methods in the three models occurs for the lowest values of the silver volume fraction when the silver particles are just contacting each other in the respective unit cells (i.e. the radius of spherical silvers for SC model is 0.5, for BCC model is √ 3=4 ≈ 0:433, for FCC model is √ 2=4 ≈ 0:354; and the corresponding volume fraction of silver for SC model is 0.555, for BCC model is 0.713, for FCC model is 0.779). This is caused by the inherent property of the voxel meshing which has poor capabilities for approximating curved surfaces like those of the spherical silver particles. Voxel meshing produces sawtooth shaped surfaces when mesh density is small. By increasing the mesh density, the values of FDM gradually tend to those of FEM. It is to say that the FDM is more sensitive to lower volume fractions of silver, especially when the spheres are merely contacting.

(3) In contrast, from Fig. 4.71b,d,f, one could not find the similar sensitivities. It is because for the conjugated models, the silvers regions are connected in a different way.

The conjugated connections ensure that the silver phase forms always a contiguous part.

Stochastic granular models

To further test the usability of the developed FD algorithm in other situations, we additionally consider a kind of stochastic models. Fig. 4.74 illustrates the geometrical models of the stochastic structure. The structures are set to unit size and composed of two phases: silver and air. In order to make the spherical silver particles connect with each other to form a continuous entirety, an overlap is proposed herein (for example, this can mimic the neck formation between silver nanoparticles after the sintering process). The generation method of these spheres is discussed in [START_REF] Wang | Numerical modeling of low-temperature and low-pressure sintering of silver microparticles based on surface and grain boundary diffusion mechanisms[END_REF][START_REF] Benabou | Simulation of silver nanoparticles sintering at high temperatures based on theoretical evaluations of surface and grain boundary mobilities[END_REF]. Here, 56 spheres are generated, and their radii obey a normal distribution. In Figs.4.74, three cases are considered: Case 1 is a stochastic model with minimal overlap between particles, having an average particle radius of 0.139 and a standard deviation of 0.026 (silver fraction = 0.663); Case 2 is a stochastic model with a larger overlap, having an average particle radius of 0.146 and a standard deviation of 0.028 (silver fraction = 0.734); and Case 3 is the conjugated structure of Case 2 with an average particle radius of 0.059 and a standard deviation of 0.011, i.e., keeping the silver fraction constant (silver fraction = 0.734) and replacing the silver phase with air and the air phase with silver (this can mimic the situation of air bubbles in solids).

To ensure the periodicity of these models, spheres that go out from one surface will be enforced to re-enter from the opposite surface. The simulation results of the equivalent thermal conductivity are shown in Table 4.29. By comparing Case 1 and Case 2, we find that the equivalent thermal conductivity obtained by the two methods are very similar, and the finite difference method is closer to the finite element method as the overlap portion increases. This result also verifies our when DOF is less than 30 6 , the rate of the iterative solution and the direct solution is almost the same. However, as the DOF gets larger, the speed of the direct solution is slightly faster than that of the iterative solution.

Comparsion Eigen and Pardiso

First of all, both solvers are free for academics. Secondly, Pardiso uses the CSR format to input a sparse matrix. CSR matrix compression: Three arrays are used to represent a two-dimensional matrix. The three one-dimensional arrays of CSR are represented by A, IA, and JA, respectively. NNZ represents the number of non-zero elements of the sparse matrix M (mxn). A (non-zero values), IA (the extents of rows), JA (column indices). Eigen can use a variety of methods to input sparse matrices. One is to create a list of triplets with an element type and then convert it to a sparse matrix. The so-called triplets are a non-zero-valued data structure defined in Eigen for storing sparse matrices. The storage format is (i, j, value). Three parameters are commonly used to create functions, namely rows, columns, and values. In addition to using Triplet, a more "simple and rude" method is to call the member function .insert() of the sparse matrix to directly insert the value.

The advantage of this is that there is no need to create a new Triplet object and list, which may have higher performance and less memory usage. It should be noted that the method of directly inserting numerical values requires pre-setting the number of non-zero entries in the matrix. Otherwise, every time a numerical value is inserted, Eigen will reapply for space to save, which consumes a lot of time.

Comparsion speed of Eigen and Pardiso

With Figure 4.81, we see that the computation speed of Pardiso is much higher than that of Eigen. Eigen is faster than Pardiso only when the DOF is less than 1e4, but since the time is too short (much less than 1 second), the difference is not large. However, as the conditions and various loading methods for Neumann boundary conditions for the purpose of resolving the final matrix singularity problem. Additionally, methods for loading periodic boundary conditions are described, as well as a variation method that permits the matrix to remain symmetrical after loading periodic boundary conditions.

Chapter 4 compares FDM, FEM, and FEM+pixel (voxel). The results demonstrate that the FDM approach beats the FEM+pixel(voxel) method in both convergence speed and accuracy. Additionally, the approach has been applied to the REV of complicated heterogeneous materials, such as sintered silver, with the following results. First, our proposed algorithm yielded comparable results with the FEM, whether dealing with the classical SC, BCC, FCC models, or complex stochastic models. Second, the calculation of the equivalent thermal conductivity of the conjugated structures is more accurate than the typical structures for all geometrical models that we used. Finally, the FDM is more effective in dealing with conjugated structures because it is less sensitive to the density of mesh. This point indicates that the developed FDM is suitable for handling this type of material. A comparison between Eigen and Pardiso is given at the end and it can be seen that Eigen is usually parallel. In contrast, Pardiso can be well parallelized and shows a much higher speed than Eigen when dealing with large solving problems.

In Chapter 5, we present complementary research, mainly consisting of dynamic models attempted to solve the matrix singularity problem, including the Euler explicit scheme, the Euler implicit scheme, the Crank-Nicolson method et ADI method. All of these schemes solve the matrix singularity problem. For example, the explicit scheme does not require a matrix solution but has a strict step size requirement. The rest of the methods are all implicit schemes, with ADI speed being the only priority. Finally, we introduce an alternating direction iteration method, which uses alternating direction iteration to solve steady-state problems, effectively increasing the solution speed and reducing the amount of computer memory consumed by the solution.

In the following, we give our perspectives. In Chapter 3,1 we give the new models developed after receiving inspiration from the Günter model, and so far the new models developed are all second-order models. Work on the extension of the new models can continue, for example, with the development of fourth-order models.

Secondly, in this area of ADI parallel computing, the construction of a 3D model has not yet been completed for this thesis, and this part will also be continued as a postdoctoral project. Using ADI parallelism, the computation of some very large problems has been realized.

Appendix A

3D isotropic scheme

+ (T i+1;j+1;k+1 + T i+1;j;k+1 + T i;j+1;k+1 + T i;j;k+1 -T i+1;j+1;k -T i+1;j;k -T i;j+1;k -T i;j;k ) ∆x ∆y 2 ) D i+1=2;j+1=2;k+1=2 +(((T i+1;j;k+1 + T i+1;j-1;k+1 -T i;j;k+1 -T i;j-1;k+1 + T i+1;j;k + T i+1;j-1;k -T i;j;k -T i;j-1;k ) ∆y 2 + (-T i+1;j;k+1 + T i+1;j-1;k+1 -T i;j;k+1 + T i;j-1;k+1 -T i+1;j;k + T i+1;j-1;k -T i;j;k + T i;j-1;k ) ∆x )∆z 2 + (T i+1;j;k+1 + T i+1;j-1;k+1 + T i;j;k+1

+ (T i;j+1;k+1 + T i;j;k+1 + T i-1;j+1;k+1 + T i-1;j;k+1 -T i;j+1;k -T i;j;k -T i-1;j+1;k -T i-1;j;k ) ∆x ∆y 2 ) D i-1=2;j+1=2;k+1=2 +(((-T i;j;k+1 -T i;j-1;k+1 + T i-1;j;k+1 + T i-1;j-1;k+1 -T i;j;k -T i;j-1;k + T i-1;j;k + T i-1;j-1;k ) ∆y 2 + (-T i;j;k+1 + T i;j-1;k+1 -T i-1;j;k+1 + T i-1;j-1;k+1 -T i;j;k + T i;j-1;k -T i-1;j;k + T i-1;j-1;k ) ∆x )∆z 2 + (T i;j;k+1 + T i;j-1;k+1 + T i-1;j;k+1 + T i-1;j-1;k+1 -T i;j;k -T i;j-1;k -T i-1;j;k -T i-1;j-1;k ) ∆x ∆y

+(((-T i;j;k-1 -T i;j-1;k-1 + T i-1;j;k-1 -T i;j;k -T i;j-1;k + T i-1;j;k + T i-1;j-1;k + T i-1;j-1;k-1 ) ∆y 2

Appendix B 3D anisotropic scheme ∇ • -→ q = (-k13 i+ 1 2 ;j+ 1 2 ;k+ (-k13 i-1 2 ;j+ 1 2 ;k-1 2 (16∆x∆z) -k13 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆z) -k12 i-1 2 ;j+ 1 2 ;k-1 2 (16∆x∆y ) + k12 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k11 i-1 2 ;j+ 1 2 ;k-(16∆x∆x) -k11 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆x) + k23 i-1 2 ;j+ 1 2 ;k-1 2 (16∆y ∆z) -k23 i-1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) + k22 i-1 2 ;j+ 1 2 ;k-1 2 (16∆y ∆y ) + k22 i-1 2 ;j-1 2 ;k-(16∆y ∆y ) + k21 i-1 2 ;j+ 1 2 ;k-1 2 (16∆x∆y ) -k21 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k33 i-1 2 ;j+ 1 2 ;k-1 2 (16∆z∆z) -k33 i-1 2 ;j-1 2 ;k-1 2 (16∆z∆z) -k32 i-1 2 ;j+ 1 2 ;k-(16∆y ∆z) + k32 i-1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k31 i-1 2 ;j+ 1 2 ;k-1 2 (16∆x∆z) -k31 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆z) )T i-1;j;k-1 +(-k13 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆z) -k12 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k11 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆x) -k23 i-1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k22 i-1 2 ;j-1 2 ;k-(16∆y ∆y ) -k21 i-1 2 ;j-1 2 ;k-1 2 (16∆x∆y ) -k33 i-1 2 ;j-1 2 ;k-1 2 (16∆z∆z) -k32 i-1 2 ;j-1 2 ;k-1 2 (16∆y ∆z) -k31 i-1 2 ;j-1 2 ;k-