# An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using CT images

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

► To cite this version:

Xiaodong Liu, Julien Réthoré, Marie-Christine Baietto, Philippe Sainsot, Antonius Lubrecht. An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using CT images. Computational Mechanics, 2020, 10.1007/s00466-020-01909-y. hal-02919095

HAL Id: hal-02919095

https://hal.science/hal-02919095

Submitted on 21 Aug 2020

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An efficient finite element based multigrid method for simulations of the mechanical behavior of heterogeneous materials using CT images

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

¹GeM, CNRS UMR 6183, Centrale Nantes, F-44321, France

²Univ Lyon, INSA-Lyon, CNRS UMR5259, LaMCoS, F-69621, France

19 August, 2020

## Abstract

X-ray tomography techniques give researchers the full access to material inner structures. With such ample information, employing numerical simulation on real material images becomes more and more common. Material behavior, especially, of heterogeneous materials, e.g. polycrystalline, composite, can be observed at the microscopic scale with the computed tomography (CT) techniques. In this work, an efficient strategy is proposed to carry out simulations on large 3D CT images. A proposed matrix free type finite element based MultiGrid method is applied to improve convergence speed and to reduce memory space requirements. Homogenization techniques are used to obtain specific operators to enhance the convergence of the MultiGrid method when large material property variations are present. Hybrid parallel computing is implemented for memory space and computing time reasons. Free edge effects in a laminate composite with more than 16 billion degrees of freedom and crack opening in a cast iron are studied using the proposed strategy.

**Keywords**— Heterogeneous materials – Finite element method – MultiGrid – Homogenization – Mechanical behavior – Tomography images

## 1 Introduction

One of the subjects of material science is to understand material behavior and then to improve material performance to respond to industrial requirements. However, to manufacture materials with such good performance, is not an easy job. Many defects are present in materials, e.g. voids, cracks or broken fibers in composite materials. CT images of such complex materials permit one to evaluate the presence of these defects. Nevertheless, the aim of researchers is to understand how these defects affect material performance. Numerical simulations are often used to analyze material performance. The main idea of this work is to carry out the calculation directly on CT images. This permits one to observe how the material structure affects its performance at the microscopic scale. However, these defects make the simulation expensive or even impractical.

<footnote>*Corresponding author: julien_rethore@ec-nantes.fr</footnote>
One of the most used methods is to homogenize materials, which means instead of focusing on the local defects, the material is supposed to be homogenized at a macroscopic scale. This technique is called homogenization. Hashin (1979) explains how to predict the macroscopic property of a fiber composite from its microscopic structure. Özdemir et al. (2008) presents the application of a homogenization method for the thermal conduction of heterogeneous materials.

However, despite many investigations in the homogenization, the drawback of this technique is that it can not represent microscopic details. One loses physical phenomena related to the local scale. Thus, CT simulations at the microscopic scale become essential to overcome these drawbacks of homogenization methods.

Finite element methods (FEM) are often used to carry out CT simulations. Lengsfeld et al. (1998) and Bessho et al. (2007) present an early and a recent CT images application on the human bone problem by FEM, e.g. hip fractures of the human femur. Ferrant et al. (1999), Michailidis et al. (2010) and Proudhon et al. (2016) analyze properties of industrial materials by applying FEM simulations on CT images. Another well-known method is the fast fourier transform (FFT). One of the most typical examples is presented in Moulinec and Suquet (1995), which introduces a FFT-based numerical method to simulate the mechanical properties of composite materials. The finite difference method (FDM) is also often used to carry out CT simulations. Gu et al. (2016) presents the application of the FDM with the MultiGrid method on simulations of 3D composite materials. All these three well-known methods perform well for certain cases. However, for FEM applications on materials with complex structures, e.g irregular grains in polycrystalline materials, the meshing step requires human interventions, which is very time consuming. The work in Nguyen et al. (2017) presents the application of FEM on crack propagation for polycrystalline materials, the meshing step for this kind of structure is very complex. For FFT-based numerical methods, the limitation is that it is not robust for no-periodic problems. Gu et al. (2016) is a good application of a FDM, however, it was developed only for small problems, furthermore, the implementation of the boundary condition in FDM is intricate.

Motivated by the limitations of the current methods, we propose an efficient and automatic strategy to perform numerical simulations directly from large scale 3D tomography images. A proposed matrix free type finite element method is used in this work. The idea of matrix free-FEM (MF-FEM) was first introduced in Carey and Jiang (1986) in 1980s. It was investigated by Augarde et al. (2006); Hughes et al. (1987); van Rietbergen et al. (1996) during the last three decades. Nowadays, it is often used in parallel computing, see Tezduyar et al. (1993). The MF-FEM was proposed to save memory space and to apply parallel computing. It is particularly favorable for problems with regular and few element types. It does a good job for problems arising from voxel conversions, where the mesh generation can be automatic. In this case, all the generated elements have exactly the same geometry and orientation. Arbenz et al. (2008) presents a human bone structures simulation using CT images by the MF-FEM. The main drawback of this technique is that it is difficult to find an efficient preconditioner. An often used preconditioner is the element-by-element (EBE) preconditioning (Hughes et al. (1987)) or the Jacobi preconditioning (diagonal scaling). In this work, a modification is proposed to the typical MF-FEM. We propose to not use any matrix but to compute the pre-conditioner degree of freedom (DoF) by DoF (see in Section 2).

However, the MF-FEM does not sufficiently efficient for complex problems like CT simulations. An often used accelerator is the well-known MultiGrid algorithms, as presented in Arbenz et al. (2008). The MF-FEM coupled with the aggregation-based MultiGrid method. However, their objective was to perform numerical simulations on complex structure, i.e. human bone, with homogeneous material. In this work, the material heterogeneity is accounted for. In recent years, Kronbichler and Ljungkvist (2019) studied the performance of MultiGrid methods for MF high-order FEM on GPUs.

The MultiGrid (MG) method enhances convergence of the FEM and the FDM, e.g. see Liu et al. (2019) for the efficiency of MG algorithms applied on the FEM. MG methods were first
introduced in the 1960s, and then, enriched by Brandt Brandt (1973, 1977) in the 1970s. It is often used with the FDM. Work of Venner and Lubrecht Venner and Lubrecht (2000) developed an efficient and robust FDM based MG solver for the elasto-hydrodynamic lubrication problem. The work presented in Watremetz et al. (2007) introduced the application of MG methods for a 2D graded problem. Boffy et al. Boffy et al. (2012); Boffy and Venner (2014) employed MG FDM for small heterogeneous material mechanical problems. The work of Gu et al. Gu et al. (2016) investigates the FDM based MG methods for the application of composite material properties. This work applies the MG algorithm to achieve a good convergence speed.

