1 - Introduction

Summary

This chapter provides a very brief summary of the types of heterogeneous materials considered in this monograph: fiber-reinforced composites, particulate composites, nanocomposites, porous composites, and so on. A succinct summary is given of analytical homogenization methods to determine the overall properties of particulate composites based on the upper and lower bounds of Hashin and Shtrikman; the Eshelby ellipsoidal inclusion theory and the Self-consistent Method of Eshelby; and the Mori-Tanaka Method and some other semi-analytical methods. Numerical methods such as the finite element method, the boundary element method, XFEM, and so on to model a representative volume element (RVE) of a heterogeneous material are reviewed, and thus the motivation for the Computational Grains method discussed in the rest of this book is presented.

Summary

This chapter provides a very brief summary of the types of heterogeneous materials considered in this book: fiber-reinforced composites, particulate composites, nanocomposites, porous composites, and so on. A succinct summary is given of analytical homogenization methods to determine the overall properties of particulate composites based on the upper and lower bounds of Hashin and Shtrikman, the Eshelby ellipsoidal inclusion theory and the self-consistent method of Eshelby, the Mori–Tanaka method, and some other semi-analytical methods. Numerical methods such as the finite element method (FEM), the boundary element method (BEM), extended finite element method (XFEM), and so on to model a representative volume element (RVE) of a heterogeneous material are reviewed. The chapter also presents the motivation for the Computational Grains (CGs) method discussed in the rest of this book.

1.1 Heterogeneous Materials and Their Applications

Many factors catalyzed the progressive development of modern composite materials. Since the 1930s, the resin materials have been the subject of extensive academic research and industrial applications. Due to the low modulus of strength of pure resin materials, glass fibers were used to reinforce the weaker resin materials to increase their mechanical properties. In 1936, the first generation of glass fiber-reinforced plastics (GFRPs) was introduced. Due to their advantages of high performance and low price, GFRPs still capture a huge market share in material industry.

In the 1960s, Texaco developed high-strength and high-modulus boron fibers through chemical vapor processes. A boron fiber has a relatively large diameter of about 100–200 μm, and it is very brittle and sensitive to surface damage. Compared with boron fibers, carbon fibers, which appeared at the same time, can be manufactured through simpler processes. The modulus of modern carbon fiber is about four times higher than that of steel, while the density is about four times lower than that of steel [1]. Due to such advantages of a simple manufacturing method, low cost, and high performance, carbon fibers have become some of the most widely used high-performance fiber materials up to now (Figure 1.1).

(a) carbon fiber and

(b) boron fiber [2]

Figure 1.1 SEM scan of

In addition to inorganic fibers, high-performance polymer fiber materials have also been rapidly developed. From the 1960s to the 1980s, high-performance polymer fibers such as Kevlar and Zylon, which had many excellent mechanical and thermal properties (low thermal conductivity, good thermal insulation, excellent toughness, and corrosion resistance), were developed. The high-performance polymer fiber materials are widely used nowadays in many fields, such as automobiles, ships, digital electronic equipment, sports equipment, body armor, aerospace, and so on [3].

Applications of materials in very high-temperature environments necessitated the development of metal matrix composites (MMCs) and ceramic matrix composites (CMCs). Metal matrix composites add particles/fibers (such as ceramics, carbon fibers, etc.) and other enhanced phases to metal or alloy matrices to increase the heat resistance, strength, and rigidity of light metals and reduce their thermal expansion coefficients [4]. Currently, the MMC technology has become increasingly mature, and is widely used in transportation, aerospace, sports equipment, and many other fields.

Compared with MMCs, CMCs have higher strength and melting points, and more stable chemical properties. In the past thirty years, advances in technologies for improving the fracture toughness of the CMCs mitigated the problem of brittleness of CMC materials, allowing them to be widely used in cutting tools, in hot sections of engines, as brake materials, and in spacecraft insulation systems (Figure 1.2).

(a) MMCs [5] and

(b) CMCs [6]

Figure 1.2 SEM scan of

