# An efficient MultiGrid solver for the 3D simulation of composite materials

Hanfeng Gu, Julien Réthoré, Marie-Christine Baietto, Philippe Sainsot, Pauline Lecomte-Grosbras, C.H. Venner, Antonius Lubrecht

► To cite this version:

Hanfeng Gu, Julien Réthoré, Marie-Christine Baietto, Philippe Sainsot, Pauline Lecomte-Grosbras, et al. An efficient MultiGrid solver for the 3D simulation of composite materials. Computational Materials Science, 2016, 112, Part A, pp.230-237. 10.1016/j.commatsci.2015.10.025. hal-01234623

HAL Id: hal-01234623

https://hal.science/hal-01234623

Submitted on 3 May 2021

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 MultiGrid solver for the 3D simulation of composite materials

Hanfeng Gu<sup>a,*</sup>, Julien Réthore<sup>a</sup>, Marie-Christine Baietto<sup>a</sup>, Philippe Sainsot<sup>a</sup>, Pauline Lecomte-Grosbras<sup>b</sup>, Cornelis H. Venner<sup>c</sup>, Antonius A. Lubrecht<sup>a</sup>

<sup>a</sup> LaMCoS, CNRS UMR 5259, INSA de Lyon, 18-20 rue des Sciences, 69621 Villeurbanne, France

<sup>b</sup> LML, CNRS UMR 8107, Ecole Centrale de Lille, 59650 Villeneuve-d'Ascq, France

<sup>c</sup> Faculty of Mechanical Engineering, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands

The subject of this paper is to present the application of MultiGrid (MG) methods in 3D composite mate-rial simulations through the solution of the elastic equations. An efficient MG solver is further developed for modeling composite structures with strong discontinuities. Different types of boundary conditions are imposed in the solver. The model is validated by comparing the MG numerical results with the theoretical results existing in the literature and Finite Element (FE) results and a good agreement is found. The potential of MG methods with respect to the homogenization process is illustrated. Then an ideal lami-nated structure is analyzed and a real topology is simulated for the first time by a MG model. The effect of fiber orientation, interface layer thickness, fiber layer thickness and the ratio of material properties on the surface displacement are investigated. MG results show the detailed local behavior and provide new insights into possible initiation of delamination.

## 1. Introduction

Composite materials have been widely used in many industrial fields in past decades due to their excellent mechanical properties as well as the adaptation under specific loading conditions. Many researchers have investigated these mechanisms and improved the performance. However, the random nature of its micro-structure, which may consist of voids, cracks or broken fibers, etc. makes the simulation of composite structures expensive or impractical. One of the currently used methods is to consider the heterogeneous material as a “homogeneous” layer and obtain the averaged material properties of each component. This kind of approximation is commonly referred to as the homogenization process [1], which takes the microscopic structure into account and is capable to predict the macroscopic mechanical behavior of a composite.

However, despite the globally correct predictions, the homoge-nization results lose the microscopic details related to the local physical mechanisms. These mechanisms may induce premature material failure, such as cracks. Meantime, this process requires fine grids to simulate the complex structure, which usually results in a heavy computational cost.

MultiGrid (MG) methods are numerical algorithms for solving differential equations using a hierarchy of discretizations. The advantage is that the cost of MG methods depends linearly on the number of unknowns, which provides a unique potential for solving large scale problems. Although they were first introduced in the early 1960s, it was not until the mid-seventies that Brandt [2,3] enriched the related theories and made them efficient in var-ious scientific and engineering fields. Lately a number of research-ers combined MG methods in their work. Venner and Lubrecht [4] used the MultiLevel techniques to develop a very efficient and stable solver for the ElastoHydrodynamic Lubrication (EHL) prob-lems. Watremetz et al. [5] applied the MG solutions to 2D graded materials and Boffy et al. [6,7] applied the MG solutions to 3D heterogeneous materials problems.

The topic of this paper is to present an efficient and stable solver for modeling composite structures. First, related MultiGrid tech-niques and simulation details are explained. Second, the model is validated from microscopic and macroscopic aspects in a simple case where a homogeneous matrix contains a spherical inhomo-geneity. Microscopic solutions are validated by comparing the local stress fields with Eshelby's analytic solution [8,9]. Macroscopic solutions are validated by comparing the homogenized elastic constants with Mori-Tanaka's prediction [10,11]. Finally, we apply the model to the laminated composite structure case and compare the results with finite element results. Then we apply the model to

<footnote>* Corresponding author.
E-mail address: hanfeng.gu@insa-lyon.fr (H. Gu).</footnote>
## Nomenclature