The aim of this work is to use a tomographic image of a real material as input, and to carry out a numerical simulation with linear elastic equations on this image to observe the material behavior under a certain loading. Due to the large size of the problem, parallel computing is used. However, the detailed analysis of the performance, in terms of memory and computational cost, is beyond the scope of the present paper which focuses on the ability of the proposed methodology to deal with real 3D images of heterogeneous materials.

The outline of this paper is: Section 2 reviews the fundamentals of linear elasticity and its discretization by the proposed MF-FEM; Section 3 introduces the strategy in detail, both the MG method with special intergrid operators and parallel computing are applied to improve efficiency; Section 4 presents the validation of proposed method and some applications; Section 5 provides discussions and conclusions.

# 2 Problem statement and theory

To solve a mechanical problem, the equations of equilibrium are the basics. The fundamentals of solid mechanic equations and their numerical solutions are presented in this section.

## 2.1 Governing Equations

Assuming a deformed domain $\Omega$ as presented in Figure 1, its boundary is denoted by $\partial\Omega$. Its exterior normal is $\vec{n}$. The second order stress tensor is described by $\sigma$. The body force in $\Omega$, which is due to gravity, magnetism, etc., is $\vec{f}$. The acceleration is presented by $\vec{a}$. The density is denoted by $\rho$. The equations of equilibrium can be described as:

$$ \nabla \cdot \sigma + \rho \vec{f} = \rho \vec{a} \qquad (1) $$

For the sake of simplicity, inertia effects are neglected in this paper. The equation of equilibrium can be simplified as:

$$ \nabla \cdot \sigma = \vec{0} \qquad (2) $$

Besides the equilibrium equations, the boundary conditions on $\partial\Omega$ are essential to solve a mechanical problem. Two often used boundary conditions are Dirichlet and Neumann boundary conditions. The Dirichlet boundary condition applied on $\partial\Omega_D$ is often presented as prescribed displacements, e.g. $\vec{u} = \vec{u}_0$ on $\Gamma_1$. The Neumann boundary condition applied on $\partial\Omega_N$ refers to the external force $\vec{f}_{ex}$, which can be described as:

$$ \vec{f}_{ex} = \sigma \cdot \vec{n} $$

e.g. applying $\vec{f}_0$ on $\Gamma_2$.

To solve equation (2), as stated at the beginning of this paper, one has many choices e.g. the FDM, the FEM, the FFT, etc.. The aim of this work is to build an automatic solver using information directly from CT images. As presented in the above section, the matrix free finite element method (MF-FEM) is suitable for the CT simulations. Instead of taking the heavy meshing step, each voxel in the CT image is supposed to be an elementary node. 8-node cuboid
<figure><img src="image_2.png" /><figcaption>Figure 1: An elastic body $\Omega$ and its boundary</figcaption></figure>

elements are therefore used to discretize $\Omega$. With this strategy, the mesh generation step becomes automatic.

Multiplying equation (2) by a test function $\vec{u}^*$ and integrating on $\Omega$, the weak form of equation (2) reads:

$$ \int_{\Omega} \nabla \cdot \sigma \vec{u}^* d\Omega = 0 \qquad (3) $$

Applying the divergence theorem and integrating by parts, it reads:

$$ \int_{\Omega} \sigma \nabla \vec{u}^* d\Omega = \int_{\partial\Omega} \sigma \cdot \vec{n} \vec{u}^* dS \qquad (4) $$

which can be described as the equilibrium of internal force $\vec{f}_{in}$ and external force $\vec{f}_{ex}$:

$$ \vec{f}_{in} = \vec{f}_{ex} \qquad (5) $$

with