Since the 1990s, composite materials reinforced with nano-sized constituents have gradually become more common. “Nanocomposites” refer to composite materials containing one or more constituent materials with a size of 1–100 nm, and their reinforcement phases consisting generally of nanoparticles (nano SiC, nano Al₂O₃), nanoflakes (graphite), and nanofibers (carbon nanotubes). Unlike conventional composite materials, the ratio of area to volume or the aspect ratio of nanocomposite reinforcement phase is very large, and therefore its physical properties can be engineered very significantly. By studying the effect of the size of the reinforcing-phase material on the mechanical properties of the composite material from an experimental perspective in detail, it has been shown that, through appropriate manufacturing processes, nano-sized reinforcing-phase materials can effectively improve the mechanical properties, thermal properties, and other physical properties of composite materials [7].

In addition to the aforementioned composite materials used for strength, stiffness, and fracture toughness-based structural applications, biomimetic composite materials, multifunctional composite materials, and other types of novel composite materials have been vigorously developed for various applications.

1.2 Analytical and Numerical Micromechanics of Heterogeneous Materials

As can be seen from the previous section, composite materials have been developed rapidly in modern times, and have become one of the most important class of structural as well as multifunctional materials. Studies of micromechanics of composite materials play a vital role in the study of stiffness, strength, and durability of such materials; so it is of great significance to establish efficient and accurate tools for the micromechanical analyses of composite materials. This chapter briefly introduces the status of the subject of micromechanics of composite materials through analytical as well as numerical methods.

1.2.1 Analytical Micromechanics Methods

The study of micromechanics of heterogeneous materials usually involves a representative volume element (RVE) for analysis. From the macroscale perspective, an RVE can be regarded as an infinitesimal point, and its macro-stress and macro-strain are uniform. From the microscale perspective, an RVE must represent a large enough volume so as to contain sufficient microstructural information, and delineate the details of the internal microscopic stress and strain fields, which are generally nonuniform. By selecting an appropriate RVE and then using analytical micromechanical methods, we can obtain the macroscopic mechanical parameters of the material. The commonly used semi-analytical micromechanical methods are (i) variational methods [8–12], (ii) homogenization methods based on the Eshelby inclusion theory [18], and (iii) differential methods [19].

The commonly used variational methods include the Voigt–Reuss methods of bounds [8], and the Hashin–Shtrikman variational method [12].

Voigt applied the condition $\overline{u}=\overline{\epsilon }\cdot x$ in the entire domain of the RVE, where $\overline{u}$ denotes the displacement vector, $\overline{\epsilon }$ denotes the applied strain tensor, and $x$ denotes the position vector. This is to assume that the strain field in the reinforcing-phase material and in the matrix are equal to the applied strain $\overline{\epsilon }$ , and the equivalent elasticity tensor ${\text{C}}_{\text{v}}$ of the composite material can be derived for isotropic matrix and isotropic inclusions as:

$C_{\text{v}}=(1-{\text{f}}_{\text{m}})C_{\text{m}}+\sum _{\alpha =1}^n{\text{f}}_{\alpha }{C}_{\alpha }$(1.1)

where $C_{\text{m}}$ and $C_{\alpha }$ denote the elasticity tensors of the matrix and the $\alpha {\mathrm{th}}^{}$ inclusion of the heterogeneous material, respectively, ${\text{f}}_{\text{m}}$ and ${\text{f}}_{\alpha }$ denote the volume fractions of the matrix and the $\alpha \mathrm{th}$ inclusion, respectively; and these notations are used consistently throughout this chapter. As summarized in [9], for the symmetric elasticity tensor $C$ of the isotropic materials, transversely isotropic materials, orthotropic materials, and general anisotropic materials, there are, respectively, two, five, nine, and twenty-one independent constants. According to the principle of minimum complementary energy [8], $C_{\text{v}}$ is an upper bound for the true equivalent elasticity tensor of composite materials. Reuss [10] applied surface boundary conditions $\overline{t}=\overline{\sigma }\cdot n$ and assumed that the average stress of the reinforcing-phase material and the matrix is equal to the applied stress $\overline{\sigma }$ , and showed that the equivalent elasticity tensor $C_{\text{R}}$ of the composite material can be derived for isotropic matrix and isotropic inclusions as:

$C_{\text{R}}={\left((1-f_{\text{m}})C_{\text{m}}^{\text{-1}}+\sum _{\alpha =1}^n{\text{f}}_{\alpha }{C}_{\alpha }^{\text{-1}}\right)}^{-1}$(1.2)

According to the principle of minimum potential energy [11], $C_{\text{R}}$ is a lower bound for the equivalent elasticity tensor of the composite material.

