Micromechanics-Based Homogenization of the Effective Physical Properties of Composites With an Anisotropic Matrix and Interfacial Imperfections

Ryu, Seunghwa; Lee, Sangryun; Jung, Jiyoung; Lee, Jinyeop; Kim, Youngsoo

doi:10.3389/fmats.2019.00021

REVIEW article

Front. Mater., 01 March 2019
Sec. Mechanics of Materials
Volume 6 - 2019 | https://doi.org/10.3389/fmats.2019.00021

Micromechanics-Based Homogenization of the Effective Physical Properties of Composites With an Anisotropic Matrix and Interfacial Imperfections

Seunghwa Ryu¹^*^†

Sangryun Lee¹^†

Jiyoung Jung¹

Jinyeop Lee²

Youngsoo Kim¹

¹Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology (KAIST), Daejeon, South Korea
²Department of Mathematical Sciences, Korea Advanced Institute of Science and Technology (KAIST), Daejeon, South Korea

Micromechanics-based homogenization has been employed extensively to predict the effective properties of technologically important composites. In this review article, we address its application to various physical phenomena, including elasticity, thermal and electrical conduction, electric, and magnetic polarization, as well as multi-physics phenomena governed by coupled equations such as piezoelectricity and thermoelectricity. Especially, for this special issue, we introduce several research works published recently from our research group that consider the anisotropy of the matrix and interfacial imperfections in obtaining various effective physical properties. We begin with a brief review of the concept of the Eshelby tensor with regard to the elasticity and mean-field homogenization of the effective stiffness tensor of a composite with a perfect interface between the matrix and inclusions. We then discuss the extension of the theory in two aspects. First, we discuss the mathematical analogy among steady-state equations describing the aforementioned physical phenomena and explain how the Eshelby tensor can be used to obtain various effective properties. Afterwards, we describe how the anisotropy of the matrix and interfacial imperfections, which exist in actual composites, can be accounted for. In the last section, we provide a summary and outlook considering future challenges.

Introduction

Composites, typically referred to materials consisting of reinforcements and a matrix, offer many advantages that cannot be gained solely with only one of its constituents to improve various materials properties, such as resistance to chemicals (Jawaid et al., 2011; Taurino et al., 2016), a high strength-to-weight ratio (Walther et al., 2010), electrical (Allaoui et al., 2002; Gojny et al., 2005; Tuncer et al., 2007) or thermal insulation properties (Wei et al., 2011; Li et al., 2016), and/or combinations of these properties (Flahaut et al., 2000; Park et al., 2012). Most advanced materials in many engineering applications are composites, such as anti-corrosion composites in the marine industry (Mouritz et al., 2001), lightweight structural carbon composites in the car (Obradovic et al., 2012; Friedrich and Almajid, 2013) and airplane industries (Immarigeon et al., 1995; EL-Dessouky and Lawrence, 2013), and electric or thermal shielding composites (Imai et al., 2006; Zheng et al., 2009) as used in relation to electric wire and heat pipe technology.

In addition to manmade synthetic composites, nature has exploited the advantages of composites to meet certain requirements related to survival and sustained living conditions, even with the limited resources and building blocks available in nature (Gibson et al., 1995; Wegst and Ashby, 2004; Launey et al., 2009; Sen and Buehler, 2011). For example, the combined high toughness and high strength levels of nacre, bone, and cone shells are attributed to their unique hierarchical composite structures ranging from the nano to macro scale, and significant effort has been expended to understand and mimic the natural composites (Tang et al., 2003; Ji and Gao, 2004; Li et al., 2012; Hu et al., 2013; Das et al., 2015; Shao and Keten, 2015; Gao et al., 2017). For the facile design and application of various composites, it is of paramount importance to understand and predict the effective properties of composites as a function of their shape, volume fraction, and spatial distribution of the reinforcements used in them.

Numerical modeling based on finite element analyses has been widely used to predict the effective properties of composites, including the mechanical, thermal, electrical, piezoelectric, thermoelectric properties (Pan et al., 2008; Wang et al., 2011; Miled et al., 2013; Lu et al., 2014; Doghri et al., 2016; Lee et al., 2019, 2018a). However, to obtain statistically meaningful results, finite element analyses require multiple evaluations of large simulation cells involving a large number of fillers, thus necessitating a much fine mesh near the boundary to serve as representative volumes and leading to the requirement of computationally expensive and time-consuming calculations (Xu and Yagi, 2004; Marcos-Gomez et al., 2010; Lee et al., 2019, 2018a). Such extensive numerical calculations are feasible for predictions of the effective properties of systems in the linear response regime. However, it becomes a formidable task numerically to predict the behavior of composites in the nonlinear response regime, where multiple linearization and convergence tasks are required (Miled et al., 2013; Doghri et al., 2016).

When composites have relatively periodic and regular arrangements, they can be modeled reasonably well by an analytic model focusing on the load transfer mechanism of a unit cell. Typical examples include synthetic composites with very long fiber reinforcements (Suh, 2005) and natural staggered platelet structures (Gao, 2006; Kim et al., 2018). In contrast, composites involving random arrangements of short-fiber fillers or spherical reinforcements must be studied via homogenization schemes, as it is infeasible directly to model the interaction among randomly/aperiodically distributed reinforcements and load transfer mechanisms in detail.

The mechanical properties of such short-fiber-reinforced composites and multicomponent alloys can be deduced by computing the strain fields in the reinforcements and the matrix (Mura, 1982). Especially for composites including well-dispersed low-volume-fraction reinforcements (below 20%), the influence of other reinforcements can be taken into account via a mean-field approach by approximating the surrounding area of each reinforcement as the matrix subjected to local strain that is identical to the average strain of the matrix within the entire specimen (Benveniste, 1987) (Mori-Tanaka method) or using a medium with the effective stiffness of the composite (Hill, 1965) (self-consistent method). In this regime, where mean-field approximation is valid, a solution to the single-inhomogeneity problem can be applied to model the effective properties while taking into account multiple reinforcements. Once such homogenization scheme has been established for the linear response regime, and it can be extended to model the non-linear response of composites based on various linear approximation schemes such as the incremental tangent method (Castaneda and Tiberio, 2000; Idiart et al., 2006; Doghri et al., 2016), the incremental secant method (Wu et al., 2013), variational approaches (Castaneda, 1991; Lahellec and Suquet, 2007a; Brassart et al., 2011), and others (Nemat-Nasser, 1999; Lahellec and Suquet, 2007b).

Although extensive efforts have been devoted to determine various effective properties of composites, relatively less attention has been paid to the modeling of the effects of imperfect matrix-reinforcement interfaces and anisotropic matrices. Moreover, while the governing equations of various phenomena are mathematically analogous, the connections among them have not been discussed in detail in the literature. Hence, in this paper, we present research which considers interfacial imperfections and anisotropy effects for the predictions of various effective physical properties and suggest a universal formalism and numerical recipe based on a mathematical analogy. Especially, for this special issue, we introduce several papers published recently from our group on this issue.

The paper is structured as follows. In section Eshelby tensor and its application to effective stiffness calculations in elasticity, we provide a brief overview on the concept of the Eshelby tensor of isotropic materials in elasticity and its application to the computation of the effective stiffness of composites. In section Mathematical analogy between governing equations in various physical phenomena, we discuss the mathematical analogy among the steady-state equations governing various physical phenomena and how the Eshelby tensor and homogenization concept can be applied to predict various effective physical properties. Section Interfacial imperfections and anisotropy of the matrix describes how anisotropy of the matrix and interfacial imperfections, present in most realistic composites, can be accounted for. In the last section, we summarize the discussion and provide our perspective on future challenges.

Eshelby Tensor and Its Application to Effective Stiffness Calculations in Elasticity

Single Inclusion Problem