* $\epsilon, \sigma$ additional disturbed strain and stress due to the presence of inhomogeneity

* $\epsilon^*$ eigenstrain

* $\epsilon^\infty, \sigma^\infty$ remote uniform strain applied and uniform stress induced

* $\Delta w$ global averaged displacement

* $\lambda, \mu$ Lamé coefficients

* **C, C*** stiffness tensor of matrix and inhomogeneity

* **I** unit tensor

* **S** Ehselby tensor

* $\mu_{\text{equivalent}}^{\text{MT}}$ equivalent shear modulus

* $\nabla$ gradient operator

* $v$ Poisson's ratio

* **E** Young's modulus

* $E_m, E_f$ Young's modulus of matrix and fiber

* $F_i$ body force

* $I_h^H, I_H^H$ intergrid operators

* **L** cube length

* $L^h, L^H$ difference operator on fine and coarse grid

* **r** radius of the sphere or fiber

* **u, v, w** displacements

* $V_f$ volume percent of the inhomogeneity in matrix

the real structure to find a similar case which has perfect cylindrical fibers. Comparable results are found between the real case and the ideal case. The possible factors which may affect the delamination process are investigated as well.

## 2. Theory

The constitutive equations to be solved are the 3D linear elastic equations expressed by the displacements $u, v$ and $w$ in Cartesian coordinates:

$$ (\lambda + \mu)u_{jji} + \mu\nabla^2 u_i + F_i = 0, \quad i, j = 1, 2, 3 \quad (1) $$

where $\nabla$ is the gradient operator:

$$ \nabla = \left( \frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z} \right) \quad (2) $$

and $\lambda, \mu$ are called the Lamé coefficients which are functions of the Young's modulus $E$ and the Poisson's ratio $\nu$ in 3D space $\Omega$:

$$ \lambda(x, y, z) = \frac{E(x, y, z)v(x, y, z)}{(1 + v(x, y, z))(1 - 2v(x, y, z))} \quad (3) $$

$$ \mu(x, y, z) = \frac{E(x, y, z)}{2(1 + v(x, y, z))} \quad (4) $$

and $F_i$ is the body force. In this study, we assume it to be zero and $\Omega$ is the 3D rectangular domain. Eq. (1) describes the static equilibrium in $\Omega$, while on the boundary $\partial\Omega$ displacement boundary conditions and stress boundary conditions can both be considered.

### 2.1. MultiGrid theory

Iterative solution methods like Gauss-Seidel or Jacobi iteration usually have a good convergence speed for small problems while becoming very slow for large problems. This slowness can be explained by the presence of low frequency error components. The main idea of the MG method is to accelerate the convergence speed by solving the target equations on different size grids, which

means that different grids solve different frequency errors. A smooth error on fine grids can be mighty oscillatory on coarser grids.

The first step of MG methods is to discretize the target Eq. (1) on the target domain $\Omega^h$. Here, a second order accurate difference scheme is selected (for details see Appendix A). This step follows the pioneering work of Boffy [12]. The target equations can be simply denoted as $Lu = f$.

The simplest MG algorithm contains two different size grids: one coarse grid with size $H$ and one fine grid with size $h$ where $H > h$. For the nonlinear problem solved in this work, the Full Approximation Scheme (FAS) is selected for the coarse grid correction cycle. A two level MG scheme is given by (see Fig. 1):

(1) Relax $v_1$ times on $\Omega^h$ to obtain an approximation $\bar{u}^h$ to $u^h$.

(2) Calculate residuals $r^h = f^h - L^h \bar{u}^h$, restrict $r^h$ and $\bar{u}^h$ from grid $\Omega^h$ to grid $\Omega^H$ respectively, get $r^H = I_h^H r^h$ and $\bar{u}^H = I_h^H \bar{u}^h$.

(3) Relax $v_0$ times on $\Omega^H$: $L^H \bar{u}^H = r^H + L^H \bar{u}^H$, obtain an approximation $\bar{u}^H$.

(4) Interpolate the difference $\bar{u}^H - \bar{u}^h$ from coarse grid $\Omega^H$ to fine grid $\Omega^h$, obtain the corrected approximation: $\bar{u}^h = \bar{u}^h + I_H^h (\bar{u}^H - \bar{u}^h)$.

(5) Relax $v_2$ times on $\Omega^h$.

This scheme can be used recursively for any number of levels.

#### 2.1.1. Local grid refinement

