# An efficient strategy for large scale 3D simulation of heterogeneous materials to predict effective thermal conductivity

Xiaodong Liu, Julien Réthoré, Marie-Christine Baietto, Philippe Sainsot, Antonius Adrianus Lubrecht

<figure><img src="image_1.png" /></figure>

## To cite this version:

Xiaodong Liu, Julien Réthoré, Marie-Christine Baietto, Philippe Sainsot, Antonius Adrianus Lubrecht. An efficient strategy for large scale 3D simulation of heterogeneous materials to predict effective thermal conductivity. Computational Materials Science, 2019, 166, pp.265-275. 10.1016/j.commatsci.2019.05.004. hal-02150163

HAL Id: hal-02150163

https://hal.science/hal-02150163

Submitted on 21 Oct 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An efficient strategy for large scale 3D simulation of heterogeneous materials to predict effective thermal conductivity

Xiaodong Liu<sup>a</sup>, Julien Réthoré<sup>a,*</sup>, Marie-Christine Baietto<sup>b</sup>, Philippe Sainsot<sup>b</sup>, Antonius Adrianus Lubrecht<sup>b</sup>

<sup>a</sup> Research Institute in Civil Engineering and Mechanics (GeM), CNRS UMR 6183, École Centrale de Nantes, 1 rue de la Noë, 44321 Nantes Cedex 3, France

<sup>b</sup> Univ Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France

The behavior of heterogeneous materials (e.g. composite, polycrystalline) is one of the most challenging subjects for researchers. Imaging techniques based on X-ray tomography permit one to obtain the inner structure of these materials. To account for such detailed information and understand the material behavior, the use of tomographic images as an input for numerical simulations is becoming more and more common. The main difficulty is the computational cost, mesh generation in the context of finite element simulations and the large discontinuities of the material properties, which can lead to divergence. The subject of this paper is to resolve these problems and understand the behavior of a material with high heterogeneity. The application of a MultiGrid method coupled with homogenization technique is proposed. The MultiGrid method is a well known method to increase convergence speed. The homogenization technique is applied to compute the coarse grid operators of the MultiGrid process. The combination of the two methods can efficiently deal with large property variations. Hybrid MPI/OpenMP parallel computing is used to save computational time. The influence of material heterogeneity is analyzed as well as the ratio of material properties for the thermal conduction. The effective thermal conductivity of material is illustrated.

## 1. Introduction

The use of composite materials in many industrial fields has become more and more widespread during the last decade. It is well known that many composite materials exhibit an excellent mechanical behavior. However, due to the complex structures and variable components of composite materials, it is not simple to understand their properties, which limits the application of these materials. Traditional experiments permit to globally characterize some of their properties. Nevertheless, it can just reduce this problem but not solve it fully. The characterization of the thermal properties of a composite material is an example.

Homogenization techniques have been developed to obtain effective thermal and mechanical properties of heterogeneous materials [8,26,14]. The limitation of this method is that except for periodic materials, for which (under periodic boundary conditions) the prediction is robust, only upper and lower bounds of the effective properties can be obtained. For instance, for materials with high heterogeneity or with an irregular structure, this method has a poor performance. More importantly, the manufacturing process can neither ensure a defect-free material nor an exact geometry of the phases e.g. cast iron. In addition,

many materials, e.g. polycrystalline materials, have an intricate structure. To sum up, one needs the real inner structure of the material to understand its behavior.

Fortunately, imaging techniques based on X-ray tomography show the inner structure of materials [13]. With these details, one can better understand the material behavior. Motivated by the secret of the material behavior, e.g. mechanical and thermal properties, using real tomographic images as an input to perform numerical simulations is under development. Much work has been devoted to this subject. The work of Lengsfeld et al. [19] and Bessho et al. [2] presented the numerical simulation of bone tomography. They studied mechanical problems e.g. hip fractures of the human femur, using Finite Element Methods (FEM). Ferrant et al. [11], Michailidis et al. [20] and Proudhon et al. [23] also applied FEM simulations to tomographic images of industrial materials to analyze their properties. As we saw previously, FEM is widely used for this kind of simulations. Nevertheless, the mesh generation of FEM needs human intervention, which is time consuming. The work of Gu et al. [15] that introduced a 3D simulation of the elastic behavior of a laminated composite material, proposed to use a Finite Difference Method (FDM), to take one voxel per

<footnote>* Corresponding author.
E-mail address: julien.rethore@ec-nantes.fr (J. Réthoré).</footnote>
grid to avoid heavy human work in the meshing step, and to use the MultiGrid (MG) method to accelerate the convergence speed. Besides numerical simulation, for instance, another choice to understand the material behavior, is to use tomographic techniques during the experiment. e.g. the work of [25] studied the heat conductivity for materials with complex structures. However, to experimentally measure the heat conductivity of materials is still difficult, especially, for materials with a complicated structure. e.g. composite materials. It is still a challenge due to several kinds of problems e.g. thermal diffusion, thermal expansion during the experiment, material anisotropy.

Motivated by these practical considerations, the development of a standard process to obtain the effective thermal conductivity of heterogeneous materials received considerable attention. This work takes the tomographic image as an input to a thermal conduction simulation to study the material thermal behavior and to obtain its effective conductivity.

The remainder of this paper is organized as follows. Section 2 reviews the background, equations and some fundamentals of FEM. Section 3 presents the strategy for an efficient large scale simulation; Matrix Free FEM, MG method, homogenization techniques and parallel computing are illustrated in this section. Applications and results of this work are described in Section 4. A discussion and conclusion section is provided in Section 5.

## 2. Problem statement

### 2.1. Governing equations

Thermal conduction can be treated by a heat equation according to the first law of thermodynamics (i.e. conservation of energy):

$$ \rho c_p \frac{\partial T}{\partial t} - \nabla \cdot (\alpha \nabla T) = \dot{q}_v \quad (1) $$

where:

* $\rho$ is mass density of material

* $c_p$ is specific heat capacity

* $T$ is temperature

* $t$ is time

* $\nabla$ denotes the gradient operator

* $\alpha$ is the thermal conduction coefficient which is a second order tensor

* $\dot{q}_v$ is the volumetric heat source.

Since the focus of this work is thermal conductivity, it is considered that there is no extra source and the thermal field does not depend on time. The heat Eq. 1 becomes a typical Poisson equation:

$$ \nabla \cdot (\alpha \nabla T) = \operatorname{div}(\alpha \nabla T) = 0 \quad (2) $$

The tomographic images used in this study are cubic $\Omega \in \mathbb{R}^3$. Two kinds of boundary conditions are applied on $\partial\Omega$ as illustrated in Fig. 1:

