Micromechanics-based homogenization of the effective physical properties of composites with an anisotropic matrix and interfacial imperfections

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, 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


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., 2018a;Lee et al., 2018b).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., 2018a;Lee et al., 2018b).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 modelled 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-volumefraction 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.
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, we introduce several papers published recently from our group on this issue.
The paper is structured as follows.In section 2, 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 3, 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 4 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.

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  identical to that of the matrix, while an inhomogeneity refers to an embedded material with a different stiffness ′ .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 J. Dvorak, 1992;Tirry and Schryvers, 2009).The Eshelby tensor is defined as the fourth-order tensor  which links the constrained strain within the inclusion  to the eigenstrain  * , as    * (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    in elastostatics, which indicates the displacement in the i th direction at point  by the unit body force in the j th direction at point , in an infinite elastic medium: Here,  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,  can be represented by two independent elastic constants (Section 4 includes a discussion of anisotropy).Green's function for the isotropic material is available in a closed form (Mura, 1982), as follows, where  and  are the shear modulus and the Poisson's ratio of the material, respectively, and | | denotes the standard norm of vector  .
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  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 (), the inclusion exerts traction of   on the matrix (step 4).Due to the constraining effect of the matrix, the inclusion deforms by the constrained strain  , which is different from the eigenstrain  * .The constrained strain field can be expressed by employing Green's function, as follows: (3) Because the eigenstrain and constrained strain are symmetric, the Eshelby tensor has minor symmetry (   ) but no major symmetry   .In an isotropic medium, for an ellipsoidal inclusion with three different semi-axes  ,  , and  , 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    , 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 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 4. We will also describe the Eshelby tensor for an anisotropic medium in Section 4.

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  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 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., 2018a;Lee et al., 2018d).
In the Mandel notation, the strain vector  ⃗, stress vector ⃗, stiffness matrix 〈〉, and Eshelby matrix 〈〉 are defined as follows: The prefactors √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 ,  denote a fourth-order tensor with minor symmetry (    ), and 〈〉, 〈〉 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

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.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 (   ,     :  * , , and   , respectively) are related as follows: :      : , where  :     (where   :    :     ) is the strain concentration tensor for the Mori-Tanaka method. converges to the strain concentration tensor for the single inhomogeneity problem,  , in the limit of very small inhomogeneity volume fraction,  .The effective stiffness can be obtained from the relationship between the volume-averaged stress and strain within the composite (          :     :   and        , respectively) as follows: 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, 1965a;b;Qu and Cherkaoui, 2007).Instead, it requires an iterative numerical solution.

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  in a medium in the absence of an internal source can be written using the following equations, 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   ,   and the second-order tensor   have the physical interpretations given in Table I.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 nondissipative transport phenomena use a second-order derivative with respect to time, leading to electro-magnetic or acoustic wave propagation (Kreyszig, 2006).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 I to obtain the effective electrical and thermal conductivities, the effective dielectric and magnetic permittivities, the effective diffusivity, and the effective fluid permeability.
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 where  is the electrical displacement vector,  is the electric field vector,  is the second rank tensor of the dielectric permittivity constants under constant strain with   , and  is the third rank tensor of the piezoelectric constants under constant stress  with   .We can express Eq. ( 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 (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 where  is obtained by the differentiate  , which is given by In a similar manner, we can simplify the stress and electric displacement to one matrix using the notation Then, the electroelastic moduli and constitutive equation are obtained using the equations below.

𝛴 𝐶 𝑍
On the other hand, for a thermoelectric material, similarly, the constitutive equation is arranged in tensor form, as follows: By using the notation introduced for a piezoelectric material (Jung et al., 2018), the expression is as follows, where Here,   is the electrical current density,   / is the entropy flux,  is the heat flux,  is the temperature,   is the electrical conductivity,  is the Seebeck coefficient,    •   •  is the heat conductivity under a zero electric field, and  is the heat conductivity at zero current.As shown in Eq. ( 21),  contains , which depends on the position  and therefore does not have mathematical similarity with other physics in general.
For a small temperature difference across the hot and cool sides    0 , we can assume that the temperature variable  in the constitutive equation is a constant average temperature  , allowing us to obtain a linear constitutive equation.The validity of this assumption has been studied for wide range of  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  ,  , and  , by simplifying Eq. ( 3), the elastic Eshelby tensor can be expressed as (Mura, 1982), where    ̅  ̅   .Here,    •  •  and  are Green's function and a vector in the Fourier space, respectively.Because   is a homogeneous function of degree 2,     is identical to  ̅  ̅   , where  is a normalized vector.Eq. ( 22) can be evaluated by using the integral variables  and  ; thus,  ̅ 1  cos ,  ̅ 1  sin , and  ̅  given the semi-axes of an ellipsoidal inclusion  .
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), Here, ℎ   ̅  ̅  where    ̅  ̅  .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 Eq. ( 11) and physical property tensor  in Table I, the Eshelby tensor is obtained from 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: 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., 2018c): Or it can be written in an integral involving two variables: where   ̅  ̅  .For a spherical inclusion in an isotropic matrix, the Eshelby tensor can be simplified as   .For a thermoelectric material, with the coupling of electrical and heat conduction, the Eshelby tensor can be written as follows (Jung et al., 2018), where   ̅  ̅  .The Eshelby tensor for a spherical inclusion in an isotropic matrix becomes    .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  , where  is a   matrix and the input (output) field   has  components, the linear operator in elasticity  , conduction (or dielectric)  , piezoelectricity  , and thermoelectricity  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 Eqs.( 30) and ( 31), respectively: As shown by Eqs. ( 30) and ( 31), 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 and 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 1GC 10 C for the charge to avoid numerical problems due to the very large order of the magnitude difference between the elastic constant (~10 9 N/m 2 ) and the coupling constant (~C/N).
For the physical phenomena represented by matrix equation   in the Mandel notion, the effective property can be obtained by           , where          and  is the   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.

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           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., 2018d) 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 Eqs. ( 22), ( 23), ( 27), and (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., 2018b).For piezoelectricity, Jin H. Huang and W. S. Kuo (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 et al., 2018a;Lee et al., 2018d;Lee and Ryu, 2018), 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., 2018b).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;Wang et al., 2014b).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., 2018a;Lee et al., 2018b;Lee et al., 2018d).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;Wang et al., 2014b;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           and          by simply replacing the Eshelby tensor  with a modified Eshelby tensor   accounting for the interfacial imperfection (Qu, 1993;Barai and Weng, 2011;Yanase and Ju, 2012;Pan et al., 2013;Wang et al., 2014a;Wang et al., 2014b;Shokrieh et al., 2016;Lee and Ryu, 2018).However, in addition to the Eshelby tensor, different expressions for  and  should be used because the derivation in Section 2.2. and 2.3.(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 Eq. ( 34) in the form of a second-order tensor, as follows: The constitutive equations and traction  equilibrium equation at the interface are expressed as follows, where ∂Ω and ∂Ω denote the interface on the matrix and the inclusion side, respectively.Formulating the Eshelby inclusion problem by adopting the interfacial condition in Eq. ( 35), the constrained strain is written as (Othmani et al., 2011;Lee et al., 2018a Eq. ( 36) is shown to reproduce the perfect interface case with zero spring compliance; i.e.,  0.
Because Eq. ( 36) is an implicit integral equation involving the constrained strain  , it is difficult to determine the relationship between  and  * , i.e., the modified Eshelby tensor.Zhong, Z and Meguid, S. A. (Zhong and Meguid, 1997) showed that, when   ≡  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, where the fourth-order tensor  in Eq. ( 37) is defined by The constrained strain field then can be expressed by a tensor algebraic equation as follows (Othmani et al., 2011;Lee et al., 2018d), where  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 To confirm the validity of Eq. ( 43), 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., 2018d).
We now turn our attention to other physical phenomena where interfacial imperfections play an important role, specifically the conduction problem.The interfacial thermal resistance  (or the Kapitza resistance) is defined as (Kapitza, 1941;Quang et al., 2011) where  and  refer to the temperatures on the outer and inner surfaces of the interface, respectively, and  is the heat flux at the interface (see Figure 7).The SI unit of the interfacial 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;Lee et al., 2018b), 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., 2018b), 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,  is the second-order identity tensor and  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: We can also obtain the modified localization tensor (subjected to the constant flux boundary condition   ) as where   is the heat flux within a single inhomogeneity and     ; 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  of the composite as follows (Quang et al., 2011;Lee et al., 2018b), and the theoretical predictions are well matched with the FEA results (see Figure 9).
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;Wang et al., 2014b).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.

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 matrixinhomogeneity 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;Pan et al., 2013;Shin et al., 2013;Yang et al., 2013a;Yang et al., 2014;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;Pan et al., 2013;Yang et al., 2013b;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;Brassart et al., 2012;Miled et al., 2013;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.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.

Figures and captions Figure 1 .Figure 2 .
Figures and captions

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 Cureinforced Bi2Te3 composite at 300K.All properties are normalized with respect to the

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.

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

Figure 6 .
Figure 6.Effective Young's modulus of a composite with respect to the normalized interface damage parameter  / .The Poisson's ratio of the two phases is 0.25.

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