In a composite structure simulation, local grid refinements can be very useful when a locally much finer grid is required. It can reduce the computing time and allocated memory dramatically compared with a global fine grid. The extension from a basic MG

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Two level MultiGrid scheme.</figcaption></figure>
scheme could be direct: the fine grid residual $r^h$ is only calculated where fine grids exist and excludes points on the “boundary” between fine and coarse grid. Likewise, the relaxation step is only carried out where the fine grid exists.

## 2.2. Operators

The strong material discontinuity existing in a composite structure, such as the interface between the fibers and the matrix, makes the standard MG operators inefficient. For example, Boffy et al. [6] observed divergence using standard MG operators in modeling the heterogeneous material for elastic contact applications.

Alcouffe et al. [13] have found that the MG relaxation should smooth the error in the product of the strongly discontinues coefficients and the derivative of the solution rather than the error in the derivative only, and proposed several alternative methods to solve the diffusion equation with strongly discontinuous coefficients in 2D. Based on this work, Hoekema et al. [14] developed the method into a 3D potential field for nerve simulations and Boffy et al. [7] developed the method into the 3D heterogeneous material field. For detailed information about the construction of intergrid operators ($I_h^H$ and $I_h^h$) see [7,13,14].

Another operator to improve is the coarse grid operator: $L^H$. Since the simple coarse grid operator cannot reflect the information correctly from the fine grid operator ($L^h$) in strong discontinuous cases, a more precise operator is required using Galerkin coarse grid approximation:

$$L^H = I_h^H L^h I_h^h \qquad (5)$$

Eq. (5) results in a 27-point operator instead of a 7-point operator for the simple coarse grid operator. The most important disadvantage of the 27-point operator is the required larger allocated memory compared to the simpler case.

## 2.3. Boundary conditions

According to the problem studied, the boundary conditions can be divided into 3 types: displacements specified on $\partial\Omega$, stress specified on $\partial\Omega$ or displacements specified on a portion $\partial\Omega_1$ and stress specified on a portion $\partial\Omega_2$. For the homogenization process in a composite material, there are two interesting boundary conditions extended from the above boundary conditions: one called Homogeneous Strain Boundary Condition (HSBC) and another called Periodic Boundary Condition (PBC).

HSBC [15] is suitable for the case where the representative volume element (RVE) of the composite structure is not periodic. The boundary condition on $\partial\Omega$ is:

$$u_i = \bar{\epsilon}_{ij}x_j, \quad i, j = 1, 2, 3 \qquad (6)$$

If the applied strain $\bar{\epsilon}_{ij}$ is given, then all the displacements on $\partial\Omega$ can be specified according to Eq. (6).

For a periodic RVE, PBC is more suitable [16,17]. The boundary condition on $\partial\Omega$ is:

$$u_i = \bar{\epsilon}_{ij}x_j + u_i^*, \quad i, j = 1, 2, 3 \qquad (7)$$

where $u_i^*$ is the unknown periodic displacement on the boundary $\partial\Omega$ and is dependent on the given global strain field $\bar{\epsilon}_{ij}$. Hence for a given strain $\bar{\epsilon}_{ij}$, the displacement $u_i$ on $\partial\Omega$ can also be determined by Eq. (7).

<figure><img src="image_3.png" /><figcaption>Fig. 2. Cubic matrix containing a spherical inhomogeneity.</figcaption></figure>

## 3. Validation

In this section, the validation of the MultiGrid results are implemented in both microscopic and macroscopic aspects. This simulation considers a common case where a spherical inhomogeneity is embedded in a homogeneous matrix (see Fig. 2). Let $\Omega$ be the matrix domain with the elastic stiffness tensor $C$ and $D$ be the spherical domain with the elastic stiffness tensor $C^*$. The accuracy of the MG method can be tested by comparing the perturbed local stress field and global stress-strain response induced by the presence of the inhomogeneity with existing theories.

### 3.1. Microscopic: Eshelby's solution

The local stress field inside the elliptical inhomogeneity embedded in the homogeneous matrix which is subjected to a remote uniform strain field can be obtained from Eshelby's "equivalent inclusion method" [8,9]:

$$\sigma^{\infty} + \sigma = C^{*}(\epsilon^{\infty} + \epsilon) = C(\epsilon^{\infty} + \epsilon - \epsilon^{*}) \quad \text{in } D \qquad (8)$$

where $\epsilon^{\infty}$ is the remote uniform strain field applied on the matrix $\Omega$ and $\sigma^{\infty}$ is the uniform stress obtained if there is no inhomogeneity inside the matrix, which means $\sigma^{\infty} = C\epsilon^{\infty}$, $\epsilon$ and $\sigma$ are the additional disturbed strain and stress respectively. $\epsilon^*$ is called eigenstrain and is related to the additional disturbed strain $\epsilon$:

$$\epsilon = S\epsilon^* \qquad (9)$$

$S$ is the Eshelby tensor, which is determined by the shape of the inhomogeneity. The details of the determination of the Eshelby tensor components $S_{ijkl}$ can be found in [18] and Appendix B lists the spherical case.

According to Eqs. (8) and (9), the analytic local stress field inside the inhomogeneity can be determined.

In this simulation, we assume the spherical inhomogeneity to be located at the center of the cube. In order to obtain a uniform stress field on the cube boundary, we set the ratio between the sphere radius $r$ and the cube length $L$ to be $r/L = 1/20$. The Young's modulus of the sphere is 10 times larger than that of the matrix. The matrix is subjected to a homogeneous traction strain $\epsilon_{zz}^{\infty} = 0.01$ on the boundary. Local grid refinement is applied and 409 points are used per sphere diameter. According to Eshelby's theory, the stress field produced in the inhomogeneity is uniform.
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Stress component comparison between the Eshelby solution and the MultiGrid solution.</figcaption></figure>

Fig. 3 shows a very good agreement for the inner points of sphere e.g. at the center $z = 0$:

$$ \sigma_{xx}^{\text{Eshelby}} = 0.582\sigma_{zz}^{\infty}, \quad \sigma_{xx}^{\text{MG}} = \sigma_{xx}^{\infty} + \sigma_{xx} = 0.584\sigma_{zz}^{\infty} $$

$$ \sigma_{yy}^{\text{Eshelby}} = 0.582\sigma_{zz}^{\infty}, \quad \sigma_{yy}^{\text{MG}} = \sigma_{yy}^{\infty} + \sigma_{yy} = 0.584\sigma_{zz}^{\infty} $$

$$ \sigma_{zz}^{\text{Eshelby}} = 1.662\sigma_{zz}^{\infty}, \quad \sigma_{zz}^{\text{MG}} = \sigma_{zz}^{\infty} + \sigma_{zz} = 1.665\sigma_{zz}^{\infty} $$

except for the adjacent points near the boundary between sphere and matrix, where these points are coupled with the points in the matrix.

## 3.2. Macroscopic: homogenization–Mori Tanaka

As pointed out above, the difficulty of simulating the real composite structure prompted researchers to develop a number of homogenization methods to predict the global stress–strain response [19]. A classic homogenization method is the Mori–Tanaka model.

In the Mori–Tanaka’s theory [10,11], the global effective property tensor $C^{\text{eff}}$ can be derived from:

$$ C^{\text{eff}} = C + V_f[(C^* - C)T] [(1 - V_f)I + V_f T]^{-1} \quad (10) $$

where $T$ is:

$$ T = [I + SC(C^* - C)]^{-1} \quad (11) $$

where $S$ is the Eshelby tensor, $I$ denotes the unit tensor, $V_f$ is the volume percent of the inhomogeneity, $C$ and $C^*$ are the stiffness tensor of matrix and inhomogeneity respectively.

<caption>**Table 1**</caption>

Comparison of elastic constants of the stiffness matrix between the Mori–Tanaka (MT) prediction and MultiGrid (MG) calculations.

<table><thead><tr><th colspan="7">(C<sup>eff</sup><sub>ijkl</sub>)<sub>MT</sub></th></tr></thead><tbody><tr><td>ijkl</td><td>11</td><td>22</td><td>33</td><td>12</td><td>23</td><td>13</td></tr><tr><td>11</td><td>4.2488</td><td>1.7899</td><td>1.7899</td><td>0</td><td>0</td><td>0</td></tr><tr><td>22</td><td>1.7899</td><td>4.2488</td><td>1.7899</td><td>0</td><td>0</td><td>0</td></tr><tr><td>33</td><td>1.7899</td><td>1.7899</td><td>4.2488</td><td>0</td><td>0</td><td>0</td></tr><tr><td>12</td><td>0</td><td>0</td><td>0</td><td>1.3158</td><td>0</td><td>0</td></tr><tr><td>23</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1.3158</td><td>0</td></tr><tr><td>13</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1.3158</td></tr></tbody><thead><tr><th colspan="7">(C<sup>eff</sup><sub>ijkl</sub>)<sub>MG</sub></th></tr></thead><tbody><tr><td>ijkl</td><td>11</td><td>22</td><td>33</td><td>12</td><td>23</td><td>13</td></tr><tr><td>11</td><td>4.2695</td><td>1.7852</td><td>1.7852</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr><tr><td>22</td><td>1.7852</td><td>4.2695</td><td>1.7852</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr><tr><td>33</td><td>1.7852</td><td>1.7852</td><td>4.2695</td><td>0.0000</td><td>0.0000</td><td>0.0000</td></tr><tr><td>12</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>1.2260</td><td>0.0000</td><td>0.0000</td></tr><tr><td>23</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>1.2260</td><td>0.0000</td></tr><tr><td>13</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>0.0000</td><td>1.2260</td></tr></tbody></table>