However, there is a large difference between the upper bound of Voigt and the lower bound of Reuss, when the elastic moduli of the matrix material and those of the reinforcement phase materials differ substantially. Later, Hashin and Shtrikman [12] resolved this problem by using variational principles and extremum methods for the first time, leading to what are popularly known as the Hashin–Shtrikman variational bounds. Since the aforementioned methods did not take into account the internal microstructure of the composite material, it is generally used only to make a general approximate prediction of the equivalent moduli of the material, but cannot provide information on the detailed stress and strain fields within the RVE.

Another set of analytical micromechanics methods is based on the celebrated work of J. D. Eshelby. By employing the Kelvin’s fundamental solution for an isotropic material Eshelby [13] gave an analytical solution for an infinite medium containing an ellipsoidal inclusion which is subjected to a uniform characteristic strain field (eigenstrain) ${\epsilon }^{*}$ , which is briefly summarized in here.

First consider the problem of an eigenstrain field ${\epsilon }^{*}\left(X\right)$ prescribed in the ellipsoidal inclusion $W$ and vanishing outside, where the isotropic elasticity tensor for the entire solid medium is $C$ [14]:

${\epsilon }^{*}\left(\text{x}\right)=\{\begin{array}{c}{\epsilon }^{*}\left(\text{x}\right),\text{ for }\text{x}\text{i}\text{n}\text{W}\\ 0\text{ , }\text{o}\text{t}\text{h}\text{e}\text{r}\text{w}\text{i}\text{s}\text{e}\end{array}$(1.3)

Let $u^e$ be the displacement field produced in the infinite isotropic medium by ${\epsilon }^{*}$ , and denote by ${\epsilon }^e$ and ${\sigma }^e$ the corresponding strain and stress fields. These strain and stress fields are expressed as:

${\epsilon }^e=\frac{1}{2}\left\{\left(\nabla u^e\right)+{\left(\nabla u^e\right)}^T\right\}$(1.4)

${\sigma }^e=C:{\epsilon }^e-C:{\epsilon }^{*}$(1.5)

We know that the displacement field $u$ produced by a distributed body force $f$ should be subject to the following governing equation:

$\nabla \cdot \left(C:\nabla \otimes \text{u}\left(\text{x}\right)\right)+\text{f}\left(\text{x}\right)=0$(1.6)

By substituting eq. (1.4) and eq. (1.5) into the equations of equilibrium of a solid body with negligible body force, we get:

$\nabla \cdot {\text{σ}}^e=\nabla \cdot \left(C:\nabla \otimes u^e\right)+\nabla \cdot \left(-C:{\epsilon }^{*}\right)=0$(1.7)

wherein the term $\nabla \cdot \left(-C:{\epsilon }^{*}\right)$ acts like a distributed body force in the governing equations for the displacement field, shown as eq. (1.6).

The displacement field $u^e$ corresponding to the distributed body forces as in eq. (1.7) can be expressed in terms of the Green’s function $G$ (Kelvin’s solution for an isotropic material), as:

$u_i^e\left(x\right)={\int }_W{G}_{ki.l}^{\infty }\left(y-x\right)C_{klmn}{\epsilon }_{mn}^{*}\left(y\right)\text{d}{\text{V}}_{\text{y}}$(1.8)

where ${\left(\right)}_{.l}$ denotes a derivative with respect to $x_l$ . The Green’s function $G_{ij}^{\infty }\left(x-y\right)$ from the Kelvin solution shown in eq. (1.8) gives the displacement component in the ${\text{x}}_i$ -direction at point $x$ , produced by a unit point force applied in the ${\text{x}}_j$ -direction at point $y$ , which can be written as:

$G_{ij}^{\infty }\left(x-y\right)=\frac{1}{4\pi \mu }\frac{{\delta }_{ij}}{\left|x-y\right|}-\frac{1}{16\pi \mu (1-\nu )}\frac{{\partial }^2}{\partial x_i\partial x_j}\left|x-y\right|$(1.9)

where $\mu$ and $\nu$ denote the Lamé constant and the Poisson’s ratio of the isotropic material, $\left|x-y\right|$ denotes the distance between the point $x$ and point $y$ , and ${\delta }_{ij}$ denotes the Kronecker delta.

The strain field ${\epsilon }_{ij}^{}\left(x\right)$ in the infinite isotropic medium with an isotropic elasticity tensor $C_{mnkl}$ thereafter becomes:

${\epsilon }_{ij}\left(x\right)=S_{ijkl}^{\infty }\left(x;W\right){\epsilon }_{kl}^{*}$(1.10)

where the tensor field $S_{ijkl}^{\infty }\left(x;W\right)$ is referred to as the Eshelby tensor, derived as:

$S_{ijkl}^{\infty }\left(x;W\right)\equiv \left\{{\int }_W{\Gamma }_{ijmn}^{\infty }\left(x-y\right)\text{d}{\text{V}}_{\text{y}}\right\}{C}_{mnkl}$(1.11)

The tensor field ${\Gamma }_{}^{\infty }$ is derived for the Green’s function in eq. (1.9), which is limited to only an isotropic material:

${\Gamma }_{ijkl}^{\infty }=-\frac{1}{4}\left\{G_{ik.jl}^{\infty }+G_{jk.il}^{\infty }+G_{il.jk}^{\infty }+G_{jl.ik}^{\infty }\right\}$(1.12)

where ${\left(\right)}_{.jk}$ denotes the second-order derivative with respect to $x_j$ and $x_k$ .

It can be shown that the strain inside the inclusion is uniform, unlike that outside of the inclusion. On this basis, homogenization methods based on the Eshelby inclusion theory, for isotropic matrices and isotropic inclusions only, have been widely developed. The commonly used ones are as follows.

The dilute approximation [13,15] assumes that the strain field of each inclusion in the RVE is equal to the strain field due to a single inclusion in an infinite solid and employs the Eshelby equivalent inclusion method, which is to consider the elasticity tensor of the isotropic inclusion $C_c$ as the elasticity tensor of the isotropic matrix $C_m$ with the eigenstrain field ${\epsilon }_{}^{*}$ , as shown as Figure 1.3.

Figure 1.3 Schematic diagram of the Eshelby equivalent inclusion method [ $C_m$ is the isotropic elasticity tensor of the matrix, $C_c$ is the isotropic elasticity tensor of the inclusion, and ${\epsilon }_{}^{*}$ is the eigenstrain field shown in eq. (1.3)]

The strain in the inclusions can be written as:

${\epsilon }_{\alpha }={\left[I+S_{\alpha }:C_{\text{m}}^{-1}:\left(C_{\alpha }-C_{\text{m}}\right)\right]}^{-1}:\overline{\epsilon }$(1.13)

where ${\epsilon }_{\alpha }$ denotes the strain field due to the $\alpha \mathrm{th}$ inclusion, $S_{\alpha }$ denotes the Eshelby tensor, $C_{\text{m}}^{}$ and $C_{\alpha }^{}$ denote the isotropic elasticity tensors of the matrix and the $\alpha \mathrm{th}$ inclusion, respectively, and $\overline{\epsilon }$ denotes the applied uniform strain. Therefore, we can obtain the equivalent elasticity tensor as [15]:

$C_{\text{dilute}}^{}=C_{\text{m}}+\sum _{\alpha =1}^N{\text{f}}_{\alpha }\left(C_{\alpha }-C_{\text{m}}\right){\left[I+S_{\alpha }:C_{\text{m}}^{-1}:\left(C_{\alpha }-C_{\text{m}}\right)\right]}^{-1}$(1.14)

where $f_{\alpha }$ denotes the volume fraction of the $\alpha \mathrm{th}$ inclusion. Equation (1.14) has been derived mainly for an isotropic matrix with an elasticity tensor $C_m$ , and isotropic inclusions with isotropic elasticity tensors $C_{\alpha }$ .

The dilute approximation shown in eq. (1.14) ignores the arbitrary distribution of inclusions and the interactions between the inclusions, and is only suitable for the case of inclusions with small-volume fractions. The Mori–Tanaka method [16] considers such interactions between inclusions and shows that the strains in the inclusions should be written as:

${\epsilon }_{\alpha }={\left[I+S_{\alpha }:C_{\text{m}}^{-1}:\left(C_{\alpha }-C_{\text{m}}\right)\right]}^{-1}:{\epsilon }_{\text{m}}$(1.15)

where ${\epsilon }_{\text{m}}$ denotes the strain in the matrix. By using the equation $〈\epsilon 〉{\text{=f}}_{\text{m}}〈{\epsilon }_{\text{m}}〉\text{+}\sum _{\alpha \text{=1}}^N{\text{f}}_{\alpha }〈{\epsilon }_{\alpha }〉$ , we can find the strain concentration tensor of each inclusion:

$A_{\alpha }=T_{\alpha }:{\left[{\text{f}}_{\text{m}}{T}_{\text{m}}+\sum _{n=1}^N{\text{f}}_n{T}_n\right]}^{-1}$(1.16)

where $T_n\text{=}{\left[I+S_n:C_{\text{m}}^{-1}:\left(C_n-C_{\text{m}}\right)\right]}^{\text{-1}}$ . By eq. (1.16), the equivalent elastic tensor of the composite can be obtained as:

$C_{\text{MT}}^{}=C_{\text{m}}+\sum _{\alpha =1}^N{\text{f}}_{\alpha }\left(C_{\alpha }-C_m\right):T_{\alpha }:{\left[{\text{f}}_{\text{m}}{T}_{\text{m}}+\sum _{n=1}^N{\text{f}}_n{T}_n\right]}^{-1}$(1.17)

The Mori–Tanaka method is generally applicable to the case where the volume fractions of the inclusions are less than 30 percent.

The self-consistent method [18] assumes that the inclusions are buried in the reference material, whose elastic tensor is equal to the equivalent elastic tensor of the composite material, as shown in Figure 1.4(a).

(a) self-consistent method and

(b) generalized self-consistent method

Figure 1.4 Model of

The strain in each inclusion can be expressed as:

${\epsilon }_{\alpha }={\left[I+{\overline{S}}_{\alpha }:{\left(C_{\text{sc}}^{}\right)}^{-1}:\left(C_{\alpha }^{}-C_{\text{sc}}^{}\right)\right]}^{-1}:E$(1.18)

where ${\overline{S}}_{\alpha }$ denotes the Eshelby tensor of the composite material and $C_{\text{sc}}^{}$ denotes the isotropic elasticity tensor of the composite material obtained by the self-consistent method.

We can derive the equivalent elasticity tensor as follows [17]:

$C_{\text{sc}}^{}=C_{\text{m}}+\sum _{\alpha =1}^N{\text{f}}_{\alpha }\left(C_{\alpha }-C_m\right):{\left[I+{\overline{S}}_{\alpha }:{\left(C_{\text{sc}}^{}\right)}^{-1}:\left(C_{\alpha }^{}-C_{\text{sc}}^{}\right)\right]}^{-1}$(1.19)

The self-consistent method [17] gives an implicit expression of the equivalent modulus tensor, so that eq. (1.19) needs to be solved iteratively. The generalized self-consistent method assumes that there is a layer of matrix between the inclusion and the matrix, where the ratio of matrix to inclusion is the same as its ratio in the composite material. This method can lead to more accurate predictions.

The differential method [18] consists in adding the inclusion materials one by one to the matrix gradually, then regarding the resulting homogenized heterogeneous material as a new uniform matrix material, and calculating its elastic constants, and then repeating this process until the volume fraction of inclusions in the RVE reaches the volume fraction in the composite (see Figure 1.5).

Figure 1.5 Schematic diagram of the differential method

Taking a two-phase material as an example, assuming that the matrix volume is $\left|{\Omega }_{\text{m}}\right|$ and the inclusion volume is $\left|{\Omega }_{\text{c}}\right|$ , the equivalent elastic tensor of the composite can be expressed as:

$C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)=C_{\text{m}}+{\text{f}}_{\text{c}}\left(C_{\text{c}}-C_{\text{m}}\right):A\left(L_0\right)$(1.20)