$$ \left\{ \begin{array}{l} \vec{f}_{in} = \int_{\Omega} \sigma \nabla \vec{u}^* d\Omega \\ \vec{f}_{ex} = \int_{\partial\Omega} \sigma \cdot \vec{n} \vec{u}^* dS \end{array} \right. \qquad (6) $$

Applying finite element discretization, the displacements $\vec{u}$ can be described by:

$$ \vec{u}^c \approx \sum_{i=1}^{N} \psi_i^c \hat{u}_i^c \qquad (7) $$

where $i$ denotes the node number, and $N$ is the number of nodes. The index $c$ represents three directions with $c = \{1, 2, 3\}$. The displacement at the elementary node $i$ is $\vec{u}_i$. The component $c$ of the displacement at the node $i$ is described by $\hat{u}_i^c$. The shape function $\psi$ is the typical shape function of a 8-node cubic element. The test function $\vec{u}^*$ reads:

$$ \vec{u}^* = \psi_j^c $$

where $j$ is the node number.

The internal force for node $j$ and component $c$ reads:

$$ f_{in}^c = \sum_m \sum_e \sum_{g=1}^{8} w_g \sigma_{cm}(\vec{u}) \nabla_m \psi_j^c \qquad (8) $$
with $m = \{1, 2, 3\}$. The Gauss integration point number is $g$, and $w_g$ is its weight. 8 Gauss integration points are used in each element, and $e$ is the element number. The relationship between $\sigma$ and $\vec{u}$ is implied by $\sigma(\vec{u})$, which can be obtained by the constitutive law (The constitutive law of a linear elastic material is presented in the next part).

To solve equation (4), the typical FEM process is to compute the stiffness matrix and to use direct or iterative solvers. However, a typical FEM process is almost impossible for a large scale CT image simulation because of memory space limitations, e.g. for a problem with more than 18 billion DoF. In this work, a modified MF-FEM is used.

The often used matrix free methods are: the power method, the conjugate gradient method, the locally optimal block preconditioned conjugate gradient method and etc.. However, in these methods, small local matrix are presented. We propose a new matrix free type FEM to not reserve any local matrix, but to compute the unknowns node-by-node (DoF-by-Dof) as presented in Liu et al. (2019) for the thermal conduction problem.

In this work, it can be described by equation (9).

$$ \hat{i}^{\vec{t}e+1}_{\vec{u}} = \hat{i}^{\vec{t}e}_{\vec{u}} + \omega \mathbf{S}(\vec{f}_{ex} - \vec{f}_{in}) \quad (9) $$

where $\hat{i}^{\vec{t}e+1}_{\vec{u}}$ and $\hat{i}^{\vec{t}e}_{\vec{u}}$ are the displacements of one node in the current iteration and previous iteration. The relaxation coefficient is denoted by $\omega$, for $0 < \omega < 1$, it refers to damped Jacobi, i.e. under-relaxation, which is used in this work. The diagonal matrix $\mathbf{S}$ is the compliance of the diagonal value of the stiffness matrix at this node. It reads:

$$ \mathbf{S} = \begin{bmatrix} \frac{1}{K_u} & 0 & 0 \\ 0 & \frac{1}{K_v} & 0 \\ 0 & 0 & \frac{1}{K_w} \end{bmatrix} \quad (10) $$

where $K_u$, $K_v$ and $K_w$ are the diagonal values of the stiffness matrix at each DoF. To obtain $\mathbf{S}$, the stiffness matrix is not necessary, only its diagonal terms are used.

With the proposed MF-FEM iterative solver, the large scale problem can be solved node by node.

## 2.2 The linear elasticity

In this work, the above strategy is applied to a linear elastic problem.

For a linear elastic problem, the constitutive equation reads:

$$ \boldsymbol{\sigma} = K \operatorname{tr}(\boldsymbol{\varepsilon}) \boldsymbol{I} + 2G \left( \boldsymbol{\varepsilon} - \frac{1}{3} \operatorname{tr}(\boldsymbol{\varepsilon}) \right) \quad (11) $$

where $\boldsymbol{\varepsilon}$ is the second order strain tensor, $K$ and $G$ refer to the bulk modulus and the shear modulus, respectively. The unit tensor is denoted by $\mathbf{I}$. Using with the Young's modulus $E$ and the Poisson ratio $\nu$, $K$ and $G$ read:

$$ K = \frac{E}{3(1 - 2\nu)}, \quad G = \frac{E}{2(1 + \nu)} $$

Strain $\boldsymbol{\varepsilon}$ can be described by displacement $\vec{u}$:

$$ \boldsymbol{\varepsilon} = \frac{1}{2} (\nabla \vec{u} + \nabla^T \vec{u}) \quad (12) $$

Its trace being defined as

$$ \operatorname{tr}(\boldsymbol{\varepsilon}) = \operatorname{tr} \left( \frac{1}{2} (\nabla \vec{u} + \nabla^T \vec{u}) \right) = \operatorname{tr}(\nabla \vec{u}) $$
Combining equation (11) and (12), equation (11) reads:

$$ \sigma = \left(K - \frac{2G}{3}\right) \mathrm{tr}(\nabla \vec{u}) \mathbf{I} + G(\nabla \vec{u} + \nabla^T \vec{u}) \quad (13) $$

Applying the finite element discretization with equation (7), the stress reads:

$$ \sigma_{cm}(\vec{u}) = \sum_{i} \left[ \left(K - \frac{2G}{3}\right) \left(\sum_{k} \nabla_{k} \psi_{i}^{k} \hat{u}_{i}^{k}\right) \delta_{cm} + G \left(\nabla_{m} \psi_{i}^{c} \hat{u}_{i}^{c} + \nabla_{c} \psi_{i}^{m} \hat{u}_{i}^{m}\right) \right] \quad (14) $$

where $k = \{1, 2, 3\}$. $\delta_{cm}$ is the Kronecker delta. It reads:

$$ \delta_{cm} = \begin{cases} 1 & c = m \\ 0 & c \neq m \end{cases} $$

Combining equation (8) and (14), one obtains the component $c$ of the internal force on node $j$:

$$ f_{in}^{c} = \sum_{e} \sum_{g=1}^{8} \sum_{i} w_{g} \left\{ \left(K^{g} - \frac{2G^{g}}{3}\right) \left(\sum_{k} \nabla_{k} \psi_{i}^{k} \hat{u}_{i}^{k}\right) \nabla_{c} \psi_{j}^{c} + \sum_{m} \left[G^{g} \left(\nabla_{m} \psi_{i}^{c} \hat{u}_{i}^{c} + \nabla_{c} \psi_{i}^{m} \hat{u}_{i}^{m}\right) \nabla_{m} \psi_{j}^{c}\right] \right\} \quad (15) $$

where $K^g$ and $G^g$ are the material property at a Gauss point. As mentioned before, the material property is assigned to each node from voxel information. An interpolation is used to obtain the material property at each Gauss point. It reads:

$$ K^g = \sum_{i=1}^{8} \psi_i K_i \quad G^g = \sum_{i=1}^{8} \psi_i G_i. $$

where $K_u$ is the first diagonal term of $\mathbf{S}$ in equation (10). It reads:

$$ K_u = \sum_{e} \sum_{g} w_g \left[ \left(K^g - \frac{2G^g}{3}\right) \nabla_c \psi_j^c \nabla_c \psi_j^c + \sum_m G^g \nabla_m \psi_j^c \nabla_m \psi_j^c + G^g \nabla_c \psi_j^c \nabla_c \psi_j^c \right] \quad (16) $$

where $c = 1$. $K_v$ and $K_w$ are computed with $c = 2$ and $c = 3$, respectively. Equation (9) can thus be solved DoF by DoF.

## 3 MultiGrid

As presented in Liu et al. (2019), using only the MF-FEM Jacobi solver, the convergence rate can be very slow after several iterations. The MG method is therefore applied to the iterative solver to improve convergence.

The principle of the MG method is that relaxations on the fine grid can eliminate high frequency errors, the low frequency errors can be eliminated by relaxing on the coarse grid. The first step of MG is therefore to construct several levels of grid. In this work, the grid size on level $l$ is two times larger than that on level $l+1$, e.g, for a 128³-grid problem with 3 levels, the number of points on level is: 128³ on level 3, 64³ on level 2 and 32³ on level 1.
Once different grids are constructed, three important operators are needed to employ a MG scheme: the coarse grid operator, the restriction operator **R** and the prolongation operator **P**. As presented in Liu et al. (2019), the standard MG scheme can not deal efficiently with problems with large variations of the material properties. Special intergrid operators can be generated as proposed in Liu et al. (2019).

The process to obtain the $\vec{F}_{in}$ and the **S** indicates that for a representative coarse grid operator, one needs the material property and the shape function on the coarse grid. The shape function of an element on each coarse grid is simple to obtain with the finite element theorem. For the material property on the coarse grid, the standard method is to compute the arithmetic average, i.e Voigt approximation, of the material property on the coarse grid. However, for heterogeneous materials with large property variations, a simple average can not represent the material on the fine grid. Many researchers have tried to apply MG algorithms for problems with large material property variations. The first work to improve MG convergence for problems with large variations is the work of Alcouffe et al. Alcouffe et al. (1981) for a 2D problem. Hoekema et al. Hoekema et al. (1998) extended it to 3D. This work proposes to apply the homogenization method proposed in Liu et al. (2019), which applies the homogenization method on MG algorithms to obtain the coarse grid operator for a thermal conduction problem.

The idea is to first compute the Voigt and Reuss approximations of **K** and **G** on all the coarse grids recursively as follow:

$$ \left\{ \begin{array}{ll} \text{Voigt: } K^{HV} = \frac{\sum_{1}^{N_h} K^h}{N_h}, & G^{HV} = \frac{\sum_{1}^{N_h} G^h}{N_h} \\ \text{Reuss: } K^{HR} = \frac{\sum_{1}^{N_h} \frac{1}{K^h}}{\sum_{1}^{N_h} \frac{1}{K^h}}, & G^{HR} = \frac{\sum_{1}^{N_h} \frac{1}{G^h}}{\sum_{1}^{N_h} \frac{1}{G^h}} \end{array} \right. $$

where $K^{HV}$, $G^{HV}$, $K^{HR}$ and $G^{HR}$ are Voigt and Reuss approximations of the bulk modulus and shear modulus on the coarse grid, The number of nodes on level $l$ is described by $N_h$, which has the same volume as one element on level $l-1$.

According to the work on homogenization techniques, the effective material property lies between the Voigt and Reuss approximations, which is called Voigt-Reuss (VR) bounds. VR bounds are not the most accurate bounds, but one can compute them recursively from the finest level. The effective material property on the coarse grid can be obtained by the following equation:

$$ \left\{ \begin{array}{l} K^H = \frac{\left( \frac{K^{HV}+K^{HR}}{2} + \sqrt{K^{HV}K^{HR}} \right)}{\frac{2}{2}} \\ G^H = \frac{\left( \frac{G^{HV}+G^{HR}}{2} + \sqrt{G^{HV}G^{HR}} \right)}{\frac{2}{2}} \end{array} \right. \qquad (17) $$

which means the effective $K^H$ and $G^H$ on the coarse grid is the mean of the arithmetic and geometric average of the VR bounds on this level. According to the comparison of different homogenization methods in Liu et al. (2019), the proposed homogenization scheme is robust and efficient.

With the correct material property and correct shape function on the coarse grid, one can finally define the coarse grid operator. For the restriction operator **R** and the prolongation operator **P**, the process of the prolongation of corrections can be described with one element on the coarse grid $l-1$ and its eight elements on the fine grid $l$ as presented in Figure 2. The subject is to bring the corrections, i.e displacement correction $\vec{e}$, from level $l-1$ to level $l$. The displacement corrections on level $l-1$ at black points is known. Instead of computing a simple average, one proposes to account for the material property at each node. The bulk modulus **K** is used in this strategy.

For the black points on level $l$, one performs an injection from level $l-1$ to $l$, e.g. point A1:

$$ e_{A1}^{\vec{l}} = e_{A1}^{\vec{l}-1}. $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Prolongation scheme</figcaption></figure>

For red points, one computes them from the $\vec{e}$ at black point of level $l$, e.g. point B1:

$$e_{B1}^{\vec{l}} = \frac{K_{A1}e_{A1}^{\vec{l}} + K_{A2}e_{A2}^{\vec{l}}}{K_{A1} + K_{A2}}.$$

For blue points, one computes them from the $\vec{e}$ at red points on level $l$, e.g. point C1:

$$e_{C1}^{\vec{l}} = \frac{K_{B1}e_{B1}^{\vec{l}} + K_{B2}e_{B2}^{\vec{l}} + K_{B3}e_{B3}^{\vec{l}} + K_{B4}e_{B4}^{\vec{l}}}{K_{B1} + K_{B2} + K_{B3} + K_{B4}}.$$

For yellow points, one computes them from the $\vec{e}$ at blue points on level $l$, e.g. point O:

$$e_{O}^{\vec{l}} = \frac{1}{K_{C1} + K_{C2} + K_{C3} + K_{C4} + K_{C5} + K_{C6}} \cdot (K_{C1}e_{C1}^{\vec{l}} + K_{C2}e_{C2}^{\vec{l}} + K_{C3}e_{C3}^{\vec{l}} + K_{C4}e_{C4}^{\vec{l}} + K_{C5}e_{C5}^{\vec{l}} + K_{C6}e_{C6}^{\vec{l}}).$$

For the restriction process, one can use $R = P^T$ to obtain the restriction operator.

After constructing several levels and its operators to go up and go down, one has to choose a MG cycle. The best-known MG cycles are V-cycles, W-cycles and FMG cycles. In this work, FMG cycles are performed to have a good initial solution on the finest grid. FMG cycles with 3 levels are illustrated in Figure 3. On level 1, one performs $\nu_0$ relaxations, $\nu_1$ relaxations are performed on each level going up. $\nu_2$ relaxations are performed on each level going down. $n_{cy}$ V-Cycles are used on each level. A linear interpolation of the solution of level $l$ is applied to obtain the initial solution of level $l + 1$.

With the above strategy, the construction of a matrix free FEM iterative solver based homogenized MG algorithms is finished. One can start performing numerical simulations. However, for a large image, e.g. with more than 18 billion DoF, the computational time is too long using a standard computer. Hybrid MPI-OpenMP programming is therefore applied to achieve good parallel computational performance. For more information about the process and efficiency of parallel computing on the proposed strategy, see Liu et al. (2019).

# 4 Validation and applications

In this section, the proposed strategy is validated using a spherical inclusion case. The effective elastic modulus of a spherical inclusion is computed both by the analytical homogenization
<figure><img src="image_4.png" /><figcaption>Figure 3: FMG scheme</figcaption></figure>

method and the computational homogenization method with the proposed strategy. A comparison of the effective modulus obtained by different methods is then realized. Two industrial materials applications are presented following the validation subsection.

## 4.1 Solution validation

To validate the proposed strategy, the typical spherical inclusion case is used. The cubic domain $\Omega$ is filled with a spherical inhomogeneity and a homogeneous matrix as presented in figure 4. The sphere radius is a quarter of the cube size $L$. The elastic modulus of the sphere is $C_s$, and $C_m$ denotes the elastic modulus of matrix. The ratio between these two materials is $r_e = E_m/E_s$, where $E_m$ and $E_s$ are the Young modulus of matrix and inclusion e.g. Figure 4 presents the case of $E_m = 500$ GPa and $E_s = 1$ GPa. The Poisson ratio equals to 0.3 for both materials. The objective is to compute the effective modulus $C$ of the domain $\Omega$ for different $r_e$. For the sake of simplicity, $E_s$ equals to 1 GPa.

<figure><img src="image_5.png" /><figcaption>Figure 4: Young modulus of inclusion and matrix</figcaption></figure>

Eshelby seems to have already solved the spherical inclusion. However, Eshelby analyzed it with an infinite matrix with is not for this case, and it is not possible to have an infinite domain. One of the best-known methods is the Mori-Tanaka (MT) homogenization method, see Mori and Tanaka (1973) and Benveniste (1987) for details. The effective modulus obtained by MT is:

$$C_{MT} = C_m + V_s(C_s - C_m)A$$

where

$$A = T[(V_m - V_s)I + V_sT]^{-1}$$
with

$$T = [I + S_s C_m^{-1} (C_s - C_m)]^{-1}$$

where *I* is the unit tensor, *S*<sub>s</sub> is the Eshelby tensor (see Gu et al. (2016) for details and Appendix 5 for its formulations), *V*<sub>s</sub> denotes the volume fraction of the sphere part, *C*<sub>MT</sub> is the effective elastic modulus obtained by MT homogenization.

For the computational homogenization with the proposed MG strategy, the homogeneous displacement boundary condition is applied to obtain the effective modulus *C*<sub>MG</sub>. Domain Ω is discretized by 256×256×256 elements. The under-relaxation coefficient equals to 0.5 for the Jacobi solver. The coarsest grid has 4×4×4 elements, thus there are 7 levels. The FMG cycles use *n*<sub>cy</sub> = 6, *ν*<sub>0</sub> = 50, *ν*<sub>1</sub> = 4 and *ν*<sub>2</sub> = 2. The comparison between *C*<sub>MT</sub> and *C*<sub>MG</sub> is qualified using the following equation:

$$e_{MT-MG} = \frac{1}{n_t} \sqrt{\sum_{i=1}^{n_t} \left( \frac{C_i^{MT} - C_i^{MG}}{C_i^{MT}} \right)^2}$$

where $C_i^{MT}$ is each term of $C_{MT}$. $C_i^{MG}$ is each term of $C_{MG}$. $n_t$ is the number of components in elastic modulus, which equals to 36.

The elastic modulus for $r_e = 10$ obtained by the MT and the MG is presented as the following:

$$C_{MT} = \begin{pmatrix} 11.933 & 4.997 & 4.997 & 0 & 0 & 0 \\ 4.997 & 11.933 & 4.997 & 0 & 0 & 0 \\ 4.997 & 4.997 & 11.933 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3.563 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3.563 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3.563 \end{pmatrix}$$

$$C_{MG} = \begin{pmatrix} 11.986 & 4.991 & 4.991 & 0 & 0 & 0 \\ 4.991 & 11.986 & 4.991 & 0 & 0 & 0 \\ 4.991 & 4.991 & 11.986 & 0 & 0 & 0 \\ 0 & 0 & 0 & 3.488 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3.488 & 0 \\ 0 & 0 & 0 & 0 & 0 & 3.488 \end{pmatrix}$$

From these two elastic modulus, we find that the effective material is isotropic, which confirms to this case. Table 1 presents the first item of the effective modulus obtained by the MT and the MG method for different $r_e$. $e_{MT-MG}$ is also presented in Table 1. From this table, a good agreement can be found between the MT and the MG. There is a small difference (less than 0.5%) between these two methods. The proposed method is therefore validated.

<table><caption>Table 1: $C_{1111}$ obtained by MT and MG</caption><thead><tr><th>r<sub>e</sub></th><th>C<sup>MT</sup><sub>1111</sub>/GPa</th><th>C<sup>MG</sup><sub>1111</sub>/GPa</th><th>e<sub>MT-MG</sub></th></tr></thead><tbody><tr><td>1</td><td>1.346</td><td>1.346</td><td>0</td></tr><tr><td>10</td><td>11.934</td><td>11.986</td><td>0.31%</td></tr><tr><td>50</td><td>58.220</td><td>58.557</td><td>0.40%</td></tr><tr><td>100</td><td>116.037</td><td>116.730</td><td>0.42%</td></tr><tr><td>200</td><td>231.6670</td><td>233.068</td><td>0.43%</td></tr><tr><td>500</td><td>578.549</td><td>582.073</td><td>0.43%</td></tr></tbody></table>

The convergence check is also carried out by computing the error $L_2$ norm on level $l-1$ using the converged solution on $l$ during FMG cycles. Figure 5 presents the error $L_2$ norm in
function of the element size for the spherical case $r_e = 10$. It demonstrates the linear relationship between the error $L_2$ norm and the element size $h$ in logarithmic scale. This is as expected for FE based multi-grid. The accuracy increases by a factor of 3 (4 for homogeneous materials) when the element size decreases by a factor of 2, which confirms the ability of the proposed multi-grid strategy to converge even for the large material property variation cases.

<figure><img src="image_6.png" /><figcaption>Figure 5: Convergence check</figcaption></figure>

The strain field of this spherical case with $r_e = 500$ is illustrated in Figure 6. It is the field of $\frac{Strain_{xx}}{Strain_o}$, where $Strain_{xx}$ is the xx strain, and $Strain_o$ is the strain at the macroscopic scale. The prescribed displacements to obtain this result are indicated in equation (18).

$$ \vec{u} = \{x, 0, 0\} \text{ on all the surfaces} \quad (18) $$

where $x$ is the coordinate of the X direction. Some strain raisers are found at the interface of matrix and inhomogeneity. The strain inside the sphere is smaller than in the matrix which is typical for the soft inclusions.

<figure><img src="image_7.png" /><figcaption>Figure 6: Strain field of spherical inclusion case.</figcaption></figure>
## 4.2 Efficiency analysis

The efficiency of the proposed strategy is also analyzed for this mechanical problem. The principle is to perform the numerical simulation for $r_e = 10$ of the spherical inclusion case with both the FMG scheme and the single level Jacobi iterative solver. The total cost of FMG cycles and the V-Cycle is described by the following equation (19) according to Venner and Lubrecht (2000).

$$ 
\begin{aligned}
W_{FMG} & \le \frac{n_{cy}}{1 - (H/h)^{-d}} W_{V-Cycle} \\
& \le \frac{n_{cy}(\nu_1 + \nu_2)}{(1 - (H/h)^{-d})^2} W_U
\end{aligned} 
\quad (19) $$

where $W_U$ is the cost of one relaxation on the finest grid. $W_{FMG}$ and $W_{V-Cycle}$ are the cost of FMG cycles and the V-Cycle. $H/h = 2$ in this work which presents the ratio between the grid size of coarse $l-1$ and fine grid $l$. $d$ is the problem dimension which is $d=3$ in this work. There is also the cost of transform between grids, but problem solved in this work are too large, the relaxation time is much more expensive than others routines. One accounts only the cost of relaxations on all grids.

Figure 4.2 implies that with 2000 relaxations, i.e. a cost of 2000 $W_U$, the single level Jacobi iterative solver still can not achieve the initial residual of the FMG scheme. The FMG achieves a residual of $10^{-5}$ in with 36 relaxations on the finest grid, i.e. a cost of $\frac{6(4+2)}{(1-2^{-3})^2} W_U \approx 47 W_U$. The FEM convergence slows down after only few relaxations. However, the FMG scheme remains its good convergence. The red line in Figure 7(b) demonstrates the residual at the finest level after each V-Cycle of the FMG scheme, where a good convergence can be found.

As mentioned above, the proposed method uses only the diagonal values of the stiffness matrix. For a $N$ nodes mechanical problem ( discretized using 8-node cubic elements), the memory requirement for the sparse global stiffness matrix is bounded by $(3 \times 27) \times (3 \times N) \times 8$ bytes (assuming 8 bytes floating point numbers) which gives 1944$N$ bytes. This estimate assumes that each node has 27 neighbors. However, there are also surface and edge nodes for which the number of neighbors is lower which makes the proposed estimate an upper bound of the actual memory cost. For the proposed method, the memory requirement for the diagonal of the stiffness matrix is $3N \times 8 = 24N$ bytes, which is 81 times cheaper than the global stiffness matrix. For a $l$ levels MG algorithm, the total size is:

$$ 24N \times \left[ 1 + \frac{1}{8} + \frac{1}{8^2} + \cdots + \frac{1}{8^l} \right] < 27.43N \text{ bytes.} $$

It demonstrates the memory requirement of the global sparse stiffness matrix is about 71 times more expensive than the proposed method. e.g. The largest problem in this work has more than 6 billion nodes, which means the size of the global sparse stiffness matrix is more than 10 TB. The proposed method needs only about 0.15 TB.

## 4.3 Laminated material simulations

In this subsection, a simulation on a problem with more than 18 billion DoF is presented. The numerical result is then compared to the experimental result obtained by the digital image correlation (DIC).

The CT image of a laminated material is used in this case. This CT image has 700×1300×1700 voxels with 4.5 µm/pixel. This laminated structure consists of four layers with a fiber orientation of 15°, -15°, -15° and 15°, respectively. The fiber is the E-glass fiber with a Young's modulus of 80.0 GPa. The matrix is s M9 epoxy with a Young's modulus of 3.2 GPa. The Poisson ratio of these two materials equals to 0.22. Figure 8 illustrates a section view of an interface of layers, one
<figure><img src="image_8.png" /><figcaption>Figure 7: Convergence of the Jacobi solver (a) and FMG scheme (b) on a 256³ nodes problem</figcaption></figure>

can observe two different fiber orientations crossing, and the fiber distributions are not uniform. All these defects, which can only be seen by CT, have an impact on the material mechanical behavior under certain loading. More information about this image and this material can be found in Lecomte-Grosbras et al. (2015).

The main subject of this subsection is to analyze free edge effects in the laminated structure by numerical simulation. A qualitative comparison between the numerical results and the experimental results is also presented. The free edge effect was firstly presented by Pipes and Pagano (1994), who found the strain concentrations around free edges and the ply interface. Lecomte-Grosbras et al. (2009) illustrates the free edge effects by the DIC experimental method of an unidirectional carbon fiber reinforced plastic laminated structure. For more information about free edge effects, see Pipes and Pagano (1994). In this paper, instead of carrying out DIC experiments, the CT image of the laminated material is used directly to employ numerical simulations and understand the free edge effect.
<figure><img src="image_9.png" /><figcaption>Figure 8: Sections of material</figcaption></figure>

To carry out numerical simulation, one proposes to take a part of this image which refers to the region of interest (ROI). The ROI is constituted of 577×1153×1153 voxels as presented in Figure 10. To have more voxels in the fiber section, one performs a sub-sampling, i.e. linear interpolation, on this ROI. Figure 11 shows that we have about ten pixels per fiber diameter. The final input image, i.e. domain Ω, has 1153×2305×2305 voxels which means we have more than 6 billion elements, i.e. 18 billion DoF. The coordinates of the center of Ω is (0,0,0) with a size of L×2L×2L.

The material property in each element is indispensable to perform numerical simulations. However, CT images contain the gray level (GL) which varies from 0 to 255. One has to choose a strategy to obtain material properties from CT images. The often used method is the thresholding method. As presented in Figure 12(a), 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). From the histogram of this laminated composite in Figure 12(b), the interface is not obvious neither. Instead of applying the thresholding method, one proposes to apply a continuous property.

$$ \alpha = \left(1 - e^{-\frac{|GL - 160.5|}{20}}\right) \cdot \mathrm{sign}(GL - 160.5) + 1 $$

where α is the material property variation, e.g. K,G. The parameters in this equation are obtained by the trial and error method to ensure a good agreement between the original image and material property distribution. This equation can be described by Figure 4.3.

<figure><img src="image_10.png" /><figcaption>Figure 9: Gray level and material property variations. "0" represents epoxy, "1" represents fibers</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: The gray level (GL) of the ROI opened in Paraview with colors</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Pixels in fiber section; red is fiber, blue is matrix</figcaption></figure>

The boundary conditions are given by the following equation:

$$ \begin{cases} \vec{u} = \{0, 0, -0.01L\}, & \text{on } z = -L \\ \vec{u} = \{0, 0, 0.01L\}, & \text{on } z = L \end{cases} \qquad (20) $$

where $u_z$ is the displacement in the Z direction.

There are 1153×2305×2305 elementary nodes on the finest level and 9×18×18 elements on the coarsest level with a total of 8 levels. The under-relaxation coefficient equals to 0.5 for the Jacobi solver. The parameters of the FMG cycles are: $n_{cy} = 9$, $\nu_0 = 500$, $\nu_1 = 4$ and $\nu_2 = 2$.

Figure 13 illustrates the Young modulus and the shear strain XZ on surface $y = -L$. This figure demonstrates that the shear strain field mimics the Young modulus distribution. Shear strain concentrations are found on the two interfaces. This is the so called free edge effect in laminated materials. In the area with few fibers with too much matrix, they also aches strain concentrations which can lead to preventive damage in industrial applications. The same phenomena are found by the DIC experimental method in Lecomte-Grosbras et al. (2009). The
<figure><img src="image_13.png" /><figcaption>Figure 12: Laminated composite image and its histogram</figcaption></figure>

material used in Lecomte-Grosbras et al. (2009) is a carbon fiber reinforced plastic laminated structure which is similar to the material used in this work. A good correlation between the numerical and the experimental results is found.

Figure 14 reveals the displacement in the Z direction on the two opposite surfaces $y = -L$ and $y = L$. Equally, the free edge effects can be found on the interfaces, which leads to the displacement variations on the interfaces. This displacement variation can also be found in Lecomte-Grosbras et al. (2009), the displacement curves illustrated in Figure 15, have the same tendency both for the numerical results and the DIC results. Another phenomenon that we can observe is that the displacement variations on $y = -L$ is almost anti-symmetric to $y = L$. This is because the fiber orientations on these two opposite surfaces are anti-symmetric. The Z displacements along the X axis in the center of these two surfaces are illustrated in Figure 15 where a clear anti-symmetry can be found.

<figure><img src="image_14.png" /><figcaption>Figure 13: Young modulus and strain on the surface $y = -L$; the white dashed lines are the interfaces of layers. Strain$_o$ is the macroscopic strain, which is 1%.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 14: Displacement $u_z$ on the surface $y = -L$ and on the surface $y = L$; the while dashed lines are two interfaces of layers, the green line is the $z = 0$ line which is the center of the surface</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Displacement $u_z$ along the line $z = 0$ on surface $y = -L$ and $y = L$ (see Figure 14 for position of line)</figcaption></figure>

## 4.4 Cast iron applications

In this work, we perform the proposed strategy on a nodular graphite cast iron CT image. The image used in this paper is the image obtained by Rannou et al. (2010) using X-Ray tomography. This image consists of 340×340×512 voxels with a 5.06 µm/pixel size. For more information about this CT image, please refer to Rannou et al. (2010). A ROI with 257×257×257 voxels is taken from this CT image. As illustrated in Figure 16, many carbon nodules with a random distribution, can be found at the microscopic scale of the CT image.

The objective is to perform the linear elastic simulation on the ROI, with a prescribed rectangular crack to see how carbon nodules affect the crack opening. Assuming $\Omega$ is the ROI. The center of the $\Omega$ is the origin of axis. The size of $\Omega$ is $L$. 256³ 8-node cubic elements are used to discretize $\Omega$. The rectangular crack is presented in Figure 17. The width of the crack is the size of the cube $L$. Its length is $L/3$. The crack is constructed by setting material property as 0 one three layers of nodes in Z direction. The prescribed boundary conditions are given by
<figure><img src="image_17.png" /><figcaption>Figure 16: The GL of the ROI of the CT image of cast iron</figcaption></figure>

the following equations.

$$
\begin{cases}
u_z = 0, & \text{on } Z = -\frac{L}{2} \\
u_z = 0.01L, & \text{on } Z = \frac{L}{2} \\
\vec{u} = \vec{0}, & \text{at node } (0,0, -\frac{L}{2})
\end{cases}
$$

<figure><img src="image_18.png" /><figcaption>Figure 17: The prescribed rectangular crack in cast iron CT image</figcaption></figure>

Here, the thresholding method, i.e. Figure 4.4, is used to obtain the material property due to an obvious boundary between iron and carbon nodular. To confirm this choice, the volume fraction of carbon nodular after the thresholding is equal to the value provided by the manufacturing company.

The material properties are given in Table 2. Table 2 implies that the crack is defined by setting the material property to 0 at the crack nodes. The crack domain is considered as nodes with a material property 0, which means the computational domain is the entire cubic Ω.
<figure><img src="image_19.png" /><figcaption>Figure 18: Thresholding method to obtain material property from GL in image</figcaption></figure>

To perform this numerical simulation, a specific treatment is required to deal with the infinite material property jump. This treatment involves in the following steps.

* Step A: The MF-FEM iterative solver for the relaxation step.

* Step B: The injection of solutions from the fine grid to the coarse grid.

* Step C: The restriction of residuals from the fine grid to the coarse grid.

* Step D: The prolongation of the corrections from the coarse grid to the fine grid.

<table><caption>Table 2: Material properties in cast iron</caption><thead><tr><th>Component</th><th>Young's modulus/GPa</th><th>Poisson ratio</th></tr></thead><tbody><tr><td>Iron</td><td>210</td><td>0.3</td></tr><tr><td>Carbon granules</td><td>21</td><td>0.2</td></tr><tr><td>Crack</td><td>0</td><td>0</td></tr></tbody></table>

For "Step A": during the relaxation processing, the displacements on nodes without material are not updated.

For "Step B": as mentioned above, the Reuss approximation is

$$K^{HR} = \frac{N_h}{\sum_{1}^{N_h} \frac{1}{K^h}}, \quad G^{HR} = \frac{N_h}{\sum_{1}^{N_h} \frac{1}{G^h}}$$

To avoid 0 in the denominator, one defines: if there is $K^h = 0$ and $G^h = 0$, then $K^{HR} = 0$ and $G^{HR} = 0$.

For "Step C" and "D": the restriction and the prolongation process are done only when the material property on this node and on all of its nearest neighbor nodes is not 0 at the fine grid. If 0 appears in the denominator when computing the restriction matrix, one replaces it with $10^{-6}$.

For the grids, one has 256×256×256 elements on the finest level and 4×4×4 elements on the coarsest level with 7 levels. The under-relaxation coefficient is taken to 0.5 for the Jacobi solver. The parameters of the FMG cycles are: $n_{cy} = 4$, $\nu_0 = 100$, $\nu_1 = 8$ and $\nu_2 = 4$. Since the displacement and the crack thickness are too small for the visualization, i.e. 1%, in the following figures, the amplitude of mesh deformation is multiplied by 20..
Figure 19 illustrates the Young's modulus and the strain field in the cast iron. It presents the strain concentrations on the crack front. The largest strain can be found in carbon nodules on the crack front. The largest strain is 10 times larger than the prescribed strain at the macroscopic scale, i.e. 1%. Other strain concentrations are located on carbon nodules, because carbon nodules are 10 times softer than iron.

<figure><img src="image_20.png" /><figcaption>Figure 19: The Young's modulus and the strain<sub>zz</sub> in cast iron. Strain<sub>o</sub> = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.</figcaption></figure>

Another simulation is carried out to compare the crack opening in a homogeneous material
and in a heterogeneous material. The principle is to replace all the carbon nodules in the CT
image by iron, which means the simulation is carried with a prescribed crack in the iron. The
same boundary conditions are applied on the simulation of the crack opening in iron. Figure
20 shows the strain field in the homogeneous iron and heterogeneous cast iron. The "butterfly"
strain field in the iron presents typical strain concentrations in the homogeneous materials after
the crack opening. Compared to the strain concentrations in the iron, the strain concentrations
in the cast iron are not only in the vicinity of the crack front but also in the carbon nodules
over the entire volume. The material heterogeneity spreads strain concentrations in a large
volume. The "butterfly" strain field is destroyed by carbon nodules. Figure 21 illustrates the
strain concentrations on the crack front in the iron and in the cast iron, respectively. The largest
strain is located on the carbon nodules on the crack front in the cast iron, but in iron, a uniform
strain concentration can be found on the entire crack front. This should lead to non rectilinear
crack front (if crack propagation is performed) and strong interactions between the crack and
the material microstructure.

**5 Conclusions**

In this paper, several strategies are proposed to achieve a CT simulation with more than 16
billion of degree of freedom. The application of the matrix free finite element method permits
us to carry out such large simulation with the minimum memory space. The MultiGrid method
improves convergence compared to a single level Jacobi iterative solver. The use of the proposed
<figure><img src="image_21.png" /><figcaption>Figure 20: $Strain_{zz}$ in iron and in cast iron. $Strain_o$ = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 21: Strain concentrations on the crack front in the cast iron and in the iron. $Strain_o$ = 1%, which is the macroscopic strain. The amplitude of mesh deformation is multiplied by 20.</figcaption></figure>

homogenization techniques improve the convergence stability of the MultiGrid method. The
applications presented in this paper confirms the CT image based simulations permit one to analyze the materials behavior at the microscopic scale. The good correlation between the DIC experiments and the CT simulations for the laminated structures confirms the strategy proposed in this paper. The free edge effect in the laminated materials proves the importance to carry out numerical simulations at the microscopic scale. The analysis of the crack opening application in the cast iron demonstrates how the carbon nodules affect the strain concentrations in the materials. It implies the importance to account for the carbon nodules in the numerical simulation and also in the application of materials in the industrial domain. Since the largest strain concentrations are located in the soft carbon nodules, the damage starts firstly in this areas. Both these applications imply the importance to carry out simulations at the microscopic scale by using the CT images. Performing numerical simulations with real material microstructures permits ones to observe many phenomena that one can never observe at the macroscopic scale or the theoretical microstructure.

## Acknowledgments

We gratefully thank for the support of Région Pays de la Loire and Nantes Métropole through a grant Connect Talent IDS.

## Appendix

The following Eshelby's tensor is for the spherical inclusion of a linear elastic problem in macroscopically isotropic materials.

$$ S_s = \frac{1}{15(1 - \nu)} \begin{pmatrix} 7 - 5\nu & 5\nu - 1 & 5\nu - 1 & 0 & 0 & 0 \\ 5\nu - 1 & 7 - 5\nu & 5\nu - 1 & 0 & 0 & 0 \\ 5\nu - 1 & 5\nu - 1 & 7 - 5\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 - 5\nu & 0 & 0 \\ 0 & 0 & 0 & 0 & 4 - 5\nu & 0 \\ 0 & 0 & 0 & 0 & 0 & 4 - 5\nu \end{pmatrix} $$

where $\nu$ is the Poisson's ratio of the inclusion.

## References

Alcouffe, R., Brandt, A., Dendy, Jr, J., Painter, J., 1981. The multi-grid method for the diffusion equation with strongly discontinuous coefficients. SIAM Journal on Scientific and Statistical Computing 2 (4), 430-454.

Arbenz, P., van Lenthe, G. H., Mennel, U., Müller, R., Sala, M., 2008. A scalable multi-level preconditioner for matrix-free $\mu$-finite element analysis of human bone structures. International Journal for Numerical Methods in Engineering 73 (7), 927-947.

Augarde, C., Ramage, A., Staudacher, J., 2006. An element-based displacement preconditioner for linear elasticity problems. Computers & structures 84 (31-32), 2306-2315.

Benveniste, Y., 1987. A new approach to the application of mori-tanaka's theory in composite materials. Mechanics of materials 6 (2), 147-157.

Bessho, M., Ohnishi, I., Matsuyama, J., Matsumoto, T., Imai, K., Nakamura, K., 2007. Prediction of strength and strain of the proximal femur by a ct-based finite element method. Journal of biomechanics 40 (8), 1745-1753.
Boffy, H., Baietto, M.-C., Sainsot, P., Lubrecht, A. A., 2012. An efficient 3d model of heterogeneous materials for elastic contact applications using multigrid methods. Journal of tribology 134 (2), 021401.

Boffy, H., Venner, C. H., 2014. Multigrid solution of the 3d stress field in strongly heterogeneous materials. Tribology international 74, 121–129.

Brandt, A., 1973. Multi-level adaptive technique (mlat) for fast numerical solution to boundary value problems. In: Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics. Springer, pp. 82–89.

Brandt, A., 1977. Multi-level adaptive solutions to boundary-value problems. Mathematics of computation 31 (138), 333–390.

Carey, G. F., Jiang, B.-N., 1986. Element-by-element linear and nonlinear solution schemes. Communications in Applied Numerical Methods 2 (2), 145–153.

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

Gu, H., Réthoré, J., Baietto, M.-C., Sainsot, P., Lecomte-Grosbras, P., Venner, C. H., Lubrecht, A. A., 2016. An efficient multigrid solver for the 3d simulation of composite materials. Computational materials science 112, 230–237.

Hashin, Z., 1979. Analysis of properties of fiber composites with anisotropic constituents. Journal of Applied Mechanics 46 (3), 543–550.

Hoekema, R., Venner, K., Struijk, J. J., Holsheimer, J., 1998. Multigrid solution of the potential field in modeling electrical nerve stimulation. Computers and Biomedical Research 31 (5), 348–362.

Hughes, T. J., Ferencz, R. M., Hallquist, J. O., 1987. Large-scale vectorized implicit calculations in solid mechanics on a cray x-mp/48 utilizing ebe preconditioned conjugate gradients. Computer Methods in Applied Mechanics and Engineering 61 (2), 215–248.

Kronbichler, M., Ljungkvist, K., 2019. Multigrid for matrix-free high-order finite element computations on graphics processors. ACM Transactions on Parallel Computing (TOPC) 6 (1), 1–32.

Lecomte-Grosbras, P., Paluch, B., Brieu, M., De Saxcé, G., Sabatier, L., 2009. Interlaminar shear strain measurement on angle-ply laminate free edge using digital image correlation. Composites Part A: Applied Science and Manufacturing 40 (12), 1911–1920.

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

Lengsfeld, M., Schmitt, J., Alter, P., Kaminsky, J., Leppek, R., 1998. Comparison of geometry-based and ct voxel-based finite element modelling and experimental validation. Medical engineering & physics 20 (7), 515–522.

Liu, X., Réthoré, J., Baietto, M.-C., Sainsot, P., Lubrecht, A. A., 2019. An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity. Computational materials science 166, 265–275.
Michailidis, N., Stergioudi, F., Omar, H., Tsipas, D., 2010. An image-based reconstruction of the 3d geometry of an al open-cell foam and fem modeling of the material response. Mechanics of Materials 42 (2), 142–147.

Mori, T., Tanaka, K., 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta metallurgica 21 (5), 571–574.

Moulinec, H., Suquet, P., 1995. A fft-based numerical method for computing the mechanical properties of composites from images of their microstructures. In: IUTAM symposium on microstructure-property interactions in composite materials. Springer, pp. 235–246.

Nguyen, T.-T., Rethore, J., Yvonnet, J., Baietto, M.-C., 2017. Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials. Computational Mechanics 60 (2), 289–314.

Özdemir, I., Brekelmans, W., Geers, M., 2008. Computational homogenization for heat conduction in heterogeneous solids. International journal for numerical methods in engineering 73 (2), 185–204.

Pipes, R. B., Pagano, N., 1994. Interlaminar stresses in composite laminates under uniform axial extension. In: Mechanics of composite materials. Springer, pp. 234–245.

Proudhon, H., Li, J., Wang, F., Roos, A., Chiaruttini, V., Forest, S., 2016. 3d simulation of short fatigue crack propagation by finite element crystal plasticity and remeshing. International Journal of Fatigue 82, 238–246.

Rannou, J., Limodin, N., Réthoré, J., Gravouil, A., Ludwig, W., Baïetto-Dubourg, M.-C., Buffiere, J.-Y., Combescure, A., Hild, F., Roux, S., 2010. Three dimensional experimental and numerical multiscale analysis of a fatigue crack. Computer methods in applied mechanics and engineering 199 (21-22), 1307–1325.

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

van Rietbergen, B., Weinans, H., Huiskes, R., Polman, B., 1996. Computational strategies for iterative solutions of large fem applications employing voxel data. International Journal for Numerical Methods in Engineering 39 (16), 2743–2767.

Venner, C. H., Lubrecht, A. A., 2000. Multi-level methods in lubrication. Elsevier.

Watremetz, B., Baietto-Dubourg, M., Lubrecht, A., 2007. 2d thermo-mechanical contact simulations in a functionally graded material: A multigrid-based approach. Tribology international 40 (5), 754–762.