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 [Reference Voigt8–Reference Hashin and Shtrikman12], (ii) homogenization methods based on the Eshelby inclusion theory [Reference Norris18], and (iii) differential methods [Reference Chawla, Ganesh and Wunsch19].

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

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 [Reference Fung and Tong9], 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 [Reference Voigt8], $C_{\text{v}}$ is an upper bound for the true equivalent elasticity tensor of composite materials. Reuss [Reference Reuss10] 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 [Reference Washizu11], $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 [Reference Hashin and Shtrikman12] 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 [Reference Eshelby and Peierls13] 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$ [Reference Nemat-Nasser, Hori and Datta14]:

${\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 [Reference Eshelby and Peierls13,Reference Eshelby and Peierls15] 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 [Reference Eshelby and Peierls15]:

$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 [Reference Mori and Tanaka16] 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 [Reference Norris18] 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 [Reference Hill17]:

$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 [Reference Hill17] 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 [Reference Norris18] 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) [Reference Chawla, Ganesh and Wunsch19] (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 [Reference Daggubati, Sharma and Pydi23];

(c) meshing of BEM for solving the stress fields in porous materials [Reference Song and Yin24]

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 [Reference Fries and Belytschko20] 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 [Reference Elmabrouk, Berger, Phan and Gray21]. 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 [Reference Liu, Nishimura and Otani22].

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 [Reference Dong and Atluri25, Reference Dong and Atluri26] 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