We can deduce the mechanical properties of fiber-reinforced composites by considering the strain field in inclusions and in inhomogeneities. An inclusion refers to an embedded material with an elastic stiffness tensor L_pqrs identical to that of the matrix, while an inhomogeneity refers to an embedded material with a different stiffness ${L^{'}}_{p q r s}$ . Eshelby showed that the strain field inside an ellipsoidal inclusion embedded in an infinite matrix is uniform when the inclusion is subject to uniform eigenstrain (Eshelby, 1957, 1959). Eigenstrain indicates the stress-free deformation strain associated with thermal expansion (Jun and Korsunsky, 2010), initial strain (Chiu, 1977), or phase transformation (George and Dvorak, 1992; Tirry and Schryvers, 2009). The Eshelby tensor is defined as the fourth-order tensor S_ijrs which links the constrained strain within the inclusion $ε_{i j}^{c}$ to the eigenstrain $ε_{r s}^{*}$ , as $ε_{i j}^{c} = S_{i j r s} ε_{r s}^{*}$ (Eshelby, 1957; Mura, 1982). In this section, we discuss (i) the key concept of the Eshelby tensor, (ii) how an ellipsoidal inhomogeneity subjected to an external load can be transformed into the equivalent Eshelby inclusion problem, and (iii) the mean-field method used to predict the effective stiffness of a composite based on the average theorem, considering an isotropic matrix in the absence of interfacial damage.

We begin the discussion by introducing the governing equation for Green's function G_ij(x − y) in elastostatics, which indicates the displacement in the i^th direction at point x by the unit body force in the j^th direction at point y, in an infinite elastic medium:

\begin{matrix} L_{i j k l} G_{k p, l j} (x - y) + δ_{i p} δ (x - y) = 0 & (1) \end{matrix}

Here, L_ijkl is the fourth-order elastic stiffness tensor, and the repeated indices represent the summation over all values from 1 to 3. For an isotropic material, L_ijkl can be represented by two independent elastic constants (Section Interfacial imperfections and anisotropy of the matrix includes a discussion of anisotropy). Green's function for the isotropic material is available in a closed form (Mura, 1982), as follows,

\begin{array}{l} G_{i j} (x - y) = \frac{1}{16 π μ (1 - ν) | x - y |} \\ \times [(3 - 4 ν) δ_{i j} + \frac{(x_{i} - y_{i}) (x_{j} - y_{j})}{{| x - y |}^{2}}] & (2) \end{array}

where μ and ν are the shear modulus and the Poisson's ratio of the material, respectively, and |x − y| denotes the standard norm of vector x − y.

The schematic for the single inclusion problem is depicted in Figure 1, which is solved with a four-step procedure. We assume that the inclusion can deform by the eigenstrain ε* when there is no external displacement or load (step 1). In order to maintain the original shape, the load T is applied to the inclusion (step 2). The inclusion is then plugged into a hole having the original shape and size within an infinite matrix (step 3). After removing the applied load (T), the inclusion exerts traction of F = −T on the matrix (step 4). Due to the constraining effect of the matrix, the inclusion deforms by the constrained strain ε^c, which is different from the eigenstrain ε*. The constrained strain field can be expressed by employing Green's function, as follows:

\begin{matrix} \begin{matrix} ε_{i j}^{c} = S_{i j r s} ε_{r s}^{*}, \\ where S_{i j r s} = \frac{1}{2} \int_{Ω}^{} L_{p q r s} {\frac{\partial^{2} G_{i p} (x - y)}{\partial x_{j} \partial y_{q}} + \frac{\partial^{2} G_{j p} (x - y)}{\partial x_{i} \partial y_{q}}} d y . & (3) \end{matrix} \end{matrix}

FIGURE 1

Figure 1. Schematic of the single inclusion problem in elasticity.

Because the eigenstrain and constrained strain are symmetric, the Eshelby tensor has minor symmetry (S_ijkl = S_jikl = S_ijlk) but no major symmetry (S_ijkl ≠ S_klij). In an isotropic medium, for an ellipsoidal inclusion with three different semi-axes a₁, a₂, and a₃, the expression is given in terms of standard elliptic integrals (Mura, 1982; Qu and Cherkaoui, 2007). For isotropic material having an axi-symmetric ellipsoidal inclusion (a₁ ≠ a₂ = a₃), all components of the Eshelby tensor can be obtained in a closed form, as shown in Appendix A (Mura, 1982; Lee and Ryu, 2018). The orientation average of a randomly distributed ellipsoidal inclusion was considered in one of our studies (Lee and Ryu, 2018). For a spherical inclusion, the Eshelby tensor can be compactly expressed as

\begin{matrix} S_{i j k l} = \frac{5 ν - 1}{15 (1 - ν)} δ_{i j} δ_{k l} + \frac{4 - 5 ν}{15 (1 - ν)} (δ_{i k} δ_{j l} + δ_{i l} δ_{j k}) . & (4) \end{matrix}

We note that the Eshelby tensor for the spherical inclusion does not depend on the radius. In general, the Eshelby tensor does not depend on the absolute size of the inclusion because there is no characteristic length scale in elasticity. In contrast, in the presence of an interfacial imperfection, the modified Eshelby tensor depends on the absolute size of the inclusion, which will be discussed in Section Interfacial imperfections and anisotropy of the matrix. We will also describe the Eshelby tensor for an anisotropic medium in Section Interfacial imperfections and anisotropy of the matrix.

Single Inhomogeneity Problem

Having introduced the concept of the Eshelby tensor, we discuss how an ellipsoidal inhomogeneity subjected to a uniform external load (external strain $ε_{k l}^{A}$ for the following example) can be transformed into the equivalent Eshelby inclusion problem. In the presence of a uniform external load (either stress or strain), the total strain field of the aforementioned Eshelby inclusion problem becomes $ε_{k l}^{t} = ε_{k l}^{A} + ε_{k l}^{c}$ , and the stress field becomes $σ_{i j}^{I} = L_{i j k l} (ε_{k l}^{t} - ε_{k l}^{*, Eq})$ . When the total strain field $ε_{k l}^{t'} = ε_{k l}^{A} + ε_{k l}^{c'}$ and stress field $σ_{i j}^{I'} = L_{i j k l}^{'} ϵ_{k l}^{t'}$ of an inhomogeneity subjected to an identical applied load are identical to those of the Eshelby inclusion, $ε_{k l}^{*, Eq}$ is referred to as an equivalent eigenstrain. Therefore, from the total strain and stress equality conditions, or equivalently $ε_{k l}^{A} + ε_{k l}^{c'} = ε_{k l}^{A} + ε_{k l}^{c}$ and $L_{i j k l}^{'} (ε_{k l}^{A} + ε_{k l}^{c'}) = L_{i j k l} (ε_{k l}^{A} + ε_{k l}^{c} - ε_{k l}^{*, Eq})$ , we obtain the following expression for the equivalent eigenstrain:

\begin{matrix} [(L_{i j k l}^{'} - L_{i j k l}) S_{k l m n} + L_{i j m n}] ε_{m n}^{*, Eq} = (L_{i j k l}^{'} - L_{i j k l}) ε_{k l}^{A} . & (5) \end{matrix}

In a tensor expression,

\begin{matrix} ε^{*, Eq} = {[(L^{'} - L) : S + L]}^{- 1} : (L - L^{'}) : ε^{A} & (6) \end{matrix}

where : and []⁻¹ refer to the double contraction and inverse operator, respectively. The double contraction and inverse operation of fourth-order tensors can be facilitated by adopting the Mandel notation (Helnwein, 2001), which converts all tensor operations into conventional 6 × 6 matrix operations (Lee et al., 2019, 2018b).

In the Mandel notation, the strain vector $\vec{σ}$ , stress vector $\vec{ε}$ , stiffness matrix 〈L〉, and Eshelby matrix 〈S〉 are defined as follows:

\begin{matrix} \vec{σ} : = [\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \\ \sqrt{2} σ_{23} \\ \sqrt{2} σ_{31} \\ \sqrt{2} σ_{12} \end{matrix}], \vec{ε} : = [\begin{matrix} ε_{11} \\ ε_{22} \\ ε_{33} \\ \sqrt{2} ε_{23} \\ \sqrt{2} ε_{31} \\ \sqrt{2} ε_{12} \end{matrix}], & (7) \end{matrix}

\begin{matrix} \begin{array}{l} 〈 L 〉 : = [\begin{matrix} L_{1111} L_{1122} L_{1133} & \sqrt{2} L_{1123} \sqrt{2} L_{1131} \sqrt{2} L_{1112} \\ L_{1122} L_{2222} L_{2233} & \sqrt{2} L_{2223} \sqrt{2} L_{2231} \sqrt{2} L_{2212} \\ L_{1133} L_{2233} L_{3333} & \sqrt{2} L_{3323} \sqrt{2} L_{3331} \sqrt{2} L_{3312} \\ \sqrt{2} L_{1123} \sqrt{2} L_{2223} \sqrt{2} L_{3323} & 2 L_{2323} 2 L_{2331} 2 L_{2312} \\ \sqrt{2} L_{1131} \sqrt{2} L_{2231} \sqrt{2} L_{3331} & 2 L_{2331} 2 L_{3131} 2 L_{3112} \\ \sqrt{2} L_{1112} \sqrt{2} L_{2212} \sqrt{2} L_{3312} & 2 L_{2312} 2 L_{3112} 2 L_{1212} \end{matrix}] \\ 〈 s 〉 : = [\begin{matrix} S_{1111} S_{1122} S_{1133} & \sqrt{2} S_{1123} \sqrt{2} S_{1131} \sqrt{2} S_{1112} \\ S_{2211} S_{2222} S_{2233} & \sqrt{2} S_{2223} \sqrt{2} S_{2231} \sqrt{2} S_{2212} \\ S_{3311} S_{3322} S_{3333} & \sqrt{2} S_{3323} \sqrt{2} S_{3331} \sqrt{2} S_{3312} \\ \sqrt{2} S_{2311} \sqrt{2} S_{2322} \sqrt{2} S_{2333} & 2 S_{2323} 2 S_{2331} 2 S_{2312} \\ \sqrt{2} S_{3111} \sqrt{2} S_{3122} \sqrt{2} S_{3133} & 2 S_{3123} 2 S_{3131} 2 S_{3112} \\ \sqrt{2} S_{1211} \sqrt{2} S_{1222} \sqrt{2} S_{1233} & 2 S_{1223} 2 S_{1231} 2 S_{1212} \end{matrix}] \end{array} \end{matrix}

The prefactors $\sqrt{2}$ and 2 ensure that the matrix-matrix product and the inverse coincide with the double contraction and the inverse of the fourth-order tensors, respectively. Hence, if A, B denote a fourth-order tensor with minor symmetry (A_ijkl = A_jikl = A_ijlk), and 〈A〉, 〈B〉 represent the corresponding 6 × 6 matrix following the Mandel notation, we can calculate the double contraction and inverse from the 6 × 6 matrix multiplication and inverse, respectively, as

\begin{matrix} 〈 A : B 〉 = 〈 A 〉 〈 B 〉, 〈 A^{- 1} 〉 = {〈 A 〉}^{- 1} . & (8) \end{matrix}

Mean-Field Homogenization of the Effective Stiffness

For the problem involving multiple inhomogeneities, the Mori-Tanaka method has been widely used for predicting the effective stiffness by considering the interaction among the inhomogeneities (see Figure 2). Each inhomogeneity is assumed to be embedded effectively in the matrix subjected to a spatial average strain of the matrix over the entire specimen, which is why the technique is referred to as the mean-field homogenisation method. Because the interaction between the inhomogeneities becomes intense at a high volume fraction of inhomogeneities, the Mori–Tanaka method is known to be reliable at a relatively low volume fraction of inhomogeneities (<20%). In this mean-field homogenization scheme, the average strain fields in the matrix, inhomogeneity, and composite ( $\bar{ε_{0}}, \bar{ε_{1}} (\approx \bar{ε_{0}} + S : ε^{*, Eq}),$ and $\bar{ε} (= ε^{A})$ , respectively) are related as follows:

\begin{matrix} \bar{ε_{1}} = T \bar{ε_{0}} = T : {(c_{0} I + c_{1} T)}^{- 1} \bar{ε} = A : \bar{ε}, & (9) \end{matrix}

where $A = {T : (c_{0} I + c_{1} T)}^{- 1}$ (where $T = {[I + S : L_{0}^{- 1} : (L_{1} - L_{0})]}^{- 1}$ ) is the strain concentration tensor for the Mori-Tanaka method. A converges to the strain concentration tensor for the single inhomogeneity problem, T, in the limit of very small inhomogeneity volume fraction, c₁. The effective stiffness can be obtained from the relationship between the volume-averaged stress and strain within the composite ( $\bar{σ} = c_{0} \bar{σ_{0}} + c_{1} \bar{σ_{1}} = c_{0} L_{0} : \bar{ε_{0}} + c_{1} L_{1} : \bar{ε_{1}}$ and $\bar{ε} = c_{0} \bar{ε_{0}} + c_{1} \bar{ε_{1}}$ , respectively) as follows:

\begin{matrix} \bar{σ} = L_{e f f} : \bar{ε} L_{e f f} = (c_{0} L_{0} + c_{1} L_{1} : T) : {(c_{0} I + c_{1} T)}^{- 1} & (10) \end{matrix}

FIGURE 2

Figure 2. Non-uniform σ_xx field under uniform applied strain in the x direction within a cross-section in a three-dimensional volume simulation. The material properties used for this figure are E₁ = 10E₀and ν₁ = ν₀ = 0.25 and the volume fraction is 10%.

There is another class of homogenization method called the self-consistent method, where the surrounding area is approximated as a medium with the effective stiffness of the composite. However, we omit the discussion of the self-consistent method here because the effective stiffness cannot be solved explicitly (Hill, 1965; Qu and Cherkaoui, 2007). Instead, it requires an iterative numerical solution.

Mathematical Analogy Between Governing Equations in Various Physical Phenomena

Mathematical Analogy Discussion

Having introduced the key concepts of the Eshelby tensor and the mean-field homogenization scheme, we turn our attention to the predictions of various physical phenomena via a mathematical analogy. We consider two representative classes of steady-state equations governing different physics: conduction/dielectric-related equations and coupled multi-physics equations. An extensive discussion can be found in the literature (Milton, 2002).

For the first class of the equations, referred to as conductivity/dielectric-related equations, the governing equation under a steady state at any point x in a medium in the absence of an internal source can be written using the following equations,

\begin{matrix} j (x) = K (x) \cdot e (x), \nabla \cdot j = 0, \nabla \times e = 0, & (11) \end{matrix}

which appear in physical problems concerning the conductivity of dielectrics, such as those linked to the concepts of electrical conductivity, dielectric phenomena, magnetism, thermal conduction, diffusion, and flows in porous media. In each field, the vector fields j(x), e(x) and, the second-order tensor K(x) have the physical interpretations given in Table 1. We note that the governing equations for these phenomena in a dynamic regime differ from each other; equations pertaining to dissipative transport phenomena use a first-order derivative with respect to time, leading to a homogeneous distribution of the temperature or concentration in an infinite time limit if an external driving force is not present, whereas equations pertaining to non-dissipative transport phenomena use a second-order derivative with respect to time, leading to electro-magnetic or acoustic wave propagation (Kreyszig, 1972). Because the steady-state equations are identical in terms of their mathematical descriptions, the homogenization scheme obtained for one of these equations can be applied directly to any of the other equations in Table 1 to obtain the effective electrical and thermal conductivities, the effective dielectric and magnetic permittivities, the effective diffusivity, and the effective fluid permeability.

TABLE 1

Table 1. Interpretation of conductivity/dielectric-related equations under a steady state, as governed by a single class of PDE.

As examples of the second class of equations, we consider piezoelectric and thermoelectric equations. The aforementioned conductivity/dielectric-related equations describe idealized cases. In practice, lattice distortion can cause electric polarization in some classes of anisotropic materials, and vice versa (Fu and Cohen, 2000). Moreover, conducting electrons can carry some heat with them and thus causes coupling between thermal and electric phenomena (Zhao et al., 2014). The former is referred to as piezoelectricity and the latter is termed thermoelectricity.

The constitutive equations of piezoelectricity are given as

\begin{matrix} \begin{matrix} σ_{i j} = L_{i j m n} ε_{m n} + d_{n i j} (- E_{n}) \\ D_{i} = d_{i m n} ε_{m n} - ϵ_{i n} (- E_{n}), \end{matrix} & (12) \end{matrix}

where D_i is the electrical displacement vector, E_n is the electric field vector, $ϵ_{m k} = {(\frac{\partial D_{m}}{\partial E_{k}})}_{ε}$ is the second rank tensor of the dielectric permittivity constants under constant strain with ϵ_mk = ϵ_km, and ${- d}_{i j m} = {(\frac{\partial σ_{i j}}{\partial E_{m}})}_{ε} = {- (\frac{\partial D_{m}}{\partial ϵ_{i j}})}_{E}$ is the third rank tensor of the piezoelectric constants under constant stress σ_ij with d_kij = d_kji. We can express Equation (12) as a simplified form of a linear equation, similar to the constitutive equations of elasticity, thermal conduction and electrical conduction using the notation introduced by Barnett and Lothe (1975). This notation is identical to the conventional notation, except that the repeated capital index implies summation from 1 to 4. With this notation, we can express the elastic strain and electric field as

\begin{matrix} Z_{M n} = {\begin{matrix} ε_{m n} & M = 1, 2, 3 \\ - E_{n} & M = 4 \end{matrix}, & (13) \end{matrix}

where Z_Mn is obtained by the differentiate U_M, which is given by

\begin{matrix} U_{M} = {\begin{matrix} u_{m} & M = 1, 2, 3 \\ V & M = 4 \end{matrix} . & (14) \end{matrix}

In a similar manner, we can simplify the stress and electric displacement to one matrix using the notation

\begin{matrix} Σ_{i J} = {\begin{matrix} σ_{i j} & J = 1, 2, 3 \\ D_{i} & J = 4 \end{matrix} . & (15) \end{matrix}

Then, the electroelastic moduli and constitutive equation are obtained using the equations below.

\begin{matrix} C_{i J M n} = {\begin{matrix} L_{i j m n} & J, M = 1, 2, 3 \\ d_{n i j} & J = 1, 2, 3, M = 4 \\ d_{i m n} & J = 4, M = 1, 2, 3 \\ - ϵ_{i n} & J, M = 4 \end{matrix} & (16) \end{matrix}

\begin{matrix} Σ_{i J} = C_{i J M n} Z_{M n} \end{matrix}

On the other hand, for a thermoelectric material, similarly, the constitutive equation is arranged in tensor form, as follows:

\begin{matrix} [\frac{J_{E}}{J_{S}}] = [\begin{matrix} σ & α \cdot σ \\ α \cdot σ & γ / T \end{matrix}] [\begin{matrix} E \\ e \end{matrix}] & (17) \end{matrix}

By using the notation introduced for a piezoelectric material (Jung et al., 2018), the expression is as follows,

\begin{matrix} J_{J i} = P_{J i M n} Q_{M n} & (18) \end{matrix}

where

\begin{matrix} J_{J i} = {\begin{matrix} J_{i}^{E} & J = 1 \\ J_{i}^{S} & J = 2 \end{matrix}}, & (19) \end{matrix}

\begin{matrix} Q_{M n} = {\begin{matrix} E_{n} & M = 1 \\ e_{n} & M = 2 \end{matrix}}, & (20) \end{matrix}

\begin{matrix} P_{J i M n} = {\begin{array}{l} {(σ_{e})}_{i n} & J = 1, & M = 1 \\ {(σ_{e} \cdot α)}_{i n} & J = 1, & M = 2 \\ {(α \cdot σ_{e})}_{i n} & J = 2, & M = 1 \\ γ_{i n} / T & J = 2, & M = 2 \end{array}} . & (21) \end{matrix}

Here, J^E is the electrical current density, J^S(= q/T) is the entropy flux, q is the heat flux, T is the temperature, σ_e is the electrical conductivity, α is the Seebeck coefficient, γ(= κ+Tα·σ_e·α) is the heat conductivity under a zero electric field, and κ is the heat conductivity at zero current. As shown in Equation (21), P_JiMn contains T, which depends on the position x and therefore does not have mathematical similarity with other physics in general. For a small temperature difference across the hot and cool sides (T_H − T_L = ΔT ≈ 0), we can assume that the temperature variable T in the constitutive equation is a constant average temperature $T_{a v g} (= \frac{T_{H} + T_{L}}{2})$ , allowing us to obtain a linear constitutive equation. The validity of this assumption has been studied for wide range of ΔT by comparing theoretical predictions with finite element analysis outcomes (Jung et al., 2018).

Eshelby Tensor for Various Physical Phenomena

Due to the mathematical analogy, the Eshelby tensor can be expressed in a similar form. In the most general case of a medium with arbitrary anisotropy and an ellipsoidal inclusion with three different semi-axes a₁, a₂, and a₃, by simplifying Equation (3), the elastic Eshelby tensor can be expressed as Mura (1982),

\begin{matrix} S_{i j r s} = \frac{1}{8 π} L_{p q r s} \int_{- 1}^{1} \int_{0}^{2 π} {g_{i p j q} (\bar{ξ}) + g_{j p i q} (\bar{ξ})} d θ d \bar{ζ_{3}}, & (22) \end{matrix}

where $g_{i j k l} (\bar{ξ}) = {\bar{ξ}}_{k} {\bar{ξ}}_{l} Z_{i j} (\bar{ξ})$ . Here, Z(ξ)=[(L · ξ) ·ξ]⁻¹ and ξ are Green's function and a vector in the Fourier space, respectively. Because Z_ij(ξ) is a homogeneous function of degree −2, ξ_kξ_lZ_ij(ξ) is identical to ${\bar{ξ}}_{k} {\bar{ξ}}_{l} Z_{i j} (\bar{ξ}),$ where $\bar{ξ}$ is a normalized vector. Equation (22) can be evaluated by using the integral variables $\bar{ζ_{3}}$ and θ; thus, ${\bar{ξ}}_{1} = \frac{1}{a_{1}} \sqrt{1 - {\bar{ζ_{3}}}^{2}} cos θ$ , ${\bar{ξ}}_{2} = \frac{1}{a_{2}} \sqrt{1 - {\bar{ζ_{3}}}^{2}} sin θ$ , $a n d {\bar{ξ}}_{3} = \frac{1}{a_{3}} \bar{ζ_{3}}$ given the semi-axes of an ellipsoidal inclusion (a_i). Similarly, by adding three more degrees of freedom from electric polarization and its coupling with lattice distortion, the pizeoelectric Eshelby tensor is derived as follows (Dunn and Taya, 1993),

\begin{matrix} S_{M n A b} = {\begin{matrix} \frac{1}{8 π} C_{i J A b} \int_{- 1}^{1} \int_{0}^{2 π} [h_{M J i n} + h_{n J i M}] d θ d \bar{ζ_{3}}, & M = 1, 2, 3 \\ \frac{1}{4 π} C_{i J A b} \int_{- 1}^{1} \int_{0}^{2 π} h_{4 J i n} d θ d \bar{ζ_{3}} . & M = 4 \end{matrix} & (23) \end{matrix}

Here, $h_{M J k l} (\bar{ξ}) = {\bar{ξ}}_{k} {\bar{ξ}}_{l} {(Y^{- 1})}_{M J}$ where $Y_{M J} (\bar{ξ}) = {\bar{ξ}}_{s} {\bar{ξ}}_{t} C_{s M J t}$ . Because the piezoelectric coefficients are zero for all centrosymmetric crystals, the piezoelectric Eshelby tensor for an isotropic medium does not exist.

For the conduction or dielectric phenomena described by Equation (11) and physical property tensor K in Table 1, the Eshelby tensor is obtained from

\begin{matrix} S_{i k} (x) = \frac{\partial}{\partial x_{i}} \int_{Ω}^{} \frac{\partial G (x - y)}{\partial y_{j}} d y K_{j k} . & (24) \end{matrix}

Unlike the Green's function in elasticity, the Green's function in the conduction or dielectric phenomena for a medium with arbitrary symmetry has been derived in a closed form as follows:

\begin{matrix} G (x - y) = \frac{1}{4 π \sqrt{\det (K) [{(x - y)}^{T} K^{T} (x - y)]}} . & (25) \end{matrix}

The second-order Eshelby tensor can be further simplified for an arbitrary ellipsoidal inclusion as an integral having one variable (Giordano and Palla, 2008; Lee et al., 2018a):

\begin{matrix} S = \frac{\det (a)}{2} \int_{0}^{\infty} \frac{{(a^{2} + s K)}^{- 1} \cdot K}{\sqrt{\det (a^{2} + s K)}} d s where a = [\begin{matrix} a_{1} & 0 & 0 \\ 0 & a_{2} & 0 \\ 0 & 0 & a_{3} \end{matrix}], & (26) \end{matrix}

Or it can be written in an integral involving two variables:

\begin{matrix} S_{i k} = \frac{1}{4 π} K_{i j} \int_{- 1}^{1} \int_{0}^{2 π} {\bar{ξ}}_{j} {\bar{ξ}}_{k} H^{- 1} d θ d \bar{ζ_{3}} & (27) \end{matrix}

where $H = {\bar{ξ}}_{p} {\bar{ξ}}_{q} K_{p q}$ . For a spherical inclusion in an isotropic matrix, the Eshelby tensor can be simplified as $S_{i k} = \frac{1}{3} δ_{i k}$ . For a thermoelectric material, with the coupling of electrical and heat conduction, the Eshelby tensor can be written as follows Jung et al. (2018),

\begin{matrix} S_{M n A b} = \frac{1}{4 π} P_{P i A b} \int_{- 1}^{1} \int_{0}^{2 π} {\bar{ξ}}_{i} {\bar{ξ}}_{n} {(W^{- 1})}_{M P} d θ d \bar{ζ_{3}} & (28) \end{matrix}

where $W_{M P} = {\bar{ξ}}_{a} {\bar{ξ}}_{t} P_{M s P t}$ . The Eshelby tensor for a spherical inclusion in an isotropic matrix becomes $S_{M n A b} = \frac{1}{3} δ_{M A} δ_{n b}$ . We note that the Eshelby tensor expressions for the conduction/dielectric and thermoelectric phenomena are simpler than those for elasticity and piezoelectricity, because the Green's functions for the latter consider a vector field (displacement) whereas the Green's functions for the former only involve scalar fields (such as temperature and electrical potential).

Examples of Numerical Calculations and FEA Validation

We now turn our attention to predict the effective properties of composites based on a mathematically analogous formula. In Mandel notation represented by p = Xq, where X is a N × N matrix and the input (output) field q (p) has N components, the linear operator in elasticity L_ijmn, conduction (or dielectric) κ_in, piezoelectricity C_iJMn, and thermoelectricity P_JiMn become 6 × 6, 3 × 3, 9 × 9, and 6 × 6 symmetric matrices, respectively. Although most existing studies (Dunn and Taya, 1993; Huang and Kuo, 1996; Odegard, 2004; Duschlbauer et al., 2006; Martinez-Ayuso et al., 2017) utilize the Voigt notation, we adapt the Mandel notation here. For example, if we use the Voigt notation for a piezoelectric case, the linear operator matrix (which corresponds to the material properties) and the Eshelby matrix can be expressed as Equations (29, 30), respectively:

\begin{matrix} [\begin{matrix} \begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \end{matrix} \\ \begin{matrix} σ_{23} \\ σ_{31} \\ σ_{12} \end{matrix} \\ \begin{matrix} D_{1} \\ D_{2} \\ D_{3} \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} C_{1111} & C_{1122} & C_{1133} \\ C_{2211} & C_{2222} & C_{2233} \\ C_{3311} & C_{3322} & C_{3333} \end{matrix} & \begin{matrix} C_{1123} & C_{1131} & C_{1112} \\ C_{2223} & C_{2231} & C_{2212} \\ C_{3323} & C_{3331} & C_{3312} \end{matrix} & \begin{matrix} C_{1141} & C_{1142} & C_{1143} \\ C_{2241} & C_{2242} & C_{2243} \\ C_{3341} & C_{3342} & C_{3343} \end{matrix} \\ \begin{matrix} C_{2311} & C_{2322} & C_{2333} \\ C_{3111} & C_{3122} & C_{3133} \\ C_{1211} & C_{1222} & C_{1233} \end{matrix} & \begin{matrix} C_{2323} & C_{2331} & C_{2312} \\ C_{3123} & C_{3131} & C_{3112} \\ C_{1223} & C_{1231} & C_{1212} \end{matrix} & \begin{matrix} C_{2341} & C_{2342} & C_{2343} \\ C_{3141} & C_{3142} & C_{3143} \\ C_{1241} & C_{1242} & C_{1243} \end{matrix} \\ \begin{matrix} C_{1411} & C_{1422} & C_{1433} \\ C_{2411} & C_{2422} & C_{2433} \\ C_{3411} & C_{3422} & C_{3433} \end{matrix} & \begin{matrix} C_{1423} & C_{1431} & C_{1412} \\ C_{2423} & C_{2431} & C_{2412} \\ C_{3423} & C_{3431} & C_{3412} \end{matrix} & \begin{matrix} C_{1441} & C_{1442} & C_{1443} \\ C_{2441} & C_{2442} & C_{2443} \\ C_{3441} & C_{3442} & C_{3443} \end{matrix} \end{matrix}] [\begin{matrix} \begin{matrix} ε_{11} \\ ε_{22} \\ ε_{33} \end{matrix} \\ \begin{matrix} 2 ε_{23} \\ 2 ε_{31} \\ 2 ε_{12} \end{matrix} \\ \begin{matrix} - E_{1} \\ - E_{2} \\ - E_{3} \end{matrix} \end{matrix}] & (29) \end{matrix}

\begin{matrix} [\begin{matrix} \begin{matrix} ε_{11} \\ ε_{22} \\ ε_{33} \end{matrix} \\ \begin{matrix} 2 ε_{23} \\ 2 ε_{31} \\ 2 ε_{12} \end{matrix} \\ \begin{matrix} - E_{1} \\ - E_{2} \\ - E_{3} \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} S_{1111} & S_{1122} & S_{1133} \\ S_{2211} & S_{2222} & S_{2233} \\ S_{3311} & S_{3322} & S_{3333} \end{matrix} & \begin{matrix} S_{1123} & S_{1131} & S_{1112} \\ S_{2223} & S_{2231} & S_{2212} \\ S_{3323} & S_{3331} & S_{3312} \end{matrix} & \begin{matrix} S_{1141} & S_{1142} & S_{1143} \\ S_{2241} & S_{2242} & S_{2243} \\ S_{3341} & S_{3342} & S_{3343} \end{matrix} \\ \begin{matrix} 2 S_{2311} & 2 S_{2322} & 2 S_{2333} \\ 2 S_{3111} & 2 S_{3122} & 2 S_{3133} \\ 2 S_{1211} & 2 S_{1222} & 2 S_{1233} \end{matrix} & \begin{matrix} 2 S_{2323} & 2 S_{2331} & 2 S_{2312} \\ 2 S_{3123} & 2 S_{3131} & 2 S_{3112} \\ 2 S_{1223} & 2 S_{1231} & 2 S_{1212} \end{matrix} & \begin{matrix} 2 S_{2341} & 2 S_{2342} & 2 S_{2343} \\ 2 S_{3141} & 2 S_{3142} & 2 S_{3143} \\ 2 S_{1241} & 2 S_{1242} & 2 S_{1243} \end{matrix} \\ \begin{matrix} S_{4111} & S_{4122} & S_{4133} \\ S_{4211} & S_{4222} & S_{4233} \\ S_{4311} & S_{4322} & S_{4333} \end{matrix} & \begin{matrix} S_{4123} & S_{4131} & S_{4112} \\ S_{4223} & S_{4231} & S_{4212} \\ S_{4323} & S_{4331} & S_{4312} \end{matrix} & \begin{matrix} S_{4141} & S_{4142} & S_{4143} \\ S_{4241} & S_{4242} & S_{4243} \\ S_{4341} & S_{4342} & S_{4343} \end{matrix} \end{matrix}] [\begin{matrix} \begin{matrix} ε_{11}^{*} \\ ε_{22}^{*} \\ ε_{33}^{*} \end{matrix} \\ \begin{matrix} 2 ε_{23}^{*} \\ 2 ε_{31}^{*} \\ 2 ε_{12}^{*} \end{matrix} \\ \begin{matrix} - E_{1}^{*} \\ - E_{2}^{*} \\ - E_{3}^{*} \end{matrix} \end{matrix}] & (30) \end{matrix}

As shown by Equations (29, 30), the coefficients in the transformed matrix are different for the two set of matrices. In contrast, with the Mandel notation, the two matrices can be expressed as

\begin{matrix} [\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \\ \sqrt{2} σ_{23} \\ \sqrt{2} σ_{31} \\ \sqrt{2} σ_{12} \\ D_{1} \\ D_{2} \\ D_{3} \end{matrix}] = [\begin{matrix} C_{1111} C_{1122} C_{1133} & \sqrt{2} C_{1123} \sqrt{2} C_{1131} \sqrt{2} C_{1112} & C_{1141} C_{1142} C_{1143} \\ C_{2211} C_{2222} C_{2233} & \sqrt{2} C_{2223} \sqrt{2} C_{2231} \sqrt{2} C_{2212} & C_{2241} C_{2242} C_{2243} \\ C_{3311} C_{3322} C_{3333} & \sqrt{2} C_{3323} \sqrt{2} C_{3331} \sqrt{2} C_{3312} & C_{3341} C_{3342} C_{3343} \\ \sqrt{2} C_{2311} \sqrt{2} C_{2322} \sqrt{2} C_{2333} & 2 C_{2323} 2 C_{2331} 2 C_{2312} & \sqrt{2} C_{2341} \sqrt{2} C_{2342} \sqrt{2} C_{2343} \\ \sqrt{2} C_{3111} \sqrt{2} C_{3122} \sqrt{2} C_{3133} & 2 C_{3123} 2 C_{3131} 2 C_{3112} & \sqrt{2} C_{3141} \sqrt{2} C_{3142} \sqrt{2} C_{3143} \\ \sqrt{2} C_{1211} \sqrt{2} C_{1222} \sqrt{2} C_{1233} & 2 C_{1223} 2 C_{1231} 2 C_{1212} & \sqrt{2} C_{1241} \sqrt{2} C_{1242} \sqrt{2} C_{1243} \\ C_{1411} C_{1422} C_{1433} & \sqrt{2} C_{1423} \sqrt{2} C_{1431} \sqrt{2} C_{1412} & C_{1441} C_{1442} C_{1443} \\ C_{2411} C_{2422} C_{2433} & \sqrt{2} C_{2423} \sqrt{2} C_{2431} \sqrt{2} C_{2412} & C_{2441} C_{2442} C_{2443} \\ C_{3411} C_{3422} C_{3433} & \sqrt{2} C_{3423} \sqrt{2} C_{3431} \sqrt{2} C_{3412} & C_{3441} C_{3442} C_{3443} \end{matrix}] [\begin{array}{l} ε_{11} \\ ε_{22} \\ ε_{33} \\ \sqrt{2} ε_{23} \\ \sqrt{2} ε_{31} \\ \sqrt{2} ε_{12} \\ E_{1} \\ - E_{2} \\ - E_{3} \end{array}] & (31) \end{matrix}

and

\begin{matrix} [\begin{matrix} \begin{matrix} ε_{11} \\ ε_{22} \\ ε_{33} \end{matrix} \\ \begin{matrix} \sqrt{2} ε_{23} \\ \sqrt{2} ε_{31} \\ \sqrt{2} ε_{12} \end{matrix} \\ \begin{matrix} - E_{1} \\ - E_{2} \\ - E_{3} \end{matrix} \end{matrix}] = [\begin{matrix} \begin{matrix} S_{1111} & S_{1122} & S_{1133} \\ S_{2211} & S_{2222} & S_{2233} \\ S_{3311} & S_{3322} & S_{3333} \end{matrix} & \begin{matrix} \sqrt{2} S_{1123} & \sqrt{2} S_{1131} & \sqrt{2} S_{1112} \\ \sqrt{2} S_{2223} & \sqrt{2} S_{2231} & \sqrt{2} S_{2212} \\ \sqrt{2} S_{3323} & \sqrt{2} S_{3331} & \sqrt{2} S_{3312} \end{matrix} & \begin{matrix} S_{1141} & S_{1142} & S_{1143} \\ S_{2241} & S_{2242} & S_{2243} \\ S_{3341} & S_{3342} & S_{3343} \end{matrix} \\ \begin{matrix} \sqrt{2} S_{2311} & \sqrt{2} S_{2322} & \sqrt{2} S_{2333} \\ \sqrt{2} S_{3111} & \sqrt{2} S_{3122} & \sqrt{2} S_{3133} \\ \sqrt{2} S_{1211} & \sqrt{2} S_{1222} & \sqrt{2} S_{1233} \end{matrix} & \begin{matrix} 2 S_{2323} & 2 S_{2331} & 2 S_{2312} \\ 2 S_{3123} & 2 S_{3131} & 2 S_{3112} \\ 2 S_{1223} & 2 S_{1231} & 2 S_{1212} \end{matrix} & \begin{matrix} \sqrt{2} S_{2341} & \sqrt{2} S_{2342} & \sqrt{2} S_{2343} \\ \sqrt{2} S_{3141} & \sqrt{2} S_{3142} & \sqrt{2} S_{3143} \\ \sqrt{2} S_{1241} & \sqrt{2} S_{1242} & \sqrt{2} S_{1243} \end{matrix} \\ \begin{matrix} S_{4111} & S_{4122} & S_{4133} \\ S_{4211} & S_{4222} & S_{4233} \\ S_{4311} & S_{4322} & S_{4333} \end{matrix} & \begin{matrix} \sqrt{2} S_{4123} & \sqrt{2} S_{4131} & \sqrt{2} S_{4112} \\ \sqrt{2} S_{4223} & \sqrt{2} S_{4231} & \sqrt{2} S_{4212} \\ \sqrt{2} S_{4323} & \sqrt{2} S_{4331} & \sqrt{2} S_{4312} \end{matrix} & \begin{matrix} S_{4141} & S_{4142} & S_{4143} \\ S_{4241} & S_{4242} & S_{4243} \\ S_{4341} & S_{4342} & S_{4343} \end{matrix} \end{matrix}] [\begin{matrix} \begin{matrix} ε_{11}^{*} \\ ε_{22}^{*} \\ ε_{33}^{*} \end{matrix} \\ \begin{matrix} \sqrt{2} ε_{23}^{*} \\ \sqrt{2} ε_{31}^{*} \\ \sqrt{2} ε_{12}^{*} \end{matrix} \\ \begin{matrix} - E_{1}^{*} \\ - E_{2}^{*} \\ - E_{3}^{*} \end{matrix} \end{matrix}], & (32) \end{matrix}

where the same coefficients are multiplied in the two matrixes. It is important to note that the Eshelby matrix of a piezoelectric and thermoelectric material has dimensions with N/C and V/K in the coupling term. In a piezoelectric material, instead of using SI units for all physical parameters, we can use 1nC = 10⁻⁹C for the charge to avoid numerical problems due to the very large order of the magnitude difference between the elastic constant (~10⁹N/m²) and the coupling constant (~C/N).

For the physical phenomena represented by matrix equation p = Xq in the Mandel notion, the effective property can be obtained by $X^{Eff} = (c_{0} X_{0} + c_{1} X_{1} T) {(c_{0} I + c_{1} T)}^{- 1}$ , where $T = {[I + S X_{0}^{- 1} (X_{1} - X_{0})]}^{- 1}$ and S is the N × N matrix representation of the corresponding Eshelby tensor. We then compare the effective property prediction based on the Mori-Tanaka method with a FEA analysis for the simple case involving spherical inhomogeneities embedded in a transversely isotropic medium for piezoelectric and isotropic media for others, as shown in Figure 3.

FIGURE 3

Figure 3. (A) Normalized Young's modulus of a particle-reinforced composite. The Poisson's ratios (ν) of two phases are identical at 0.25. (B) Effective heat conductivity of the composite. (C) Normalized effective piezoelectric properties of composites. The matrix and reinforcement are PZT-7A and SiC particles, respectively. (D) The effective thermoelectric properties of a Cu-reinforced Bi₂Te₃ composite at 300 K. All properties are normalized with respect to the properties of the matrix. We use the material properties for the piezoelectricity (Odegard, 2004) and the thermoelectricity (Jung et al., 2018) presented in the literature.

Interfacial Imperfections and Anisotropy of the Matrix

Effects of Anisotropy of the Matrix in Various Physical Problems

As mentioned in the previous chapters, we can obtain the effective properties from the Mori-Tanaka method because the expression $X^{Eff} = (c_{0} X_{0} + c_{1} X_{1} T) {(c_{0} I + c_{1} T)}^{- 1}$ is applicable to any arbitrarily anisotropic matrix if the Eshelby tensor is known. The Eshelby tensor for anisotropic medium in various physical problems has been extensively studied in the literature (Mura, 1982; Yu et al., 1994; Huang and Kuo, 1996; Dunn and Wienecke, 1997; Qu and Cherkaoui, 2007; Quang et al., 2011; Martinez-Ayuso et al., 2017; Lee et al., 2018b) and the Eshelby tensors for an ellipsoidal inclusion embedded in an arbitrarily anisotropic medium for elasticity, piezoelectricity, conduction/dielectric phenomena, and thermoelectricity can be obtained numerically from Equations (22, 23, 27, 28), respectively.

Although the Eshelby tensor can be obtained numerically, significant efforts have been devoted to derive the explicit expressions (either closed-form or analytical expression) for the facile application of the homogenization method and to provide better insight into the nature of the tensor. In elasticity, analytic expressions for ellipsoidal shape given in terms of a few integrals are available for spheroidal inclusion in transversely isotropic solids involving five independent elastic constants (or any medium with higher symmetry, i.e., less number of independent elastic constants) (Mura, 1982; Yu et al., 1994). In conduction/dielectric phenomena, analytic solutions have been derived for spheroidal inclusion in isotropic material (Hatta and Taya, 1986) and spherical inclusion in orthotropic and transversely isotropic material (Lee et al., 2018a). For piezoelectricity, Huang and Kuo (1996) suggested the Eshelby tensor expression of spheroidal inclusion in transversely isotropic material with a few integrals. However, for thermoelectricity, to the best of our knowledge, analytic solutions have not been obtained for an anisotropic medium.

Definition of Interfacial Imperfections

In actual composites, the interface between the matrix and an inclusion often has imperfections originating from the manufacturing process or from an inherent lattice mismatch (People and Bean, 1985; Habas et al., 2007). An interfacial imperfection in elasticity refers to debonding or slippage, i.e., a displacement jump (Qu, 1993; Lee and Ryu, 2018; Lee et al., 2019, 2018b), and an interfacial imperfection under thermal conduction (or Kapitza resistance) refers to an abrupt change in the temperature, i.e., a temperature jump (Quang et al., 2011; Lee et al., 2018a). An interfacial imperfection in electrical conduction (i.e., the electrical contact resistance) refers to an abrupt change in the electrical voltage across the interface (Giordano and Palla, 2008). The interfacial imperfection in piezoelectricity considers both displacement jump and the electric potential jump (Wang et al., 2014a,b). Similarly, the interfacial imperfection in thermoelectricity considers the abrupt discontinuities of both temperature and electric potential (Jung et al., 2018).

Out of a few characteristic methods used to describe an interfacial imperfection in elasticity (Qiu and Weng, 1991; Duan et al., 2005), in this work, we consider the interfacial spring model (Qu, 1993) owing to its (i) mathematical simplicity and (ii) mathematical analogy with the interfacial thermal (electrical) resistance in thermal (electrical) conduction. As will be shown later, due to the mathematical similarity, the expressions of the effective stiffness and effective conductivity in the presence of an interfacial imperfection are nearly identical to each other.

It is noteworthy that the interfacial imperfection can be well-treated only for a spherical inclusion, because an ellipsoidal inclusion (even with slight deviation from a sphere) introduces significantly non-uniform interior field (Qu, 1993; Qu and Cherkaoui, 2007; Lee et al., 2019, 2018a,b). In problems involving the elastic deformation, the equality between the normal and tangential interfacial compliances is additionally required to ensure the uniform field inside the inclusion (which is critical for the applicability of the Eshelby tensor and the mean field homogenization). Unfortunately, there exist several studies employing the Mori-Tanaka method to obtain effective physical properties of a composite involving ellipsoidal inhomogeneities in the presence of different normal and tangential compliances (Yang et al., 2013b; Wang et al., 2014a,b; Lee and Ryu, 2018).

Another common mistake in mean-field homogenization studies considering the interfacial damage is the application of the effective physical property expression $X^{Eff} = (c_{0} X_{0} + c_{1} X_{1} T) {(c_{0} I + c_{1} T)}^{- 1}$ and $T = {[I + S X_{0}^{- 1} (X_{1} - X_{0})]}^{- 1}$ by simply replacing the Eshelby tensor S with a modified Eshelby tensor S^M accounting for the interfacial imperfection (Qu, 1993; Barai and Weng, 2011; Yanase and Ju, 2012; Pan et al., 2013; Wang et al., 2014a,b; Shokrieh et al., 2016; Lee and Ryu, 2018). However, in addition to the Eshelby tensor, different expressions for X^Eff and T should be used because the derivation in Section Single inhomogeneity problem and Mean-field homogenization of the effective stiffness (where the problem is solved with the superposition of two sub-problems involving applied and constrained fields, respectively) does not account for the additional contribution from the interface imperfection. We will discuss these issues in detail in the next section.

Prediction of Effective Physical Properties in the Presence of Interfacial Imperfections

We adopt the interface spring model to consider the interfacial damage, as depicted in Figure 4. A displacement jump occurs at the interface due to the spring layer having a vanishing thickness between the infinite matrix and the single inclusion. The spring compliance(η) is represented by α and β for the tangential and normal directions, respectively, and is expressed by Equation (33) in the form of a second-order tensor, as follows:

\begin{matrix} η_{i j} = α δ_{i j} + (β - α) n_{i} n_{j} . & (33) \end{matrix}

FIGURE 4

Figure 4. Schematic of an interface spring model in elasticity. To visualize the interface spring in two directions, we compare undeformed and deformed states.

The constitutive equations and traction(t_i) equilibrium equation at the interface are expressed as follows,

\begin{matrix} \begin{matrix} Δ t_{i} = Δ σ_{i j} n_{j} = [σ_{i j} (\partial Ω^{+}) - σ_{i j} (\partial Ω^{-})] n_{j} = 0 \\ Δ u_{i} = [u_{i} (\partial Ω^{+}) - u_{i} (\partial Ω^{-})] = η_{i j} σ_{j k} n_{k} & (34) \end{matrix} \end{matrix}

where (∂Ω⁺) and (∂Ω⁻) denote the interface on the matrix and the inclusion side, respectively. Formulating the Eshelby inclusion problem by adopting the interfacial condition in Equation (34), the constrained strain is written as Othmani et al. (2011) and Lee et al. (2018b, 2019)

\begin{matrix} ε_{i j}^{c} = S_{i j r s} ε_{r s}^{*} + \frac{1}{2} L_{k l m n} L_{p q r s} \underset{\partial Ω}{\int^{}} η_{k p} {\frac{\partial^{2} G_{i m} (x - y)}{\partial x_{j} \partial y_{n}} + \frac{\partial^{2} G_{j m} (x - y)}{\partial x_{i} \partial y_{n}}} n_{q} (y) n_{l} (y) (ε_{r s}^{c} (y) - ε_{r s}^{*}) d y . & (35) \end{matrix}

Equation (35) is shown to reproduce the perfect interface case with zero spring compliance; i.e., η_ij = 0.

Because Equation (35) is an implicit integral equation involving the constrained strain $ε_{i j}^{c}$ , it is difficult to determine the relationship between $ε_{i j}^{c}$ and $ε_{r s}^{*}$ , i.e., the modified Eshelby tensor. Zhong and Meguid (1997) showed that, when α = β ≡ γ_s and the inclusion shape is a perfect sphere, the strain field inside the spherical inclusion is uniform, as in the case of a perfect interface. Hence, for this special case, the integral equation can be decomposed as follows,

\begin{matrix} ε_{i j}^{c} = S_{i j r s} ε_{r s}^{*} - Γ_{i j r s} (ε_{r s}^{c} - ε_{r s}^{*}) & (36) \end{matrix}

where the fourth-order tensor Γ_ijkl in Equation (36) is defined by

\begin{array}{l} Γ_{i j p q} \equiv - \frac{1}{2} γ_{s} L_{p l m n} L_{p q r s} \underset{\partial Ω}{\int^{}} {\frac{\partial^{2} G_{i m} (x - y)}{\partial x_{j} \partial y_{n}} + \frac{\partial^{2} G_{j m} (x - y)}{\partial x_{i} \partial y_{n}}} \\ \times n_{q} (y) n_{l} (y) d y . & (37) \end{array}

The constrained strain field then can be expressed by a tensor algebraic equation as follows Othmani et al. (2011) and Lee et al. (2018b),

\begin{matrix} ε^{c} = {(I + Γ)}^{- 1} : (S + Γ) : ε^{*} = S^{M} : ε^{*} & (38) \end{matrix}

where S^M is the modified Eshelby tensor. Related to this, in an earlier work by the authors, we proved that for materials with an arbitrary elastic anisotropy, the Γ tensor can be written in terms of the Eshelby tensor for a perfect interface and an elastic stiffness tensor as Lee et al. (2018b)

\begin{matrix} Γ_{i j p q} = \frac{γ_{s}}{R} (I_{i j p q} - S_{i j p q}) L_{p q r s} & (39) \end{matrix}

where R is the radius of the inclusion. Therefore, it is proven that the modified Eshelby The tensor S^M can be written as follows:

\begin{matrix} S^{M} = {[I + \frac{γ_{s}}{R} (I - S) : L]}^{- 1} : [S + \frac{γ_{s}}{R} (I - S) : L] . & (40) \end{matrix}

The numerical validation against a finite element analysis on a triclinic single crystal material (NaAlSi₃O₈) with 21 independent elastic constants (and 36 independent Eshelby tensor components) is presented in Figure 5. This leads to the simplified expression of the modified T tensor (T^M) in the Mori-Tanaka method when the matrix modulus is L₀ and the inhomogeneity modulus is L₁:

\begin{matrix} T^{M} = {(I + S : L_{0}^{- 1} : (L_{1} - L_{0}) + \frac{γ_{s}}{R} (I - S) : L_{1})}^{- 1} & (41) \end{matrix}

FIGURE 5

Figure 5. Thirty-six independent components of a modified Eshelby tensor for a wide range of interfacial damage. The triclinic material used for this calculation is NaAlSi₃O₈ (Lee et al., 2018b).

We derive the effective modulus (L^Eff) of the composite as follows Lee et al. (2019):

\begin{matrix} L^{Eff} = (c_{0} L_{0} + c_{1} L_{1} : T^{M}) : {(c_{0} I + c_{1} T^{M} + c_{1} \frac{γ_{s}}{R} L_{1} : T^{M})}^{- 1} . & (42) \end{matrix}

To confirm the validity of Equation (42), we compare the theoretical prediction with the effective Young's modulus of a composite for a wide range of interface damage as obtained by a FEA (See Figure 6). The equation can also be formulated in terms of the effective inclusion, and we concisely summarize the formulation in Appendix B. We note that many previous studies (Qu, 1993; Barai and Weng, 2011; Yanase and Ju, 2012; Yang et al., 2012; Pan et al., 2013; Shokrieh et al., 2016; Lee and Ryu, 2018) employing the effective stiffness equation in the presence of an interfacial imperfection employ a few incorrect expressions derived by violating Fubini's theorem and relying on inappropriate superposition. A detailed discussion of this issue can be found in work by Dogri's group (Othmani et al., 2011) and by our group (Lee et al., 2018b).

FIGURE 6

Figure 6. Effective Young's modulus of a composite with respect to the normalized interface damage parameter (μ₀γ/R). The Poisson's ratio of the two phases is 0.25 (Lee et al., 2019).

We now turn our attention to other physical phenomena where interfacial imperfections play an important role, specifically the conduction problem. The interfacial thermal resistance α_K (or the Kapitza resistance) is defined as Kapitza (1941) and Quang et al. (2011)

\begin{matrix} T^{out} - T^{in} = - α_{K} q \cdot n, & (43) \end{matrix}

where T^out and Tⁱⁿ refer to the temperatures on the outer and inner surfaces of the interface, respectively, and q is the heat flux at the interface (see Figure 7). The SI unit of the interfacial thermal resistance α_K is [m²K/W]. The interfacial resistance augments an additional surface integration term in the eigen-intensity problem (which corresponds to the eigenstrain problem in elasticity), as follows Quang et al. (2011) and Lee et al. (2018a),

\begin{array}{l} e_{m} (x) = \frac{\partial}{\partial x_{m}} \underset{V}{\int^{}} \frac{\partial G (x - y)}{\partial y_{i}} d y κ_{0}_{i j} e_{j}^{*} + \frac{\partial}{\partial x_{m}} α_{K} κ_{0}_{i j} κ_{0}_{s r} \\ \underset{\partial V}{\int^{}} \frac{\partial G (x - y)}{\partial y_{j}} n_{i} (y) n_{s} (y) (e_{r} (y) - e_{r}^{*} (y)) d y . & (44) \end{array}

FIGURE 7

Figure 7. (A) Schematic of the Kapitza resistance at the interface, and (B) temperature jump across the interface.

It has been found that the heat flux within a spherical inclusion is uniform in the presence of interfacial thermal resistance. In an earlier work of ours, we proved that the modified Eshelby tensor for a matrix with arbitrary anisotropy can be written as follows Lee et al. (2018a),

\begin{matrix} S^{M} = {(I + \frac{α_{K}}{R} (I - S) κ_{0})}^{- 1} (S + \frac{α_{K}}{R} (I - S) κ_{0}), & (45) \end{matrix}

which is similar to the modified Eshelby tensor in elasticity except that the double contraction for the fourth-order tensor is replaced with matrix multiplication for the second-order tensor. Here, I is the second-order identity tensor and S is the Eshelby tensor for the thermal conduction problem. This leads to the simplified expression of the modified concentration tensor (subjected to the constant temperature boundary condition) in the Mori-Tanaka method when the matrix and inhomogeneity thermal conductivity tensors are given as κ₀ and κ₁, respectively:

\begin{matrix} T^{M} = {(I + S κ_{0}^{- 1} (κ_{1} - κ_{0}) + \frac{α_{K}}{R} (I - S) κ_{1})}^{- 1} . & (46) \end{matrix}

We can also obtain the modified localization tensor (subjected to the constant flux boundary condition (q₀)) as

\begin{matrix} q_{1} = B^{M} q_{0}, & (47) \end{matrix}

where q₁ is the heat flux within a single inhomogeneity and $B^{M} = κ_{1} T^{M} κ_{0}^{- 1}$ ; subsequently, we compare the equation with the FEA results for the anisotropic matrix case (see Figure 8; the FEA with the flux boundary condition is more convenient to perform). Based on the modified concentration tensor, we derive the effective conductivity (κ^eff) of the composite as follows Quang et al. (2011) and Lee et al. (2018a),

\begin{matrix} κ^{Eff} = (c_{0} κ_{0} + c_{1} κ_{1} T^{M}) {(c_{0} I + c_{1} T^{M} + c_{1} \frac{α_{K}}{R} κ_{1} T^{M})}^{- 1}, & (48) \end{matrix}

and the theoretical predictions are well-matched with the FEA results (see Figure 9).

FIGURE 8

Figure 8. (A) Normalized heat flux within an inhomogeneity with respect to the interfacial thermal resistance in a single inhomogeneity problem. The inhomogeneity is an isotropic material for which $κ_{1} = [\begin{matrix} 10 \\ 10 \\ 10 \end{matrix}] W / mK$ (B) Heat flux distribution within a matrix and an inhomogeneity on the x, y, and z planes for three different heat flux directions (Lee et al., 2018a).

FIGURE 9

Figure 9. Effective thermal conductivity values [(A) $κ_{11}^{eff}$ , (B) $κ_{22}^{eff},$ and (C) $κ_{33}^{eff}$ ] of a particle-reinforced composite with different interfacial thermal resistances. The material properties of the matrix and particle are $κ_{0} = [\begin{matrix} 1 \\ 2 \\ 3 \end{matrix}] and κ_{1} = [\begin{matrix} 10 \\ 10 \\ 10 \end{matrix}]$ (W/mK), respectively. The radius of the particle used in (A–C) is 1 (mm). The effective thermal conductivity [(D) $κ_{11}^{eff}$ , (E) $κ_{22}^{eff}$ , and (F) $κ_{33}^{eff}$ ] as a function of the radius of the particle under a fixed volume fraction of 5% is also shown (Lee et al., 2018a).

For a thermoelectric composite, we refer to our work on micromechanics homogenization considering both the interfacial thermal and electrical resistance (Jung et al., 2018) where mathematical solutions were obtained for the modified Eshelby tensor, the modified concentration tensor, and the effective thermoelectric properties. Our homogenization study provides new insight on the effect of interfacial thermal and electrical resistance on the thermoelectric properties. In piezoelectricity, the modified Eshelby tensor has been derived by violating the Fubini's theorem by changing integral order (Wang et al., 2014a,b). Hence, the effective piezoelectric property has not been obtained correctly in the presence of interfacial imperfections considering both displacement jump and electric potential jump. With the adaptation of Mandel notation and the careful mathematical derivation, the expression for the effective piezoelectric property is expected to be similar to that of the effective thermoelectric property (Jung et al., 2018). We close the section by providing numerically obtained interior field inside ellipsoidal reinforcement (Figure 10), which clearly demonstrate the limited applicability of the homogenization method in the presence of interfacial imperfection.

FIGURE 10

Figure 10. (A) Heat flux within an ellipsoidal inhomogeneity and (B) strain field within an ellipsoidal inclusion (a_x/a_y = a_x/a_z = 2). The material properties used for the calculations are (A) κ₀ = 1, κ₁ = 10 (W/mK) (Lee et al., 2018a), and (B) E₀ = 200GPa, ν₀ = 0.25 (Lee et al., 2018b).

Summary and Outlook

In this paper, we introduce a series of recent studies which attempted to obtain the effective physical properties of composites by taking into account the interfacial damage and anisotropy of the matrix, providing a universal formulation for different physical phenomena based on their mathematical analogies. First, with the specific example of elasticity, we introduced the concept of the Eshelby tensor and discussed a how single-inclusion problem can be related to a single inhomogeneity problem, which ultimately is applied to mean-field homogenization to predict the effective properties of composites. Second, based on a mathematical analogy, we show that mean-field homogenization equations for different physical phenomena can be represented by linear equations involving symmetric matrices with different dimensions by adapting Mandel's notation. Third, we extend our discussion by taking into account two common issues in realistic applications: imperfections in a matrix-inhomogeneity interface and anisotropy of the matrix.

Although we have discussed the homogenization method while covering a wide range of issues, there remain a few issues requiring further developments. We categorize the advanced issues of homogenization theory into the six problems of (i) various physical phenomena, (ii) anisotropy of the matrix, (iii) matrix-filler interfacial imperfections, (iv) nonlinear responses, (v) time-dependent responses, and (vi) high volume fractions of reinforcements. Although issues (i)-(iii) have been investigated extensively, we nonetheless cannot take into account important non-spherical reinforcement materials (including ellipsoidal types) such as carbon nanotubes and graphene in the presence of interfacial imperfections due to the non-uniform interior field. Although there have been many scaling up attempts to couple micro/nanoscale simulations and homogenization theory to predict the properties of nanocomposites (Yang et al., 2012, 2013a, 2014; Pan et al., 2013; Shin et al., 2013; Shokrieh et al., 2016), they either rely on an incorrect solution stemming from effective stiffness expressions derived from a few problematic bases, such as violations of Fubini's theorem or inappropriate superpositions, or do not recognize the non-uniform interior field (Barai and Weng, 2011; Yang et al., 2012, 2013b; Pan et al., 2013; Shokrieh et al., 2016). There are a few studies of composites showing highly nonlinear (such as hyperelastic) and time-dependent (such as viscoelastic) responses (Friebel et al., 2006; Miled et al., 2011, 2013; Brassart et al., 2012; Doghri et al., 2016), but most lack consideration of interfacial imperfections, which can play an important role on the macroscale composites during repeated cyclic loading and/or in nanocomposites including nanoscale reinforcements with large interfacial fractions. Lastly, to account for the high reinforcement volume fraction regime where the interaction among reinforcements becomes important, several heuristic methods have been developed, such as differential schemes (Hashin, 1988) and/or matrix-filler inversion schemes (Friebel et al., 2006) as well as studies of the lower and upper bounds of the effective properties (Hashin and Shtrikman, 1961), leaving room, however, for further improvements. Studies to overcome aforementioned challenges are underway in our research group with the formation of an international collaborative research network.

Author Contributions

SR designed the structure of the review paper. SR and SL wrote the manuscript. SR and SL. contributed the sections on the micromechanics review and the treatment of anisotropy and interfacial imperfections. SR, SL, and JJ contributed to the summary of theoretical derivations, prepared figures, and discussed the mathematical analogy among different physical problems. JL contributed the theoretical derivations. YK contributed the effective inclusion method section.

Funding

This work is supported by the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (2016M3D1A1900038, 2016R1C1B2011979, and 2018K2A9A2A12000223).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmats.2019.00021/full#supplementary-material

References

Allaoui, A., Bai, S., Cheng, H. M., and Bai, J. B. (2002). Mechanical and electrical properties of a MWNT/epoxy composite. Composites Sci. Technol. 62, 1993–1998. doi: 10.1016/S0266-3538(02)00129-X

REVIEW article

Micromechanics-Based Homogenization of the Effective Physical Properties of Composites With an Anisotropic Matrix and Interfacial Imperfections

Introduction

Eshelby Tensor and Its Application to Effective Stiffness Calculations in Elasticity

Single Inclusion Problem

Single Inhomogeneity Problem

Mean-Field Homogenization of the Effective Stiffness

Mathematical Analogy Between Governing Equations in Various Physical Phenomena

Mathematical Analogy Discussion

Eshelby Tensor for Various Physical Phenomena

Examples of Numerical Calculations and FEA Validation

Interfacial Imperfections and Anisotropy of the Matrix

Effects of Anisotropy of the Matrix in Various Physical Problems

Definition of Interfacial Imperfections

Prediction of Effective Physical Properties in the Presence of Interfacial Imperfections

Summary and Outlook

Author Contributions

Funding

Conflict of Interest Statement

Supplementary Material

References

This article is part of the Research Topic

People also looked at