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

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.


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. One of the most used methods is to B Julien Réthoré julien.rethore@ec-nantes.fr 1 GeM, CNRS UMR 6183, Centrale Nantes, 44321 Nantes, France 2 INSA-Lyon, CNRS UMR5259, LaMCoS, Univ Lyon, 69621 Lyon, France 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 [13] explains how to predict the macroscopic property of a fiber composite from its microscopic structure. Özdemir et al. [25] 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. [19] and Bessho et al. [5] 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. [11], Michailidis et al. [21] and Proudhon et al. [27] analyze properties of industrial materials by applying FEM simulations on CT images. Another wellknown method is the fast fourier transform (FFT). One of the most typical examples is presented in [23], which introduces a FFT-based numerical method to simulate the mechani-cal properties of composite materials. The finite difference method (FDM) is also often used to carry out CT simulations. Gu et al. [12] 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 [24] 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. [12] 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 [10] in 1980s. It was investigated by Augarde et al. [3], Hughes et al. [15] and van Rietbergen et al. [30] during the last 3 decades. Nowadays, it is often used in parallel computing, see Tezduyar et al. [29].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. [2] 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 [15] 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 Sect. 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 [2]. In this paper, the MFFEM is 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 [16] 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 [20] for the efficiency of MG algorithms applied on the FEM. MG methods were first introduced in the 1960s, and then, enriched by Brandt [8,9] in the 1970s. It is often used with the FDM. Work of Venner and Lubrecht [31] developed an efficient and robust FDM based MG solver for the elasto-hydrodynamic lubrication problem. The work presented in [32] introduced the application of MG methods for a 2D graded problem. Boffy et al. [6,7] employed MG FDM for small heterogeneous material mechanical problems. The work of Gu et al. [12] 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: Sect. 2 reviews the fundamentals of linear elasticity and its discretization by the proposed MF-FEM; Sect. 3 introduces the strategy in detail, both the MG method with special intergrid operators and parallel computing are applied to improve efficiency; Sect. 4 presents the validation of proposed method and some applications; Sect. 5 provides discussions and conclusions.

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.