In this simulation, the model geometry is the same as in the Eshelby case but the ratio between the sphere radius $r$ and the cube length $L$ is selected to be $r/L = 1/4$, which gives $V_f = 6.54\%$. The Young’s modulus of the sphere is also 10 times larger than that of the matrix. Periodic boundary conditions are applied. Table 1 shows the results of the Mori–Tanaka prediction and the MG numerical results. Components of the stiffness tensor $C_{ijkl}$ in traction part ($C_{ijjj}$) match very well where the differences are 0.49% and 0.26%, while for the shear modulus ($C_{ijjj}$, $i \neq j$) the difference is 6.82%. Since the structure studied is isotropic, the equivalent shear modulus $\mu_{\text{equivalent}}^{\text{MT}}$ calculated in MT prediction according to the other two components ($C_{1111}^{\text{MT}}$ and $C_{1122}^{\text{MT}}$) is:

$$ \mu_{\text{equivalent}}^{\text{MT}} = \frac{C_{1111}^{\text{MT}} - C_{1122}^{\text{MT}}}{2} \quad (12) $$

Its value $\mu_{\text{equivalent}}^{\text{MT}} = 1.2295$, is very close to the MG result 1.2260 where the difference becomes 0.28%.

## 4. Application

After the validation of both microscopic and macroscopic aspects, the MG solver is used to simulate the composite structure. In this section, the solver simulates the fiber-reinforced laminated structure. This kind of structure is widely used since it is easy to obtain the specific mechanical properties by modifying the orientation of the fibers and the stacking sequence. However, it is also well known that there are free edge effects induced by the mismatch of adjacent layers [20]. This mismatch can lead to stress concentrations at the interface between two different layers [21], which may induce the initial inter-laminar failures such as delamination or matrix cracking. A displacement discontinuity on the boundary surfaces observed in residual deformation states can be interpreted [22] as a signature of this early material damage.

In this study, the MG model is first applied to an ideal fiber reinforced structure case to obtain the effective properties. Then we perform the numerical calculations with both MG and FE methods. The MG routine uses the ideal laminated structure while the FE routine uses the effective properties obtained from MG homogenization. After the comparison between these two methods, we introduce the real structure into the model, which is obtained using X-ray tomography. We compare this case with a similar ideal case which has perfect cylindrical fibers. Then we study the influence of fiber orientation, interface layer thickness, fiber layer thickness as well as the different ratios of the Young’s modulus between fibers and matrix.

### 4.1. Comparison with FE

In this part, the simulated structure contains two layers of different fiber orientations with angles [−15°/+15°] to the z axis. In each layer, the fibers are staggered, see Fig. 4. The bulk dimensions compared to the fiber radius $r$ are [±8r, ±8r, ±32r]. The fibers have a Young’s modulus 10 times larger than the matrix denoted as $E_f/E_m = 10$. In each layer, the vertical distance between any two fibers is 2.2r.

First we apply a homogeneous strain boundary condition (HSBC) to obtain the effective properties of this structure with the MG method. In this step, we can also select a representative volume element (RVE) to be analyzed, but a large scale dimension will be more convincing considering the latter real case. Then in the FE routine, we use the homogenized results under the following conditions: at the bottom surface ($z = -32r$): $u = v = w = 0$, at the top surface ($z = 32r$): $u = v = 0$, $w = 0.01 \times 64r$, and on the other four surfaces the normal displacement is kept zero. These conditions are chosen because they approximate the actual free
<figure><img src="image_5.png" /><figcaption>Fig. 4. Laminated structure composed of two layers with different fiber orientations (top layer +15°, bottom layer -15°).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6. Displacement w at the intersection between the two planes $y = -8r, z = 0$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Displacement w distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Shear stress $\tau_{xz}$ distribution (left FE result with two different homogeneous layers, right MG result with two different fiber orientation layers).</figcaption></figure>
edge conditions. In the MG routine, we apply the same boundary conditions while using the ideal laminated structure (Fig. 4).