* Dirichlet boundary condition (i.e. prescribed temperature) on $\Gamma_1$ and $\Gamma_2$

* Neumann boundary condition (i.e. prescribed heat flux) on the other four surfaces.

which can be written as Eq. 3:

$$ \left\{ \begin{array}{ll} T = T_0 & \text{on } \Gamma_1 \\ T = T_1 & \text{on } \Gamma_2 \\ \alpha \nabla T \cdot \mathbf{n} = 0 & \text{on the other surfaces} \end{array} \right. \quad (3) $$

### 2.2. Finite element discretization

FEM is one of the most common methods to discretize $\Omega$ and solve

<figure><img src="image_3.png" /><figcaption>Fig. 1. Boundary conditions applied on $\partial\Omega$.</figcaption></figure>

the governing equations. However, the images representing the inner structure of the material have a very complex shape. The use of standard meshes conforming to the phase geometry, requires much human work, as mentioned before and stated in the work of many researchers, e.g. Lengsfeld et al. [19], Bessho et al. [2], Ferrant et al. [11], Michailidis et al. [20], Proudhon et al. [23] and Nguyen et al. [21]. The strategy to use one node per voxel in the image has been chosen to avoid this difficulty. That means to assign the material property in each voxel on each elementary node. As the domain $\Omega$ is cubic and each voxel is cubic, 8-node cubic elements are chosen to discretize $\Omega$.

Multiplying Eq. (2) with a test function and integrating over $\Omega$, one obtains:

$$ \int_{\Omega} \operatorname{div}(-\alpha \nabla T) \varphi \, d\Omega = 0 \quad (4) $$

where $\varphi$ is the test function.

Applying integration by parts, the formula reads:

$$ - \int_{\partial\Omega} \alpha \nabla T \cdot \vec{n} \varphi \, dS + \int_{\Omega} \alpha \nabla T \cdot \nabla \varphi \, d\Omega = 0 \quad (5) $$

where $-\alpha \nabla T \cdot \vec{n}$ is the heat flux in the outward normal direction $\vec{n}$ on the boundary. Eq. 5 can be summarized as:

$$ \int_{\partial\Omega} \alpha \nabla T \cdot \vec{n} \varphi \, dS = \int_{\Omega} \alpha \nabla T \cdot \nabla \varphi \, d\Omega \quad (6) $$

Employing finite element discretization, one gets:

$$ \vec{T} \approx \vec{T}_h = \sum_{i=1}^{N} T_i \varphi_i \quad (7) $$

where $\vec{T}_h = \{T_1, T_2, ..., T_N\}$ is an approximate solution of $T$, $N$ denotes the number of unknowns and $\varphi_i$ is the shape functions of 8-node cubic elements, which is the same as the test function.

Finally, Eq. (6) can be written as:

$$ \sum_{i=1}^{N} T_i \int_{\Omega} \alpha \nabla \varphi_i \cdot \nabla \varphi_j \, dS = \sum_{i=1}^{N} \int_{\partial \Omega} (\alpha T_i \nabla \varphi_i) \cdot \vec{n} \, \varphi_j \, dS \quad (8) $$

After using Gauss integration on the integral above, $\alpha$ used in the FEM equation becomes the conductivity in each Gauss integration point $\alpha_g$. It is obtained by the interpolation with shape functions, and, described as:

$$ \alpha_g = \sum_{s=1}^{8} \alpha_s \cdot \varphi_{s,g} \quad (9) $$

where $s$ denotes the shape function assigned with each elementary node, $\varphi_{s,g}$ the value of shape functions at Gauss integration points and $\alpha_s$ elementary node conductivity.

From Eq. (8), one obtains a system of equations:

$$ \mathbf{L} \vec{T}_h = \vec{f} \quad (10) $$

where $\mathbf{L}$ is a matrix which is often referred to as the stiffness matrix, $\vec{T}_h$
is a vector containing all unknowns (temperature at each node) and $\vec{f}$ is the right hand side vector.

The aim is to solve Eq. (10) to obtain the unknowns $T_i$.

# 3. Implementation

## 3.1. Iterative solution

The best known solution algorithms for Eq. (10) are direct solvers and iterative solvers. For the direct solver, one has to assemble the stiffness matrix, which requires a large memory. The largest image that will be computed in this application is an image containing 2049 × 2049 × 2049 voxels, which means that the number of elements is 2048³ (i.e. more than eight billion elements). Supposing one uses cubic elements with 8 Gauss integration points, the size of global sparse matrix is 2049 × 2049 × 2049 × 27 × 8bytes ≈ 1.69TB. It is impossible to have such a huge memory space available on a normal computer. The size of the stiffness matrix does not allow one to assemble the whole matrix. It forces one to use an iterative solver without assembling the stiffness matrix, which is often called the Matrix Free Finite Element Method (MF-FEM). The work of Carey and Jiang [7] invented MF-FEM. For instance, this method is developed and widely used, especially, for parallel computing (see e.g. Tezduyar et al. [28]). Flair and Arbenz [12] also developed a matrix-free MultiGrid solver for tomographic image simulation. In this work, one proposes to use a MF-FEM with the diagonal value of the stiffness matrix to accelerate the convergence rate.

## 3.2. MultiGrid and homogenization

The iterative single level Jacobi solver applied previously can quickly decrease the high frequency components of the error, but for low frequency errors, it does a poor job. The convergence speed diminishes rapidly as presented in Fig. 5a.

It is well known that the MG method is one of the most efficient ways to increase the convergence rate. The idea of MG is to construct several levels or grids. Then, iterative relaxations are carried out at each level, high-frequency errors can be eliminated on fine grids and low-frequency errors can be eliminated on coarse grids (e.g. work of [5,6]). With this method, one can solve the slowness problem of single level Jacobi relaxation due to the presence of low-frequency errors. The work of Biboulet et al. [3] shows the efficiency of using a MG method on the FEM. However, they assembled the stiffness matrix which, is very expensive for large scale problems.

Nevertheless, a standard MG method is not adapted for problems with high heterogeneity. It has a very poor convergence performance, when large variations of the material properties are to be considered, or rather, high temperature gradients on coarse grids are involved. These variations make the linear prolongation and restriction operators almost ineffective. However, the material property on a coarse grid is unknown and it is should be chosen to avoid the poor performance of classical coarse grid and inter-grid operators.

To resolve these problems, one needs new coarse grid operators, prolongation operators and restriction operators. Several researchers have investigated this problem, such as the work of Alcouffe et al. [1] for 2D, Hoekema et al. [17] for 3D, Engquist and Luo [9] and Engquist and Luo [10]. These researchers proposed several methods to alleviate the poor convergence of the standard MG method. But the problem is that the implementation of these ideas is not simple. The computational time and memory cost are the two other limitations. Based on the work of Alcouffe et al. [1] and Sviercoski et al. [27], some new operators for the MG method are proposed in this paper.

Sviercoski et al. [27] proposed to use a Cardwell and Parsons (CP) bounds type homogenization to get the analytical coarse grid operator. The idea is to compute the upper and lower CP bounds of the material property on each coarse grid from the finest grid. After that, the average

of the arithmetic and geometric averages of the CP bounds, is supposed to be the effective property on each coarse grid. With this strategy, one can obtain the diagonal components of the material property tensor, which is already sufficient for isotropic materials. For anisotropic materials, they proposed a method to calculate the off-diagonal components of the property tensor. With the effective material property tensor, the coarse grid operator on each level can be easily obtained by the equation below:

$$L_H = \int_{\Omega} \nabla \phi_i^{H} \hat{\alpha}^{H} \nabla \phi_j^{H} d\Omega \quad (11)$$

where, $\phi_i^H$ and $\phi_j^H$ are test functions on each coarse grid, $\hat{\alpha}^H$ is the effective material property tensor of coarse grid.

The weak point of Sviercoski et al. [27] is that, the CP bounds on each coarse grid has to be computed from the finest grid. It requires too much computational time when using many grids. Therefore, one proposes to use different homogenization bounds which can be computed level by level.

Among different homogenization methods, the Voigt approximation is one of the best-known methods. In addition, it can be computed recursively. It is also referred to as the arithmetic mean:

$$\alpha_{xx}^{HV} = \frac{\sum_{k=1}^{k=N^h} \alpha_{xx}^{h}}{N^h} \quad (12)$$

where $\alpha_{xx}^{HV}$ is the diagonal component of $\hat{\alpha}^{HV}$ which is the average obtained by Voigt homogenization on coarse grid $l-1$, $\alpha_{xx}^h$ is the diagonal component of $\hat{\alpha}^h$ which is the material property on the fine grid $l$, $x = 1, 2, 3, N^h$ is the number of nodes on level $l$, which has the same volume as one element on level $l-1$.

Another approximation often used is the Reuss approximation. Equally, it can be obtained recursively, It is also known as the harmonic mean:

$$\alpha_{xx}^{HR} = \frac{N^h}{\sum_{k=1}^{k=N^h} \alpha_{xx}^{h-1}} \quad (13)$$

where $\alpha_{xx}^{HR}$ is the diagonal component of $\hat{\alpha}^{HR}$ which is the average obtained by Reuss homogenization at the grid $l-1$.

Instead of the CP bounds, the Voigt-Reuss (VR) bounds are used in this work. Thus the effective material property tensor can be obtained recursively. The assumption is that the effective value lies within the arithmetic and geometric averages of the VR bounds.

**Definition:** The effective material property tensor $\hat{\alpha}^H$ is the average of the arithmetic and the geometric average of the VR bounds. Which is presented as:

$$\alpha_{xx}^{H} = \frac{1}{2}(\alpha_{xx}^{a} + \alpha_{xx}^{g}) \quad (14)$$

where, $\alpha_{xx}^H$ is the diagonal component of $\hat{\alpha}^H$, $\alpha_{xx}^a$ is the diagonal value of the arithmetic average of the VR bounds which is defined as: $\alpha_{xx}^a = \frac{1}{2}(\alpha_{xx}^{HR} + \alpha_{xx}^{HV})$, $\alpha_{xx}^g$ is the diagonal value of the geometric average of the VR bounds which, is defined as: $\alpha_{xx}^g = \sqrt{\alpha_{xx}^{HR} \cdot \alpha_{xx}^{HV}}$.

The material used in this work is supposed to be isotropic, off-diagonal values of $\hat{\alpha}^H$ are zero. The strategy to obtain the off-diagonal components is therefore not presented in this work.

An analysis of different homogenization methods (e.g. Voigt, Reuss, Hashin-Shtrikman, Self-Consistent) has been performed, to confirm the strategy used in this work. The idea is to compute the effective thermal conductivity for a cubic structure, in which there is a spherical inclusion.

* The conductivity of the sphere is 100 W·m⁻¹·K⁻¹, 1 W·m⁻¹·K⁻¹ for the other part of the cube

* The edge length of the cube is 1 m
<figure><img src="image_4.png" /><figcaption>Fig. 2. Different homogenization methods. VR+: VR upper bound, VR -: VR lower bound, HS+: Hashin-Shtrikman upper bound, HS -: Hashin-Shtrikman lower bound, SC: Self-Consistent.</figcaption></figure>

• The radius of the sphere is between 0 m and 0.4 m

Fig. 2 illustrates the evaluation of the effective conductivity ob-
tained by different homogenization methods, when the sphere radius
varies. The effective property obtained by the VR bounds lies between
the Hashin–Shtrikman bounds. It confirms that the method used in this
work is robust and efficient.

Besides the coarse grid operators, prolongation and restriction op-
erators also need special treatment. The relation between the restriction
operators (i.e.R) and the prolongation operators (i.e.P) is:

$$
P = R^T \tag{15}
$$

The work of Boffy and Venner [4] presents the principle to derive *P* and *R*. The point is to consider the material discontinuities. The prolongation process will be briefly presented in this paper.

As illustrated in Fig. 3, the big box with solid edges represents one element on the coarse grid *l* − 1, and, the eight small boxes with dotted edges are the eight elements on the fine grid *l*. The temperature at each “black” node (e.g.A1, A2) of coarse grid is known. The goal is to obtain the temperature at all the 27 nodes of the fine grid. For the temperature at the eight “black” nodes of the fine grid coinciding with the ones on the coarse grid, one does an injection, which means:

$$
T_{A1}^{l} = T_{A1}^{l-1} \qquad (16)
$$

For the other nodes, instead of the linear prolongation of the stan-
dard MG method, the nodal material property is taken into account. It

<figure><img src="image_5.png" /><figcaption>Fig. 3. Prolongation scheme.</figcaption></figure>

reads as follows.

For the “red” nodes (Fig. 3) located at the center of the link between
two coarse grid nodes (e.g. node B1):

$$
T_{B1}^{l} = \frac{\alpha_{A1} T_{A1}^{l-1} + \alpha_{A2} T_{A2}^{l-1}}{\alpha_{A1} + \alpha_{A2}} \qquad (17)
$$

For the “blue” nodes which, are located at the center of each face of
the big box (e.g. node C1):

$$
T_{C1}^{l} = \frac{\alpha_{B1} T_{B1}^{l} + \alpha_{B2} T_{B2}^{l} + \alpha_{B3} T_{B3}^{l} + \alpha_{B4} T_{B4}^{l}}{\alpha_{B1} + \alpha_{B2} + \alpha_{B3} + \alpha_{B4}} \quad (18)
$$

For the “yellow” center node of the big box (e.g. node O):

$$
T_O^l = \frac{\alpha_{C1} T_{C1}^l + \alpha_{C2} T_{C2}^l + \alpha_{C3} T_{C3}^l + \alpha_{C4} T_{C4}^l + \alpha_{C5} T_{C5}^l + \alpha_{C6} T_{C6}^l}{\alpha_{C1} + \alpha_{C2} + \alpha_{C3} + \alpha_{C4} + \alpha_{C5} + \alpha_{C6}} \quad (19)
$$

With this strategy, the material property at each node weighs its
contribution for the prolongation. It is obvious that, if α is constant, it
becomes the linear prolongation of the standard MG method.

Finally, the application of the V-cycles MG method coupled with
homogenization technique can be illustrated by the following steps:

1. Compute Voigt and Reuss approximations on each level, re-
spectively; obtain the effective material properties ˜$\alpha^H$ for all le-
vels besides the finest level which, has the real material proper-
ties.

2. Carry out relaxations with the Jacobi solver on level *l*.

3. Transfer the temperature and restrict the residual to level
*l* − 1 and perform relaxations on this level.

4. Repeat steps 2, 3, 4 from the finest grid to the coarsest grid
$l = 1$.

5. Prolong the correction to the level *l* + 1 and relax on this
level.

6. Repeat step 5 until the finest level.

7. Loop step 2, 3, 4, 5, 6 until obtaining the required residual.

8. Output results.

To have a better performance, instead of using only V-Cycles, the
full MultiGrid (FMG) cycles are used in this work. Comparing to one V-
Cycle, one FMG cycle need more relaxations on coarse grids. However,
the computational cost on coarse grid is negligible compared to fine
grid. Further, FMG cycles provide for a good initialization solution on
the finest level, the number of relaxations on the finest grid is thus
lower than that of V-Cycles. The scheme of FMG cycles for 3 levels is
illustrated in Fig. 4. For all of the applications in this work, one starts
always by a 4 × 4 × 4 grid on level 1. The size of grids for level l + 1 is
two times smaller than that of level l, e.g. for a problem of 2048³ ele-
ments, one has 10 levels. v₀ is the number of relaxations performed on
level 1, v₁ is the number of relaxations performed on each level going
up. v₂ is the number of relaxations performed on each level going down.
For this FMG cycles, one uses ncy V-Cycles on each level. For the initial
solution of each fine level l + 1, one does a bi-linear interpolation of the
solution of level l.

Level

<figure><img src="image_6.png" /></figure>

FMG-Cycles

<caption>Fig. 4. FMG Cycles.</caption>
<table><caption>**Table 1**
Comparison between single level relaxation and a MG scheme.</caption><thead><tr><th></th><th>Single level</th><th>MG scheme</th></tr></thead><tbody><tr><td>Residual achieved</td><td>1.55 × 10<sup>-2</sup></td><td>7.89 × 10<sup>-6</sup></td></tr><tr><td>Number of relaxations on the finest level</td><td>4139</td><td>15</td></tr></tbody></table>

The performance of the MG scheme is studied below. One carries out a simulation of a spherical thermal inclusion with a conductivity of 10 Wm⁻¹·K⁻¹ for sphere and 1 Wm⁻¹·K⁻¹ for the other part of the cube. The radius is a quarter of the size of cube. The domain Ω is discretized with 128³ (more than two million) elements. The boundary conditions are set according to Eq. 3. The simulation is run on an office computer equipped with one processor “Intel(R) Core(TM)2 Quad CPU Q9650 @ 3.00 GHz”. For the single level iterative relaxation, the Jacobi relaxation is applied directly on this 128³ grid problem. For the MG scheme, one has 6 level for this 128³ grid problem, $v_0$, $v_1$, $v_2$, $n_{cy}$ are set to be 10, 1, 2 and 5, respectively. Three different values for the relaxation coefficient $\omega$, i.e. 0.5, 1.0 and 1.5, are used for both simulations.

Table 1 and Fig. 5 illustrate the performance of the MG scheme and the single level Jacobi solver. The convergence rate of the single level Jacobi solver decreased rapidly both for damped Jacobi with $\omega = 0.5$ and normal Jacobi with $\omega = 1.0$. On the other hand, the convergence rate of the MG scheme remains constant with $\omega = 0.5$ and $\omega = 1.0$. For the case of over relaxation i.e. $\omega = 1.5$, it diverges with both methods for this problem. With a $\omega = 1.0$, one has the best convergence performance for the MG scheme. After 4139 relaxations, single level Jacobi relaxation does not yet reach the initial solution of FMG cycles on the finest level, it confirms that one can have a good initial solution for the finest level with FMG cycles. If one ignores the computational cost on coarse grids, the MG scheme costs about 276, i.e. $\frac{4139}{15}$, times lower, with a residual 10 000 times smaller than that of single level Jacobi solver.

Another study is also carried out to compare the convergence performance between classical MG and the proposed method. Compared to the classical Geometric MultiGrid (GMG), the proposed method has two differences:

* Instead of the bi-linear prolongation and restriction operators, one proposes to consider material property for the prolongation and restriction operators.

* Homogenization technique is used in the proposed method to obtain the coarse grids material properties, however, in the classical MG, a simple average, i.e. Voigt approximation, is often used for the coarse grids material properties.

According to the difference between a classical MG and the

proposed MG, one proposes to carry out three simulations for the spherical case with a material property contrast of 1 000 (detailed in Section 4.2) and the composite case (detailed in section 4.4).

* **GMG-A:** Bi-linear restriction and prolongation operators, Voigt approximation material property for coarse grids

* **GMG-H:** Bi-linear restriction and prolongation operators, homogenized material property for coarse grids

* **Proposed method:** Considering material property for restriction and prolongation operators, homogenized material property for coarse grids

As illustrated in Fig. 6a and b, these three methods converge for both spherical and composite case. For spherical case, both the proposed method and GMG-H have the best performance. For the composite case, the proposed method has the best performance. For both cases, it shows that coarse grid material property has a large influence on the convergence speed, a representative material property for coarse grid is highly important to ensure good convergence. The idea to includes the material property for the prolongation and restriction operators does not always have a large improvement, for some symmetrical case, e.g. spherical inclusion, it does not have a large improvement compared to GMG-H. But for the complex case e.g. composite material, it has a good performance. As the aim of the proposed algorithm is to deal with complex materials with large property variations, the proposed algorithm is more efficient.

### 3.3. Parallel computing

As mentioned above, the goal of this work is the simulation of a domain discretized by more than eight billion elements, with a single processor, the computational time and memory are a big challenge.

These limitations lead one to consider High Performance Computing (HPC). The available machine is a supercomputer equipped with 252 nodes with 128 GB RAM for each node. Each node consists of two processors and each processor has 12 cores. Message Passing Interface (MPI) and OpenMP are the two most used methods for parallel computing. MPI is mainly for distributed memory machines, whereas, OpenMP is for machines with a shared memory. Using OpenMP only can not address our problem, as the memory of one node is not sufficient. MPI only can be a choice, but the difficulty is that, for the coarsest grid, one has only 4 × 4 × 4 elements. By consequence, for the coarsest level, one can uses only 64 MPI, which does not allow a good speedup (maximum 64). To take advantage of this machine, we uses a hybrid MPI/OpenMP parallel algorithm.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Convergence of the Jacobi solver (a) and MG scheme (b) on a 129³ nodes problem.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 6. Convergence of MG method with different intergrid operators.</figcaption></figure>

# 4. Results and applications

## 4.1. Validation

A multi-layer problem is used to validate the proposed solver. The domain $\Omega$: $1.0 \times 1.0 \times 1.0 \text{ cm}^3$ cube, consists of four uniform layers. The thermal conductivity of each layer is: $\alpha_1 = 1 \text{ W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}$, $\alpha_2 = 4 \text{ W}\cdot\text{m}^{-1}$, $\text{K}^{-1}$, $\alpha_3 = 8 \text{ W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}$ and $\alpha_4 = 4 \text{ W}\cdot\text{m}^{-1}\cdot\text{K}^{-1}$. The distribution of these four materials is presented in Fig. 7. Boundary conditions are applied as mentioned in Eq. 3. $T_1 = 0$ K and $T_2 = 1$ K are applied.

The analytical solution of this problem can be described as:

$$ T(z) = \begin{cases} \frac{32}{13}z & \text{for } 0 \le z < 0.25 \text{ cm} \\ \frac{8}{13} + \frac{8}{13}(z - 0.25) & \text{for } 0.25 \text{ cm} \le z < 0.5 \text{ cm} \\ \frac{10}{13} + \frac{4}{13}(z - 0.50) & \text{for } 0.5 \text{ cm} \le z < 0.75 \text{ cm} \\ \frac{11}{13} + \frac{8}{13}(z - 0.75) & \text{for } 0.75 \text{ cm} \le z \le 1.0 \text{ cm} \end{cases} \quad (20) $$

where $z$ is the coordinate of axis Z.

To validate the numerical solution, the error between the analytical and the numerical solution is computed. However, the FEM has a discretization error, the number of elements used to obtain the numerical solution affects the error between the analytical and the numerical solution. As a result, on one side, the error between the analytical and the numerical solution is analyzed; on the other side, the FEM discretization error is studied.

For the numerical solution, one discretizes $\Omega$ into $128 \times 128 \times 128$ (i.e. more than two million) cubic elements, applies the same boundary conditions as for the analytical solution and carries out the simulation. The coarsest level has $4 \times 4 \times 4$ grids, for each finer level, the grids size is devised by two i.e. 6 levels for a $128^3$ problem. Fig. 8 shows the temperature variations of the analytical and the numerical solution, respectively, along the Z direction. The temperature obtained by the numerical simulation is almost the same as the one obtained by the analytical solution. The $L_2$ error norm between the analytical solution and the numerical solution is 0.0027.

To analyze the influence of the element size, one discretizes $\Omega$ with $16 \times 16 \times 16$, $32 \times 32 \times 32$ and $64 \times 64 \times 64$ elements and computes the $L_2$ error norm compared to the analytical solution, respectively. Fig. 9 shows the $L_2$ error norm as a function of element size. The $L_2$ norm error decreases almost linearly in log-log scale with the element size.

The analytical and numerical solution show that the strategy of using the MG method coupled with homogenization technique can deal with problems with varying coefficients.

## 4.2. Spherical inclusion with large variation

The stability of the proposed method is analyzed in this subsection, when it handles a spherical inclusion problem with large material property variations. The parallel performance of our program is also studied for this application.

The domain $\Omega$ is a cube, which has two materials, as presented in

<figure><img src="image_9.png" /><figcaption>Fig. 7. Boundary conditions and conductivity.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 8. Temperature variation on Z direction of the analytical (Red) and the numerical (Blue) solution.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. $L_2$ norm error analysis.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Spherical inclusion.</figcaption></figure>

Fig. 10a. The radius of the inclusion is $R = L/4$, where, L is the size of cube. The thermal conductivity of the material in the sphere is 1 000 W·m⁻¹·K⁻¹, whereas it is set to 1 W·m⁻¹·K⁻¹ in the other part. The contrast between these two materials is 1 000. One discretizes the cube with 128 × 128 × 128 cubic elements. The boundary conditions is presented in Eq. 3. The coarsest starts always by a 4 × 4 × 4 grid. $v_0, v_1, v_2, n_{cy}$ are set to be 100, 4, 8 and 10, respectively.

Fig. 10b shows the temperature gradient in Ω. The large variation of

the conductivity on the interface explains the large variation of the temperature gradient around the interface.

This application confirms the good stability of this strategy in case of large variations of material properties.

### 4.3. Effective conductivity of cast iron

Cast iron is a well-known and widely used material in the industrial domain. The prediction of the conductivity of cast iron is a significant difficulty for researchers. Several papers investigate the conductivity, e.g. Helsing and Grimvall [16] regarded cast iron as a composite material and created a model to predict its conductivity. Nevertheless, since the distribution of carbon grains in cast iron affects its conductivity, the property of cast iron is different for different manufacturers. One proposes to use X-ray tomographic techniques to obtain the carbon grains distribution in an image format. The numerical simulation is then employed on this image to analyze the influence of carbon granules and to obtain the effective conductivity of cast iron.

The original tomographic image of cast iron is an image with 512 × 340 × 340 voxels [24]. The voxel size is 5.06µm. The region of interest (ROI) in this work is a part of this image. This part has 257 × 257 × 257 (more than 16 million) voxels. Each voxel is supposed to be one elementary node of the FEM discretization. A conductivity is assigned to each node. Fig. 11a illustrates the conductivity of the two components in cast iron, where, the black granules in this image are the carbon grains. The carbon conductivity is 129.0 W·m⁻¹·K⁻¹, for the other part, one takes the conductivity of iron which, is 80.4 W·m⁻¹·K⁻¹.

To obtain the effective conductivity of the cast iron, the homogenization method is used. As presented in the work of Özdemir et al. [22], the idea is to consider Ω to be one element, which is also referred to be a Representative Volume Element (RVE). The theory of RVE homogenization is briefly presented below.

The well known Fourier's law is described as:

$$- \mathbf{A} \cdot \nabla \theta = Q \qquad (21)$$

where, A is the effective thermal conductivity at macroscopic scale, θ is temperature at macroscopic scale, ∇θ is its gradient and Q is the total heat flux at macroscopic scale, which can be computed from the local heat flux:

$$Q = \frac{1}{V} \int_V q \, dv = - \frac{1}{V} \int_V \alpha \nabla T \, dv \qquad (22)$$

Three simulations are carried out with:

$$ \begin{cases} \nabla \theta_x = 1 \, \text{K} \cdot \text{m}^{-1} \\ \nabla \theta_y = 0 \\ \nabla \theta_z = 0 \\ T = \nabla \theta_x x \text{ on } \partial \Omega \end{cases} \quad \begin{cases} \nabla \theta_x = 0 \\ \nabla \theta_y = 1 \, \text{K} \cdot \text{m}^{-1} \\ \nabla \theta_z = 0 \\ T = \nabla \theta_y y \text{ on } \partial \Omega \end{cases} \quad \begin{cases} \nabla \theta_x = 0 \\ \nabla \theta_y = 0 \\ \nabla \theta_z = 1 \, \text{K} \cdot \text{m}^{-1} \\ T = \nabla \theta_z z \text{ on } \partial \Omega \end{cases} \qquad (23) $$

<figure><img src="image_13.png" /><figcaption>Fig. 11. Cast iron.</figcaption></figure>
respectively, where $\nabla\theta_x$, $\nabla\theta_y$ and $\nabla\theta_z$ are the temperature gradient in direction X, Y and Z, respectively. With each boundary condition, one column of A can be obtained.

For the FMG Cycles, 7 levels of grid are used, $v_0$, $v_1$, $v_2$, $n_{cy}$ are set to be 10, 2, 1 and 5, respectively. Fig. 11b shows the distribution of the temperature gradient in the case of $\nabla\theta_z = 1$, in which there are many inclusions. The location of the inclusions coincides with the location of the carbon grains.

The effective conductivity obtained is:

$$A = \begin{cases} 82.4311 & 0.0020 & 0.0040 \\ 0.0020 & 82.4223 & -0.0026 \\ 0.0040 & -0.0026 & 82.4277 \end{cases} W \cdot m^{-1} \cdot K^{-1}$$

Up to two significant digits, A can be described as:

$$A = \begin{cases} 82.43 & 0.00 & 0.00 \\ 0.00 & 82.42 & 0.00 \\ 0.00 & 0.00 & 82.43 \end{cases} W \cdot m^{-1} \cdot K^{-1}$$

which means that cast iron of this manufacturer is almost isotropic regardless of the random carbon distribution.

## 4.4. Effective conductivity of a layered composite material

Cast iron is almost isotropic, one may measure its conductivity experimentally. However, for layered composite materials, which are also widely used in the industrial domain due to its good performance, can be an anisotropic material. To carry out an experimental measurement, several external factors have to be observed, which is not simple and sometimes, not possible. Employing numerical simulations directly on tomographic images can be good alternative to know the composite properties.

The image used in this work is the image of a laminate composite material consisting of unidirectional E-glass fibers and a M9 epoxy matrix. It is a Glass Fiber Reinforced Polymer (GFRP) manufactured by the Hexcel Company. Its mechanical properties have been studied by Lecomte-Grosbras et al. [18]. The details of this image can be found in the work of Lecomte-Grosbras et al. [18]. In this work, the heat transfer in this GFRP is studied to obtain its effective conductivity.

The original image of this GFRP is an image consisting of 700 × 1300 × 1700 voxels, As mentioned in the work of Lecomte-Grosbras et al. [18], this material is designed with four layers, the orientation of fibers is + 15°, -15°, -15° and + 15°, respectively, for each layer. The idea is to take a cubic domain from the part which has the same fiber orientation. One takes 129 × 129 × 129 voxels from the part with a fiber orientation of - 15°, as the ROI (see Fig. 12). As presented in 12, the interface between the E-glass fiber and M9 epoxy matrix is not extraordinarily sharp. It is difficult to distinguish between these two phases (matrix and fiber). Instead of applying two discontinuous phases, one proposes to apply a continuous conductivity between 0.150 W·m⁻¹·K⁻¹ (epoxy) and 1.30 W·m⁻¹·K⁻¹ (E-glass fiber). One chooses to smooth the image gray level before it is used to compute the local material property at each voxel. It can be described as:

$$\alpha = 0.575 \left( \left( 1 - e^{-\frac{|GL-160.5|}{20}} \right) \operatorname{sign}(GL - 160.5) + 1 \right) + 0.15 \quad (24)$$

where GL is the original value of each voxel obtained by X-Ray tomography, which is an integer between 0 and 255. Except for the problem of the allocation of the conductivity, another problem is that the diameter of fiber is too small to have enough voxels in it. Sub-sampling i.e. linear interpolation, is therefore applied to this ROI to have more voxels in each fiber. The FEM discretization error therefore needs to be analyzed, to obtain the number of voxels needed for each section. A simulation with $\nabla\theta_z = 1 W\cdot m^{-1}\cdot K^{-1}$, $\nabla\theta_x = \nabla\theta_y = 0$ and $T = \nabla\theta_z$ on $\partial\Omega$ is performed. One time sub-sampling (case I) and two times sub-sampling (case II) are applied to the ROI, respectively. Fig. 13 illustrates the conductivity of each node in this ROI after one time sub-sampling. For

<figure><img src="image_14.png" /><figcaption>Fig. 12. ROI of the GFRP.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. GFRP conductivity.</figcaption></figure>

the FMG Cycles, $v_0$, $v_1$, $v_2$, $n_{cy}$ are set to be 10, 2, 1 and 5, respectively.

The third column of the effective property tensor $A_e$ is computed for each case.

For case I (7 levels i.e.257³ nodes):

$$A_e^3 = \{-0.001559 \quad 0.025400 \quad 0.744922\}$$

For case II (8 levels i.e.513³ nodes):

$$A_e^3 = \{-0.001607 \quad 0.026223 \quad 0.745158\}$$

It means that about up to three significant digits, the third column of the effective conductivity tensor is the same, or rather, one can take three significant digits for the $A_e$ obtained by one time sub-sampling, which is sufficient for industrial applications. The temperature gradient is also computed, as presented in Figure 14.

Similar to the previous cast iron application, other two simulations with boundary conditions of Eq. (23) are performed. The effective
<figure><img src="image_16.png" /><figcaption>Fig. 14. Temperature gradient of E-glass fibers in composite.</figcaption></figure>

conductivity of the ROI of the GFRP is:

$$A_c = \begin{cases} 0.625386 & 0.002162 - 0.001559 \\ 0.002162 & 0.628834 - 0.025400 \\ -0.001559 & 0.025400 - 0.744922 \end{cases} \text{W/(mK)}$$

With up to three significant digits, it reads:

$$A_c = \begin{cases} 0.625 & 0.002 - 0.002 \\ 0.002 & 0.629 - 0.025 \\ -0.002 & 0.025 - 0.745 \end{cases} \text{W/(mK)}$$

which confirms that GFRP is an orthotropic material.

The effective property tensor obtained above, is for the fibers with an orientation of -15°, for that of the +15° orientation, one can derive it directly.

## 4.5. Large simulation from a X-ray tomographic image and HPC performance analysis

The applications introduced above reveal that, the effective conductivity can be obtained by numerical simulation directly from an X-ray tomographic image, without any human intervention. The current tomographic images have 2048 × 2048 × 2048 voxels or more than 8 billion elements. The final application for this work it to carry out the numerical simulation with such a large image.

The image used in this case is the GFRP image of the previous application. One takes a part from the original image, the ROI consists of 513 × 513 × 513 voxels. As presented in Fig. 15, it consists of four layers with different E-glass fiber orientations. One employs a two times subsampling to obtain an image consisting of 2048 × 2048 × 2048 elements. The smoothing process on gray level is also applied and the material property has been assigned to each node as presented in Fig. 15. The same boundary conditions are applied as for the spherical inclusion case.

For the FMG cycles in this simulation, 10 levels of grids are used, $v_0$, $v_1$, $v_2$, $n_{cy}$ are set to be 10, 2, 1 and 5, respectively. 768 cores (64 MPI, 12 OpenMP/MPI) are used simultaneously. The calculation time is about 3.16 h. Figure 16 illustrates the residual evolution with the number of V-Cycles on level 10. Regardless of the size of the problem, the convergence remains very good. To achieve a residual of $10^{-6}$, only 5 V-Cycles on the finest level are needed. It means that the number of relaxations on the finest level is only 15. It confirms the efficiency of the

<figure><img src="image_17.png" /><figcaption>Fig. 15. E-glass fiber orientation in each layer and conductivity at each elementary node.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16. Convergence.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Temperature gradient.</figcaption></figure>

strategy used in this work. The temperature gradient is presented in Fig. 17. Fig. 17 and 18 illustrate the correspondence between conductivity and temperature gradient.

Since the aim of this work is the computation from large tomographic images, the parallel computing performance is analyzed with a
<figure><img src="image_20.png" /><figcaption>Fig. 18. Conductivity (Left) and temperature gradient (Right).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19. HPC performance.</figcaption></figure>

<table><caption>**Table 2**
Comparison between different configurations.</caption><thead><tr><td>Configurations</td><td>Time/s</td></tr></thead><tbody><tr><td>32 MPI, 12 threads</td><td>480</td></tr><tr><td>48 MPI, 8 threads</td><td>726</td></tr><tr><td>64 MPI, 6 threads</td><td>463</td></tr></tbody></table>

problem of 1024³ elements. A one-time sub-sampling is applied on the 512³ ROI to obtain an image of 1024³ elements. The aim is to analyze the parallel performance, so, instead of doing FMG cycles until final convergence, one proposes to do only one V-cycle on the finest grid.

For the performance of hybrid programming, there are two parameters to investigate: the number of MPI used and the number of threads per MPI task. As mentioned before, each processor has 12 cores for the available machine, each node has 2 processors. In order to have a good efficiency, the number of threads per MPI is limited by 12 to avoid the use OpenMP between two sockets, since OpenMP suffers from poor data access patterns when using two sockets. The number of cores maximum that can be used is 1 000 fixed by the owner.

As presented in Fig. 19, one simulation with 1 core and 1 MPI has been carried out to have a reference to compute the speedup. The number of MPI is set to 1,2,4,8,16,32 and 64; 12 threads are used in each MPI task. This curve illustrates that for a 1024³ problem, even with 768 cores, good speedup is obtained, within this number of cores, the speedup increases linearly as the number of cores increasing. Further, the speedup is close to be optimal, the speedup rate being close to 0.8.

Besides the number of cores that has a big influence on the parallel performance, the number of MPI and the number of OpenMP (for the same number of cores) can also have an influence: i.e. for 384 cores, one has the following configurations: 32 MPI with 12 threads, 48 MPI with

8 threads, 64 MPI with 6 threads, 96 MPI with 4 threads, 128 MPI with 2 threads.

Simulations for three different configurations, i.e. 32 MPI with 12 threads, 48 MPI with 8 threads and 64 MPI with 6 threads, are carried out for a problem for 1024³ elements. As presented in Table 2, with 48 MPI and 8 threads per MPI, one has a poor performance since in one node, there are 3 MPI which means there is at least one MPI taking cores from two sockets. With 64 MPI and 6 threads per MPI, one obtains a better performance than that with 32 MPI and 12 threads for each, but the difference is only $\frac{480-463}{463} \approx 3.67\%$. It confirms that nevertheless the current program does not allow to use as more as MPI tasks that we want, but with 64 MPI and 12 cores, a sufficient performance is obtained. For some large problems like 2024³ grids, using 12 cores for each of the 64 MPI, i.e. 768 cores in total, one can already finish the simulation in about 3 h.

## 5. Discussion and conclusions

The goal of this paper is to show that one can employ numerical simulations directly on tomographic images. To perform the simulations with such a large number of elements and such large variations of materials properties requires dedicated algorithms and hardware. The strategy to use the MG method coupled with a homogenization technique, permits one to deal with this kind of problems. The applications and numerical comparison presented above demonstrate the efficiency of the MG method. The homogenization technique shows its capacity to increase the convergence performance of the MG scheme, when large variations of materials properties exits. The Matrix Free FEM demonstrate its good performance for large problems up to 8 billion of elements. The strategy to apply one voxel per elementary node avoids human intervention. The effective material property can be automatically obtained by using the large X-ray tomographic image, as an input, without complex experimental measurement. The Hybrid MPI/OpenMP programming shows its good feasibility and performance for the MG method.

The thermal conductivity is analyzed in this work. In future work, the mechanical property of materials will be analyzed. The property at each node is supposed to be isotropic, additional research is needed for the anisotropic case. For the parallel computing of MG method, the re-partitioning and load balancing need to be applied for the future work if necessary to have the best performance.

## Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

## CRediT authorship contribution statement

**Xiaodong Liu:** Writing - original draft, Software, Visualization, Methodology, Investigation. **Julien Réthoré:** Funding acquisition, Supervision, Validation, Methodology, Investigation. **Marie-Christine Baietto:** Writing - review & editing. **Philippe Sainsot:** Writing - review & editing. **Antonius Adrianus Lubrecht:** Writing - review & editing. Supervision, Validation, Methodology, Investigation.

## Acknowledgments

The support of Région Pays de la Loire and Nantes Métropole through a grant Connect Talent IDS is gratefully acknowledged.

## References

[1] R. Alcouffe, A. Brandt, J. Dendy Jr, J. Painter, The multi-grid method for the diffusion equation with strongly discontinuous coefficients, SIAM J. Sci. Stat. Comput. 2 (1981) 430-454.
[2] M. Bessho, I. Ohnishi, J. Matsuyama, T. Matsumoto, K. Imai, K. Nakamura, Prediction of strength and strain of the proximal femur by a ct-based finite element method, J. Biomech. 40 (2007) 1745-1753.

[3] N. Biboulet, A. Gravouil, D. Dureisseix, A. Lubrecht, A. Combescure, An efficient linear elastic fem solver using automatic local grid refinement and accuracy control, Finite Elem. Anal. Des. 68 (2013) 28-38.

[4] H. Boffy, C.H. Venner, Multigrid solution of the 3d stress field in strongly heterogeneous materials, Tribol. Int. 74 (2014) 121-129.

[5] A. Brandt, Multi-level adaptive solutions to boundary-value problems, Math. Comput. 31 (1977) 333-390.

[6] Brandt, A., Livne, O.E., 2011. Multigrid techniques: 1984 guide with applications to fluid dynamics, vol. 67, SIAM.

[7] G.F. Carey, B.N. Jiang, Element-by-element linear and nonlinear solution schemes, Commun. Appl. Numer. Methods 2 (1986) 145-153.

[8] H. Cheng, S. Torquato, Effective conductivity of periodic arrays of spheres with interfacial resistance, Proc. R. Soc. London A: Math. Phys. Eng. Sci. (1997) 145-161. The Royal Society.

[9] B. Engquist, E. Luo, New coarse grid operators for highly oscillatory coefficient elliptic problems, J. Comput. Phys. 129 (1996) 296-306.

[10] B. Engquist, E. Luo, Convergence of a multigrid method for elliptic equations with highly oscillatory coefficients, SIAM J. Numer. Anal. 34 (1997) 2254-2273.

[11] M. Ferrant, S.K. Warfield, C.R. Guttmann, R.V. Mulkern, F.A. Jolesz, R. Kikinis, 3d image matching using a finite element based elastic deformation model, International Conference on Medical Image Computing and Computer-Assisted Intervention, Springer. (1999) 202-209.

[12] C. Flaig, P. Arbenz, A highly scalable matrix-free multigrid solver for µfe analysis based on a pointer-less octree, International Conference on Large-Scale Scientific Computing, Springer, 2011, pp. 498-506.

[13] B.P. Flannery, H.W. Deckman, W.G. Roberge, K.L. D'AMICO, Three-dimensional x-ray microtomography, Science 237 (1987) 1439-1444.

[14] C. Gruescu, A. Giraud, F. Homand, D. Kondo, D. Do, Effective thermal conductivity of partially saturated porous rocks, Int. J. Solids Struct. 44 (2007) 811-833.

[15] H. Gu, J. Réthoré, M.C. Baietto, P. Sainsot, P. Lecomte-Grosbras, C.H. Venner, A.A. Lubrecht, An efficient multigrid solver for the 3d simulation of composite materials, Comput. Mater. Sci. 112 (2016) 230-237.

[16] J. Helsing, G. Grimvall, Thermal conductivity of cast iron: models and analysis of experiments, J. Appl. Phys. 70 (1991) 1198-1206.

[17] R. Hoekema, K. Venner, J.J. Struijk, J. Holsheimer, Multigrid solution of the potential field in modeling electrical nerve stimulation, Comput. Biomed. Res. 31 (1998) 348-362.

[18] P. Lecomte-Grosbras, J. Réthoré, N. Limodin, J.F. Witz, M. Brieu, Three-dimensional investigation of free-edge effects in laminate composites using x-ray tomography and digital volume correlation, Exp. Mech. 55 (2015) 301-311.

[19] M. Lengsfeld, J. Schmitt, P. Alter, J. Kaminsky, R. Leppek, Comparison of geometry-based and ct voxel-based finite element modelling and experimental validation, Med. Eng. Phys. 20 (1998) 515-522.

[20] N. Michailidis, F. Stergioudi, H. Omar, D. Tsipas, An image-based reconstruction of the 3d geometry of an al open-cell foam and fem modeling of the material response, Mech. Mater. 42 (2010) 142-147.

[21] T.T. Nguyen, J. Rethore, J. Yvonnet, M.C. Baietto, Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials, Comput. Mech. 60 (2017) 289-314.

[22] I. Özdemir, W. Brekelmans, M. Geers, Computational homogenization for heat conduction in heterogeneous solids, Int. J. Numer. Methods Eng. 73 (2008) 185-204.

[23] H. Proudhon, J. Li, F. Wang, A. Roos, V. Chiaruttini, S. Forest, 3d simulation of short fatigue crack propagation by finite element crystal plasticity and remeshing, Int. J. Fatigue 82 (2016) 238-246.

[24] J. Rannou, N. Limodin, J. Réthoré, A. Gravouil, W. Ludwig, M.C. Baietto-Dubourg, J.Y. Buffière, A. Combescure, F. Hild, S. Roux, Three dimensional experimental and numerical multiscale analysis of a fatigue crack, Comput. Methods Appl. Mech. Eng. 199 (2010) 1307-1325.

[25] M. Schneebeli, S.A. Sokratov, Tomography of temperature gradient metamorphism of snow and associated changes in heat conductivity, Hydrol. Process. 18 (2004) 3655-3665.

[26] Y.S. Song, J.R. Youn, Evaluation of effective thermal conductivity for carbon nanotube/polymer composites using control volume finite element method, Carbon 44 (2006) 710-717.

[27] R.F. Sviercoski, P. Popov, S. Margenov, An analytical coarse grid operator applied to a multiscale multigrid method, J. Comput. Appl. Math. 287 (2015) 207-219.

[28] T. Tezduyar, S. Aliabadi, M. Behr, A. Johnson, S. Mittal, Parallel finite-element computation of 3d flows, Computer 26 (1993) 27-36.