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

• An efficient Multigrid solver for simulating composites.
• A real topology is simulated for the first time by the MultiGrid methods.
• MultiGrid application in composite homogenization process.
• Factors influencing on surface displacement variation are studied.

The subject of this paper is to present the application of MultiGrid (MG) methods in 3D composite material 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 laminated 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.


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 microstructure, 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 homogenization 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 various scientific and engineering fields. Lately a number of researchers 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) problems. 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 techniques 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 inhomogeneity. 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 ⇑ Corresponding author.
E-mail address: hanfeng.gu@insa-lyon.fr (H. Gu). 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 effect the delamination process are investigated as well.

Theory
The constitutive equations to be solved are the 3D linear elastic equations expressed by the displacements u; v and w in Cartesian coordinates: where r is the gradient operator: and F i is the body force. In this study, we assume it to be zero and X is the 3D rectangular domain. Eq. (1) describes the static equilibrium in X, while on the boundary @X displacement boundary conditions and stress boundary conditions can both be considered.

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 X 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  This scheme can be used recursively for any number of levels.

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

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: 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.

Boundary conditions
According to the problem studied, the boundary conditions can be divided into 3 types: displacements specified on @X, stress specified on @X or displacements specified on a portion @X 1 and stress specified on a portion @X 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 @X is: If the applied strain ij is given, then all the displacements on @X can be specified according to Eq. (6).
For a periodic RVE, PBC is more suitable [16,17]. The boundary condition on @X is: where u Ã i is the unknown periodic displacement on the boundary @X and is dependent on the given global strain field ij . Hence for a given strain ij , the displacement u i on @X can also be determined by Eq. (7).

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 X 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.

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]: where 1 is the remote uniform strain field applied on the matrix X and r 1 is the uniform stress obtained if there is no inhomogeneity inside the matrix, which means r 1 ¼ C 1 . and r are the additional disturbed strain and stress respectively. Ã is called eigenstrain and is related to the additional disturbed strain : 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 1 zz ¼ 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.

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 eff can be derived from: where T is: 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.
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 iijj ) match very well where the differences are 0.49% and 0.26%, while for the shear modulus (C ijij ; ij) the difference is 6.82%. Since the structure studied is isotropic, the equivalent shear modulus l MT equivalent calculated in MT prediction according to the other two components (C MT 1111 and C MT 1122 ) is: Its value l MT equivalent ¼ 1:2295, is very close to the MG result 1.2260 where the difference becomes 0.28%.

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.

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 ½AE8r; AE8r; AE32r. 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):     Fig. 6. Displacement w at the intersection between the two planes y ¼À8r; z ¼ 0. 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 s xz . Fig. 7 shows the shear stress s xz distribution on the boundary surfaces. The shear stress s xz shows concentration where two fibers with different orientations are close. Fig. 8 shows the comparison of the local shear stress s xz at two different positions in the interface x ¼ 0, which shows that s 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.

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 v . 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  thickness = 2r). The variation of displacement w along the x direction can also be observed in of the real structure case. Fig. 11 shows the comparison of Dw between the real case and two similar ideal cases, where Dw 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.

Parameter study for the ideal case
In this section, we investigate several parameters influencing the Dw.
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 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 Dw 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 s 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 theses 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.

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.