Fig. 5 shows the displacement w distribution on the boundary surfaces, where the left one shows the results of the FE calculation using a homogeneous structure and the right one shows the results of the MG calculation using the ideal laminated structure. A

displacement variation along the x direction is found at the plane $y = -8r$, the plane visible in Fig. 5, because of the mismatch induced by the different fiber orientations. If we select a line on this plane, such as $z = 0$ i.e. in the middle section of the sample, Fig. 6 shows a good agreement of the local displacement w distribution along the x direction between FE and MG calculations. Obviously, the displacement w distribution in FE results is smooth due to the homogenization while in the MG results the distribution varies according to the fiber location.

The variation of the displacement w near the interface induces a stress concentration in $\tau_{xz}$. Fig. 7 shows the shear stress $\tau_{xz}$ distribution on the boundary surfaces. The shear stress $\tau_{xz}$ shows concentration where two fibers with different orientations are close. Fig. 8 shows the comparison of the local shear stress $\tau_{xz}$ at two different positions in the interface $x = 0$, which shows that $\tau_{xz}$ increases sharply when approaching to the boundary. A good agreement between FE and MG is found for both positions while MG results shows the detailed local stress distribution.

From the above comparisons, we find that the global effects of both components match very well between FE and MG calculations. However, the MG simulation allows one to capture all the strain and stress heterogeneities generated by the micro-structure.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Shear stress $\tau_{xz}$ comparison at the interface $x = 0$.</figcaption></figure>

## 4.2. Real case

In this simulation, the real structure (Fig. 9) captured using X-ray tomography is modeled. These data were obtained in [22]. Because of the large number of points describing the structure and the associated computational cost, it is difficult to use traditional numerical methods, such as FE. The real structure contains three main phases: one for the fibers, one for the matrix and the other one for the voids. The ratio of the Young's modulus between them is $E_f = 10E_m = 100E_\nu$. It has two layers with different fiber orientation with angles [−15°/+15°] as analyzed previously. The boundary conditions are also applied as in the previous case. The coarse grid is meshed with [9, 9, 33] points and 6 levels in total are used with [257, 257, 1025] points on the finest grid. Each fiber has at least 20 points along its diameter.

Since the random nature exists in the real structure case and its interface layer thickness is between $2r$ and $3r$, the thickness of interface layer in two similar ideal cases with perfect cylindrical fibers and definite interface layer thickness are selected to be $2r$ and $3r$ respectively. Fig. 10 shows the structure and displacement w distribution at the boundary surface: $y = -32r$ for both real structure case and one ideal structure case (interface layer

<figure><img src="image_10.png" /><figcaption>Fig. 9. Measured structure: yellow-fibers, light blue-matrix, dark blue-voids (top layer +15°, bottom layer -15°). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Structure and displacement w distribution at the boundary surface: $y = -32r$ for the real case and the ideal case (interface layer thickness = $2r$) respectively (from left to right: real structure, displacement w distribution for the real case, ideal structure, displacement w distribution for the ideal case).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11. Displacement $\Delta w$ comparison between the real case and the ideal case.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Influence of fiber orientation on $\Delta w$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Influence of interface layer thickness on $\Delta w$.</figcaption></figure>

thickness = $2r$). The variation of displacement $w$ along the x direction can also be observed in the real structure case.

Fig. 11 shows the comparison of $\Delta w$ between the real case and two similar ideal cases, where $\Delta w$ is the average displacement $w$ along the z direction minus its average along the x direction and describes the global displacement variation. These results are very similar.

### 4.3. Parameter study for the ideal case

In this section, we investigate several parameters influencing the $\Delta w$.

The first parameter is the fiber orientation. We model three different cases: [-10°/+10°], [-15°/+15°] and [-20°/+20°]. In all these cases, the thickness of the interface layer and fiber layer is kept at

<figure><img src="image_15.png" /><figcaption>Fig. 14. Influence of fiber layer thickness on $\Delta w$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Influence of the ratio of the Young's modulus on $\Delta w$.</figcaption></figure>

$2r$ and $6r$ respectively and the ratio of the Young's modulus between fiber and matrix is 10. Fig. 12 shows that the displacement variation $\Delta w$ increases with the angle. This confirms what is reported in the literature and indicates that a more severe mismatch of the properties between two layers yields a larger displacement variation. The mismatch depends on the angle and there is actually a maximum mismatch for an angle between 0° and 45°. The value of the angle for the maximum mismatch depends on the mechanical properties of the individual plies and should be obtained from the homogenized FE model by varying the orientation between the two plies. However, this investigation is beyond the scope of the paper.