where ${\text{f}}_{\text{c}}\text{=}\left|{\Omega }_{\text{c}}\right|\text{/}\left(\left|{\Omega }_{\text{m}}\right|\text{+}\left|{\Omega }_{\text{c}}\right|\right)$ denotes the volume fraction of the inclusion, and $A$ denotes the strain concentration tensor. After adding inclusions $\left|\Delta {\Omega }_{\text{c}}\right|$ in this “homogeneous” material, the equivalent elastic tensor of the newly obtained material is:

(1.21)$C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}+\Delta {\text{f}}_{\text{c}}\right)=C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)+\frac{\left|\Delta {\Omega }_c\right|}{\left|{\Omega }_m\right|+\left|{\Omega }_c\right|+\left|\Delta {\Omega }_c\right|}\left(C_{\text{c}}-C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)\right):A\left(C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)\right)$

Making $\Delta {\text{f}}_{\text{c}}\to 0$ , we can get a first-order differential equation to solve the equivalent elastic modulus, which can be written as:

$\frac{\text{d}{C}_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)}{{\text{df}}_{\text{c}}}=\frac{1}{1-{\text{f}}_{\text{c}}}\left(C_{\text{c}}-C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)\right):A\left(C_{\text{ds}}^{}\left({\text{f}}_{\text{c}}\right)\right)$(1.22)