Governing equations
Assuming a deforming domain Ω as presented in Fig. 1, its boundary is denoted by ∂Ω. Its exterior normal is n. The second order stress tensor is described by σ . The body force in Ω, which is due to gravity, magnetism, etc., is f. The acceleration is presented by a. The density is denoted by ρ. The equations of equilibrium can be described as: For the sake of simplicity, inertia effects are neglected in this paper. The equation of equilibrium can be simplified as: Besides the equilibrium equations, the boundary conditions on ∂Ω are essential to solve a mechanical problem. Two often used boundary conditions are Dirichlet and Neumann boundary conditions. The Dirichlet boundary condition applied on ∂Ω D is often presented as prescribed displacements, e.g. u = u 0 on 1 . The Neumann boundary condition applied on ∂Ω N refers to the external force f ex , which can be described as: To solve Eq. (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 elements are therefore used to discretize Ω. With this strategy, the mesh generation step becomes automatic.
Multiplying Eq. (2) by a test function u * and integrating over Ω, the weak form of Eq. (2) reads: Applying the divergence theorem and integrating by parts, it reads: which can be described as the equilibrium of internal force f in and external force f ex : Applying finite element discretization, the displacements u can be described by: 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û. The component c of the displacement at the node i is described byû c i . The shape function ψ is the typical shape function of a 8-node cubic element. The test function u * reads: where j is the node number.
The internal force for node j and component c reads: 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 σ and u is implied by σ (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 Eq. (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 [20] for the thermal conduction problem.
In this work, it can be described by Eq. (9).
whereû ite+1 andû ite are the displacements of one node in the current iteration and previous iteration. The relaxation coefficient is denoted by ω, for 0 < ω < 1, it refers to damped Jacobi, i.e. under-relaxation, which is used in this work. The diagonal matrix S is the compliance of the diagonal value of the stiffness matrix at this node. It reads: where K u , K v and K w are the diagonal values of the stiffness matrix at each DoF. To obtain 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.

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: where ε 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 I. Using with the Young's modulus E and the Poisson ratio υ, K and G read: Strain ε can be described by displacement u: Its trace being defined as Combining Eqs. (11) and (12), (11) reads: Applying the finite element discretization with Eq. (7), the stress reads: where k = {1, 2, 3}. δ cm is the Kronecker delta. It reads: Combining Eqs. (8) and (14), one obtains the component c of the internal force on node j: 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: where K u is the first diagonal term of S in Eq. (10). It reads: 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.

MultiGrid
As presented in [20], 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 3 -grid problem with 3 levels, the number of points on level is: 128 3 on level 3, 64 3 on level 2 and 32 3 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 [20], 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 [20].
The process to obtain the 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. [1] for a 2D problem. Hoekema et al. [14] extended it to 3D. This work proposes to apply the homogenization method proposed in [20], 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: where K H V , G H V , K H R and G H R 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: 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 [20], 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 Fig. 2. The subject is to bring the corrections, i.e. displacement correction 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: For red points, one computes them from the e at black point of level l, e.g. point B1: For blue points, one computes them from the e at red points on level l, e.g. point C1: For yellow points, one computes them from the e at blue points on level l, e.g. point O: 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 bestknown 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 Fig. 3. On level 1, one performs ν 0 relaxations, ν 1 relaxations are performed on each level going up. ν 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 [20].

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 method and the computational homogenization method with the proposed strategy. A comparison of the effective modulus obtained by different methods is then

Solution validation
To validate the proposed strategy, the typical spherical inclusion case is used. The cubic domain Ω is filled with a spherical inhomogeneity and a homogeneous matrix as presented in Fig. 4. The sphere radius is a quarter of the cube size L. The elastic tensor of the sphere is C s , and C m denotes the elastic tensor of matrix. The stiffness ratio etween 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. Fig. 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 tensor C of the domain Ω for different r e . For the sake of simplicity, E s equals to 1 GPa.
Eshelby seems to have already solve the spherical inclusion. However, Eshelby analyzed it with an infinite matrix what is not exact for this case because it is not possible to have an infinite domain. One of the best-know methods is the Mori-Tanaka (MT) homogenization method, see [4,22] for details. The effective modulus obtained by MT is: where I is the unit tensor, S s is the Eshelby tensor (see [12] for details and "Appendix" section for its formulations), V s denotes the volume fraction of the sphere part, C MT is the effective elastic tensor obtained by MT homogenization. For the computational homogenization with the proposed MG strategy, the homogeneous displacement boundary condition is applied to obtain the elastic tensor C M G . Domain Ω is discretized by 256 × 256 × 256 elements. The underrelaxation 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 cy = 6, ν 0 = 50, ν 1 = 4 and ν 2 = 2.
The comparison between C MT and C M G is qualified using the following equation: is each term of C M G . n t is the number of components in elastic tensor, which equals to 36.
The elastic tensor for r e = 10 obtained by the MT and the MG is presented as the following: 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. 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. The strain field of this spherical case with r e = 500 is illustrated in Fig. 6. It is the field of Strain xx Strain o , where Strain x x is the xx strain, and Strain o is the strain at the macroscopic scale. The prescribed displacements to obtain this result are indicated in Eq. (18). u = {x, 0, 0} on all the surfaces (18) where x is the coordinate in 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.

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 Eq. (19) according to Venner and Lubrecht [31].
where W U is the cost of one relaxation on the finest grid. W F MG 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 7 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 6(4+2) The FEM convergence slows down after only few relaxations. However, the FMG scheme remains its good convergence. The red line in Fig. 7b 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 × 27) × (3 × N ) × 8 bytes (assuming 8 bytes floatting point numbers) which gives 1944N 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 × 8 = 24N bytes, which is 81 times cheaper than the global stiffness matrix. For a l levels MG algorithm, the total size is: 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.

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 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 [18].
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 [26], who found the strain concentrations around free edges and the ply interface. Lecomte-Grosbras et al. [17] 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 [26]. 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 (Fig. 12).
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 Fig. 9. To have more voxels in the fiber section, one performs a sub-sampling, i.e. linear interpolation, on this ROI. Figure 10 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 (G L) 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 Fig. 11a, 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 Fig. 11b, the interface is not obvious neither. Instead of applying the thresholding method, one proposes to apply a continuous property.
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 Fig. 12.
The boundary conditions are given by the following equation: 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, ν 0 = 500, ν 1 = 4 and ν 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. Fig. 12 Gray level and material property variations. "0" represents epoxy, "1" represents fibers 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 [17]. The material used in [17] 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 [17], the displacement curves illustrated in Fig. 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 Fig. 15 where a clear anti-symmetry can be found.

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. [28] 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 [28]. A ROI with 257 × 257 × 257 voxels is taken from this CT image. As illustrated in Fig. 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 and strain concen- Assuming Ω is the ROI. The center of the Ω is the origin of axis. The size of Ω is L. 256 3 8-node cubic elements are used to discretize Ω. The rectangular crack is presented in Fig. 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 the following equations. Here, the thresholding method, i.e. Fig. 18, 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 Ω. 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.  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 To avoid 0 in the denominator, one defines: if there is K h = 0 and G h = 0, then K H R = 0 and G H R = 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, ν 0 = 100, ν 1 = 8 and ν 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.
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 per-formed) and strong interactions between the crack and the material microstructure.

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 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 using the theoretical microstructure.