The second parameter is the thickness of the interface layer. Four different thicknesses of the interface layer are simulated: 0r where two different fiber orientation layers touch each other at the interface, 1r, 2r and 5r. In all these cases, the thickness of the fiber layer is kept at 6r, the fiber orientation is set at [-15°/+15°] and the ratio of the Young's modulus between fiber and matrix is 10. Fig. 13 shows that a thicker interface has a larger amplitude of the displacement variation but has a smaller slope. The amplitude of the displacement variation between the 2r case and the 5r case varies little while the slope varies much, which indicates that the 2r case is more likely to induce the delamination than the 5r case from this point of view.

The third parameter is the thickness of the fiber layer. Two different thicknesses of fiber layer are simulated: 5.5r and 8r. The thickness of the interface layer is kept to be 5r, the fiber orientation is set at [-15°/+15°] and the ratio of the Young's modulus between fiber and matrix is 10. Fig. 14 shows that the laminated structure has the same displacement variation for different thicknesses of the fiber layer.
The fourth parameter is the ratio of the Young's modulus between fibers and matrix. In this comparison, the thickness of interface layer and fiber layer is kept at 2r and 6r respectively and the fiber orientation is set at [-15°/+15°]. Fig. 15 shows that the tendency of the displacement variation is steeper and larger when increasing the ratio of the Young's modulus, which indicates that a stronger discontinuity between fiber and matrix is more likely to induce delamination.

## 5. Conclusion

The efficiency of Multigrid techniques allows solution of the problem with dense grids and thus very local description of actual material topologies. This offers great opportunities for future research for material characterization, derivation of large scale (FEM) element models from local scale simulations and material optimization. The main work in this paper can be concluded as follows:

(1) An efficient and stable MG solver is further developed to solve the 3D Lamé equations for composite materials. The solver adopted the method of Alcouffe et al. to solve problems involving large material gradients. Local grid refinement techniques are applied to reduce the computational time and allocated memory.

(2) The solver is validated with respect to Eshelby's solution and the Mori-Tanaka prediction. A good agreement is found. Finally, the solver is applied to fiber reinforced materials, including theoretical material distributions and a measured distribution. These models use as much as $10^8$ DOFs and are computed efficiently. Parameters influencing the displacement variation on the boundary surfaces are studied.

(3) Fiber orientation, interface layer thickness and Young's modulus ratio are shown to be important parameters concerning the displacement variation. It is obtained that the fiber layer thickness does not influence the displacement variation. Conversely, it is shown that the interface layer thickness has a strong influence.

(4) A good agreement is found for the MG homogenization, which shows the potential of MG methods for the homogenization process.

## Acknowledgments

The first author would like to express heartfelt thanks to Dr. Boffy for his helpful directions during the author's study on Multi-Grid at the University of Twente. Financial support for this research was provided by the Chinese Scholarship Council (CSC).

## Appendix A. Discretize equations

The first equation in (1) can be written as:

$$ \underbrace{\frac{\partial}{\partial x}((\lambda + 2\mu)\frac{\partial u}{\partial x})}_{(1)} + \underbrace{\frac{\partial}{\partial y}(\mu\frac{\partial u}{\partial y})}_{(2)} + \underbrace{\frac{\partial}{\partial z}(\mu\frac{\partial u}{\partial z})}_{(3)} + \underbrace{\frac{\partial}{\partial x}(\lambda\frac{\partial v}{\partial y})}_{(4)} + \underbrace{\frac{\partial}{\partial y}(\mu\frac{\partial v}{\partial x})}_{(5)} + \underbrace{\frac{\partial}{\partial z}(\lambda\frac{\partial w}{\partial z})}_{(6)} + \underbrace{\frac{\partial}{\partial x}(\mu\frac{\partial w}{\partial x})}_{(7)} = 0 $$

Using second order accurate discretization:

$$ (1) \frac{(\lambda + 2\mu)_{i+1/2,j,k} u_{i+1,j,k} - ((\lambda + 2\mu)_{i+1/2,j,k} + (\lambda + 2\mu)_{i-1/2,j,k}) u_{i,j,k} + (\lambda + 2\mu)_{i-1/2,j,k} u_{i-1,j,k}}{h_x^2} $$

$$ (2) \frac{\mu_{i,j+1/2,k} u_{i,j+1,k} - (\mu_{i,j+1/2,k} + \mu_{i,j-1/2,k}) u_{i,j,k} + \mu_{i,j-1/2,k} u_{i,j-1,k}}{h_y^2} $$