The aforementioned analytical homogenization methods based on the Eshelby ellipsoidal inclusion theory are limited to materials with isotropic matrices and isotropic inclusions. In addition to the aforementioned commonly used methods, composite sphere models, progressive methods, and other methods are also sometimes used in composite material performance prediction.

1.2.2 Commonly Used Computer Modeling Methods of Microstructures of Materials

The analytical micromechanical methods for RVEs shown in Section 1.2.1 can obtain the overall macroscopic mechanical properties of composite materials quickly and effectively, but often do not consider the details of the geometry and random distribution of inclusions in the matrix, and cannot give the detailed stress and strain fields (such as the fields at the interfaces of the matrix and the inclusions/voids) inside the RVE. In order to address this problem, computational simulation methods are often used not only to directly compute the homogenized mechanical properties of heterogeneous materials but also to compute the detailed stress and strain fields at the interfaces, which are necessary for the study of damage and crack initiation in such materials.

The FEM is one of the most popular numerical simulation methods, and has been widely used to predict the equivalent elastic or elastic–plastic moduli of micromechanical RVEs, as well as the macroscopic mechanical response of composite structures. However, the FEM uses simple locally supported polynomials for the displacement trial functions, and it is often necessary to perform very detailed element meshing of the composite material when dealing with stress concentration problems (such as at the interfaces of inclusions, pores, and microcrack) [19] (see Figure 1.6(a)). This has led to the traditional FEM being usually used only to build RVEs with only a few inclusions, and ignoring the influence of the statistical distribution of heterogeneities in real situations. The simulation of real random microstructures with traditional finite elements requires complicated preprocessing efforts and a very high computational cost.

Figure 1.6

(a) Al/SiC unit cell model established by ABAQUS;

(b) meshing of the single SiCf/SiC fiber model by XFEM [23];

(c) meshing of BEM for solving the stress fields in porous materials [24]

In recent years, developments in the Extended Finite Element Method (XFEM) have effectively reduced the complexity of preprocessing. Unlike the traditional FEM, the core idea of XFEM is to introduce a jump function as the enrichment of the shape function to simulate the discontinuity of fields. Therefore, the description of the discontinuous field is completely independent of the mesh, and the mesh boundary does not need to coincide with the material interface, as shown in Figure 1.6(b). However, though XFEM can reduce the complexity of preprocessing, it cannot reduce the requirements of memory and computational time of traditional finite elements [20] when modeling a large number of heterogeneities.

In FEM, the elastic properties of the matrix and the inclusions can be in general anisotropic while evaluating the individual element stiffness matrices, and while assembling the stiffness matrix of the RVE. Compared with the FEM, the BEM only needs to discretize the mesh on the surface of the object, in general for only isotropic materials (as shown in Figure 1.6(c)), and the number of degrees of freedom of the discrete system is much lower than that of the FEM, and the solution accuracy of the stress field is higher, which is more suitable for solving the stress concentration problem of composite materials [21]. However, discretization by the BEM leads to unsymmetric and densely populated matrices; therefore, the computational cost of BEM will be significantly increased when the number of total degrees of freedom of the numerical model increases. Therefore, it is difficult to simulate complex microstructures with a lot of inclusions, and one needs to use fast solvers [22].

Moreover, to obtain the elasticity tensor $\overline{C}$ of a generally anisotropic composite material from the relationship of the stress tensor and the strain tensor, which is shown as:

$\left\{\begin{array}{c}〈{\sigma }_{11}〉\\ 〈{\sigma }_{22}〉\\ 〈{\sigma }_{33}〉\\ 〈{\sigma }_{12}〉\\ 〈{\sigma }_{23}〉\\ 〈{\sigma }_{31}〉\end{array}\right\}=\left[\begin{array}{cccccc}{\overline{C}}_{11}& {\overline{C}}_{12}& {\overline{C}}_{13}& {\overline{C}}_{14}& {\overline{C}}_{15}& {\overline{C}}_{16}\\ {\overline{C}}_{21}& {\overline{C}}_{22}& {\overline{C}}_{23}& {\overline{C}}_{24}& {\overline{C}}_{25}& {\overline{C}}_{26}\\ {\overline{C}}_{31}& {\overline{C}}_{32}& {\overline{C}}_{33}& {\overline{C}}_{34}& {\overline{C}}_{35}& {\overline{C}}_{36}\\ {\overline{C}}_{41}& {\overline{C}}_{42}& {\overline{C}}_{43}& {\overline{C}}_{44}& {\overline{C}}_{45}& {\overline{C}}_{46}\\ {\overline{C}}_{51}& {\overline{C}}_{52}& {\overline{C}}_{53}& {\overline{C}}_{54}& {\overline{C}}_{55}& {\overline{C}}_{56}\\ {\overline{C}}_{61}& {\overline{C}}_{62}& {\overline{C}}_{63}& {\overline{C}}_{64}& {\overline{C}}_{65}& {\overline{C}}_{66}\end{array}\right]\left\{\begin{array}{c}〈{\epsilon }_{11}〉\\ 〈{\epsilon }_{22}〉\\ 〈{\epsilon }_{33}〉\\ 2〈{\epsilon }_{12}〉\\ 2〈{\epsilon }_{23}〉\\ 2〈{\epsilon }_{31}〉\end{array}\right\}$(1.23)

we need to analyze at least six load cases. One can apply unit effective strains and calculate the resultant effective stresses, and the visa verse. For example, if we apply effective strains $〈{\epsilon }_{11}〉\text{=}1$ and $〈{\epsilon }_{22}〉\text{=}〈{\epsilon }_{33}〉\text{=}\text{2}〈{\epsilon }_{12}〉\text{=}\text{2}〈{\epsilon }_{23}〉\text{=}\text{2}〈{\epsilon }_{31}〉\text{=}0$ , then the computed effective stresses will be $〈{\sigma }_{11}〉\text{=}{\overline{C}}_{11},〈{\sigma }_{22}〉\text{=}{\overline{C}}_{12},〈{\sigma }_{33}〉\text{=}{\overline{C}}_{13},〈{\sigma }_{12}〉\text{=}{\overline{C}}_{14},〈{\sigma }_{23}〉\text{=}{\overline{C}}_{15},〈{\sigma }_{31}〉\text{=}{\overline{C}}_{16}.$ If we do this with each of the six strains, the corresponding six rows for the elasticity matrix $\overline{C}$ in eq. (1.23) can be obtained. However, in this way, one cannot guarantee that the computed $\overline{C}$ is symmetric. In another way, if we use the computed strain energy of the RVE to derive the elasticity tensor $\overline{C}$ , which will be shown in Chapter 2, we need to consider twenty-one load cases to calculate the twenty-one independent elasticity constants in general anisotropic materials. In this way, the symmetry of $\overline{C}$ is ensured. Either of these two strategies will further increase the burden of modeling microstructures of materials using FEM, XFEM, or BEM, making the direct modeling of a large number of inclusions, fibers, pores, or crack almost impossible. Such a dilemma has motivated the authors to develop Computational Grains (CGs), a highly efficient and accurate tool that facilitates the direct numerical simulation (DNS) of complex microstructures of heterogeneous materials.

1.2.3 Motivation for the Development of Computational Grains

To overcome the shortcomings and lack of general applicability of analytical micromechanics methods and computational methods such FEM, XFEM, and BEM, Dong and Atluri [25, 26] started to develop special 3D Trefftz Polyhedral Voronoi elements, each with an embedded spherical and ellipsoidal inclusion or void, for the modeling of composite and porous materials, which later became the first of a series of different CGs for different materials. As shown in Figure 1.7, such a CG is a polyhedral element with an embedded inclusion or pore, which can be seen as the generalization of Eshelby’s solution in a finite domain. It uses Voronoi diagram to create a set of elements, with each containing an embedded heterogeneity, which thus eliminates the workload of preprocessing to create a finite FEM mesh for complex microstructures. And because of the use of special trial functions and specifically designed boundary-type variational principles, the stiffness matrix of each of the CGs (with nodal displacements as degrees of freedom (DoFs)) can be directly computed and assembled. As can be seen from Figure 1.7(a–c), in this way, the unit cell modeled by tens of thousands of simple finite elements can be now accurately simulated by only one CG, and complex microstructures with many heterogeneities can be easily modeled by assembly of a number of CGs (see Figure 1.7(d)).