$$ (3) \frac{\mu_{i,j,k+1/2} u_{i,j,k+1} - (\mu_{i,j,k+1/2} + \mu_{i,j,k-1/2}) u_{i,j,k} + \mu_{i,j,k-1/2} u_{i,j,k-1}}{h_z^2} $$

$$ (4) \frac{\lambda_{i+1,j,k}(v_{i+1,j+1,k} - v_{i+1,j-1,k}) - \lambda_{i-1,j,k}(v_{i-1,j+1,k} - v_{i-1,j-1,k})}{4h_x h_y} $$

$$ (5) \frac{\mu_{i,j+1,k}(v_{i+1,j+1,k} - v_{i-1,j+1,k}) - \mu_{i,j-1,k}(v_{i+1,j-1,k} - v_{i-1,j-1,k})}{4h_x h_y} $$

$$ (6) \frac{\lambda_{i+1,j,k}(w_{i+1,j,k+1} - w_{i+1,j,k-1}) - \lambda_{i-1,j,k}(w_{i-1,j,k+1} - w_{i-1,j,k-1})}{4h_x h_z} $$

$$ (7) \frac{\mu_{i,j,k+1}(w_{i+1,j,k+1} - w_{i-1,j,k+1}) - \mu_{i,j,k-1}(w_{i+1,j,k-1} - w_{i-1,j,k-1})}{4h_x h_z} $$

## Appendix B. Eshelby's tensor components for a spherical inhomogeneity

$$ S_{1111} = S_{2222} = S_{3333} = \frac{7 - 5\nu}{15(1 - \nu)} $$

$$ S_{1122} = S_{2233} = S_{3311} = S_{1133} = S_{2211} = S_{3322} = \frac{5\nu - 1}{15(1 - \nu)} $$

$$ S_{1212} = S_{2323} = S_{3131} = \frac{4 - 5\nu}{15(1 - \nu)} $$

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

## References

[1] Z. Hashin, J. Appl. Mech. 46 (3) (1979) 543–550.

[2] A. Brandt, 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 in Springer, 1973, pp. 82–89.

[3] A. Brandt, Mathematical Software III (1977) 273–314.

[4] C. Venner, A. Lubrecht, Multi-level Methods in Lubrication, vol. 37, Elsevier, 2000.

[5] B. Watremetz, M.C. Baietto, A. Lubrecht, Tribol. Int. 40 (5) (2007) 754–762.

[6] H. Boffy, M.C. Baietto, P. Sainsot, A. Lubrecht, J. Tribol. 134 (2) (2012) 021401.

[7] H. Boffy, C. Venner, Tribol. Int. 74 (2014) 121–129.

[8] J. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 241, The Royal Society, 1957, pp. 376–396.

[9] J. Eshelby, Prog. Solid Mech. 2 (1) (1961) 89–140.

[10] T. Mori, K. Tanaka, Acta Metall. 21 (5) (1973) 571–574.

[11] Y. Benveniste, Mech. Mater. 6 (2) (1987) 147–157.

[12] H. Boffy, Techniques multigrilles et raffinement pour un modèle 3d efficace de milieux hétérogènes sous sollicitations de contact, Ph.D. Thesis, INSA de Lyon, 2012.

[13] R. Alcouffe, A. Brandt, J. Dendy Jr., J. Painter, SIAM J. Sci. Stat. Comput. 2 (4) (1981) 430–454.

[14] R. Hoekema, C. Venner, J. Struijk, J. Holsheimer, Comput. Biomed. Res. 31 (5) (1998) 348–362.

[15] Z. Hashin, B. Rosen, J. Appl. Mech. 31 (2) (1964) 223–232.

[16] K. Terada, M. Hori, T. Kyoya, N. Kikuchi, Int. J. Solids Struct. 37 (16) (2000) 2285–2311.

[17] Z. Xia, Y. Zhang, F. Ellyin, Int. J. Solids Struct. 40 (8) (2003) 1907–1921.

[18] T. Mura, Micromechanics of Defects in Solids, vol. 3, Springer Science & Business Media, 2012.

[19] B. Klusemann, H. Böhm, B. Svendsen, Eur. J. Mech. A/Solids 34 (2012) 21–37.

[20] R. Pipes, N. Pagano, J. Compos. Mater. 4 (4) (1970) 538–548.

[21] A. Nosier, M. Maleki, Int. J. Mech. Sci. 50 (10) (2008) 1435–1447.

[22] P. Lecomte-Grosbras, J. Réthoré, N. Limodin, J. Witz, M. Brieu, Exp. Mech. 55 (1) (2015) 301–311.