Figure 1.7

(a) Al/SiC unit cell model established by ABAQUS;

(b) schematic diagram of a typical computational grain (a polyhedron with a spherical inclusion);

(c) comparison of the stress by a CG and the analytical solution;

(d) RVE of 1,000 CGs with spherical inclusions

Such a work is later followed by the development of a series of different CGs for particulate composites, fiber composites, three-phase materials, nanocomposites, multifunctional composites, and so on. And the basic idea is that we use one grain/superelement to model the mechanical behavior of the smallest topological building block of various types of composite materials. Thus, each CG can either be a homogeneous material or include an embedded particle, fiber, void, or crack (as shown in Figure 1.8). However, different trial functions and variational principles or integral equations are used to develop each type of the CG, when different topologies, physics, constitutive relations, interfaces, and types of damages are to be modeled. Thus, “Computational Grains” is not a specific numerical algorithm, but a series of various computational models to achieve the same goal for different materials: DNS of complex microstructures of materials by the assembly of the most basic topological units modeled by specially designed elements named as Computational Grains. Various algorithms for different types of CGs and their applications in modeling heterogeneous materials will be presented in the following chapters of this book, in detail.

Figure 1.8 Schematic diagram of typical computational grains for homogeneous material/composites with spherical/ellipsoid/fiber/coated inclusions/microcrack

1.3 Purpose and Layout of This Book

This book is aimed at presenting high-performance computational tools for the modeling of RVEs for studying not only the homogenized properties of heterogeneous materials but also the stress/strain fields at the microlevel in detail by using CGs. It is intended for researchers in academia, industry, and government laboratories to understand the issues, challenges, and technologies involved in modeling heterogeneous materials as well as structures made of such materials in an integrated system. With advances in quantum computing, it will be possible to model entire structural components using the presently described micromechanical CGs directly rather than using the currently popular structural finite elements, in the near future. Thus damage precursors at the microlevel in a structural component can more easily be identified and modeled; and strength, degradation, and the life of a structural component can be more precisely predicted. The topics covered encompass elastic, nonelastic, and multi-physics problems for solid materials, particle-reinforced composites, fiber-reinforced composites, nano-size composites, porous materials, and so on. Several numerical examples are provided to validate and demonstrate the power of the developed methodologies. MATLAB codes for various CGs are also provided, so that the readers can implement them in their own computer programs or interface them with the off-the-shelf computer programs such as ABAQUS/ANSYS.

The rest of this book is organized as follows: Chapter 2 mainly introduces the basic concepts of the composite material homogenization method using micromechanics and the realization of attendant numerical algorithm. Chapter 3 introduces the basic principles of the computational grain method, and elaborates its entire process from the algorithm for describing the random distribution of inclusions; Voronoi meshing; finite element format construction; and the calculation of the equivalent moduli of the composite material; to the multi-scale simulation based on the computational grain method. Chapter 4 introduces the details of the derivation of Trefftz trial functions in 2D and 3D for various geometrical shapes of inclusions. Chapter 5 introduces the Trefftz CGs for particulate composites and porous materials in detail, with the corresponding boundary-type variational principles and its algorithmic implementation. Chapter 6 introduces a different boundary-type variational principle for fiber composites, and the algorithmic implementation of computational grain method for fiber composites. Chapter 7 introduces the formulation of CGs for nanocomposites considering the interface stress effect. Chapter 8 introduces a boundary-type variational principle for three-phase composites, and the algorithmic implementation of CGs with coated inclusions. Chapter 9 introduces the governing equations of viscoelasticity for the matrix and inclusions, and the formulation of CGs for viscoelastic composites. Chapter 10 introduces the CGs for piezoelectric composites/porous materials. Chapter 11 introduces the CGs with embedded microcracks in the matrix and inclusions. Chapter 12 introduces the multi-scale algorithm based on the computational grain method, which is used to quickly and accurately simulate the macro-stress state for structures and micro-stress states for materials.