We derive an analytical expression to predict the effective pro perties of a particulate- reinforced piezoelectric composite with interfacial imperfectio ns using a micromechanics- based mean–field approach. We co rrectly derive the analytical f ormula of the modified Eshelby tensor, the modified conc entration tensor, and the effe ctive property equations based on the modified Mori–Tanaka method in the presence of interfaci al imperfections. Our results are validated against finite element analyses (FEA) for the ent ire range of interfacial damage levels, from a perfect to a completely disconnected and insulat ed interface. For the facile evaluation of the nontrivial ten sorial equations, we adopt the Mandel notation to perform tensor operations with 9ൈ9 symmetric matrix operations. We apply the method to predict the effective properties of a rep resentative piezoelectric comp osite consisting of PVDF and SiC reinforcements. Keywords Piezoelectric composite, Eshelby t ensor, Interfacial damage, Ef fective modulus 1. Introduction Piezoelectricity refers to the electric charge accumulation in response to an applied mechanical loading, or, conversely, the mechanical strain gener ated from an applied electrical field (Cady, 1964). Such a characteristic has been extensively exploited in many important engineering applications such as high power generation (Lee et al., 2014; Yang et al., 2009), lightweight transducers (Guerin et al., 2019; Stamatellou and K alfas, 2018), precision positioning stages of various microscopes (Li et al., 2019; Wan g et al., 2019), vibration cancellation modules (Lin and Liu, 2006; Tan et al., 2019), and artificial skin of wearable devices (Chorsi et al., 2019; Yua n et al., 2019). Piezoelectric ity occurs in a limited set of materials such as crystals with no inversion symmetry (Jella et al., 2019) and polymer materials with high dipolar moment (Ramadan et al., 2014). For example, one of the widely used piezoelectric polymers, polyvinylidene fluoride (PVDF), ha s five different crystalline phases: 𝛼, 𝛽, 𝛾, δ , and 𝜀. However, the 𝛼 and 𝜀 phases have zero net dipole moment owing to their antiparallel chain structure, and thus, the 𝛽 phase, which has the highest net dipole moment, has been primarily used (Ramadan et al., 2014). Owing t o the limited numbers of piezoelectric materials, it often becomes difficult to obtain a piezoelectric material that simultaneously satisfies the multiple requirements of a given a pplication. For example, the pure PVDF polymer has a small elastic stiffness (~3 GPa), which limits its applicability (Pei and Zeng, 2011). Hence, various piezoelectric composites have been fabricated su ch as the lead zirconate titanate (PZT)–polymer composite (Almusallam et al., 2015; Venkatragavaraj et al., 2001) for an enhanced dielectri c and hydrostatic piezoelectric coefficient and polyvinylidene difluoridein (PVDF)–SiC compos ite (Zhang et al., 2015; Zhou et al., 2015) for improved dielectric permittivity and elastic stiffness. Further, physica l reinforcements as well as voids are known to improve some properties, i.e., the hydrostatic pie zoelectric figure of merit increases with porosity (Zhang et al., 2017). Therefore, for th e efficient design of a piezoelectric composite, it is essential to predict the effecti ve properties of the reinforced piezoelectric composites in ter ms of its size, shape, and volum e fraction of reinforcements. The Mori–Tanaka method, a micromechanics-based mean–field homog enization approaches, has been used to predict various effective physical properties of a reinforced composite, such as the elastic c onstant (Benveniste, 1987; Lee et al., 2019a; Lee and Ryu, 2018), thermal conductivity (Lee et al., 2018; Quang et al., 20 11), dielectric constant (Giordano and Palla, 2008), and thermoelectric (Jung et al., 20 18) and piezoelectric properties (Odegard, 2004; Wang et al., 2014a; Wang et al., 201 4c). The Mori–Tanaka method was originally developed to predict the effective elasti c properties of particulate- reinforced composite by using the Eshelby tensor for a composit e involving multiple inhomogeneities with the mean–fie ld approximation (Benveniste, 1987). The method has also been extended to predict other physical properties based on the mathematical analogy with steady-state governing equations (Ryu et al., 2019). Dunn and T aya were the first to adapt the Mori–Tanaka approach to model the piezoelectric composite by as suming a perfect interface between the reinforcement and t he matrix (Dunn and Taya, 1993a, c). However, a realistic piezoelectric composite interface may exhibit imperfections suc h as the displacement jump (Zhao et al., 2012) or the electri c compliance (Kim et al., 201 1). Such interfacial imperfections can significantly affect the effective piezoelect ric properties of a nanocomposite with nanoscale reinforcements that has a signific antly higher interfacial area compared to a composite involving macroscale reinforcements wit h an identical volume fraction. To consider these interfacial imp erfections, previous studies a ttempted to extend the Mori–Tanaka method by deriving the modified Eshelby tensor and the modified concentration tensor in the presence of interfacial imperfectio ns (Wang et al., 2014a; Wang et al., 2014b; Wang et al., 2014c). However, as shown in the later part of this article, we find that the analytical expressions do not match with the FEA becau se the modified Eshelby tensor was derived by violating the Fubini–Tonelli theorem (Folland, 2013), and th e modified concentration tensor was deriv ed without considering the additi onal displacement jump and the discontinuity in electric potential. Here, we correctly derive analytical expressions of the modifie d Eshelby tensor and the modified concentration tensor of a piezoelectric material t hat are applicable to a material with arbitrary symmetry, and use them to predict the effective properties of piezoelectric composites. For the facile evaluation of the nontrivial tensor equations, we adopt the Mandel notation (Lee et al., 2019b), rather than the V oigt notation, t o perform tensor operations with 9ൈ9 symmetric matrix operations. The analytical results are valida ted against the FEA results. We applied our results to predict the effective elasti c modulus, piezoelectric constant, dielectric constant, and figure of merit map of the PVDF-SiC co mposite. The rest of this paper is organized as follows. First, in Secti on 2, we define the piezoelectric modulus and derive t he modified Eshelby tensor fo r the piezoelectric composite in the presence of interfacial imperfections. We also show that the Mori–Tanaka method is not applicable for non-spheri cal inclusions, comprising elliptical inclusions with an aspect ratio other than 1, because of the non-uniform interior strain and electric field. In Section 3, we derive the modified concentration tensor in terms of the pie zoelectric moduli of matrix and reinforcement, the degree of interfacial imperfection, and the Eshelby tens or. In Section 4, e derive the expression for the effective piezoelectric modulus and obtain the effective dielectric constant, effective elastic modulus, effective piezo electric coefficient, and effective figure of merit for a representative piezoelectric composite co nsisting of PVDF and SiC. We present the concludi ng remarks in Section 5. 2. Modified Eshelby tensor for piezoelectric properties 2.1 Constitutive equation and gov erning equation for piezoelect ric material We choose the elastic strain ൫𝜀 ௜௝൯ and electric field ሺ𝐸௜ሻ as two independent input variables and relate them with t he output variables, mechanical s t r e s s ൫𝜎௜௝൯ and electric displacement ሺ𝐷௜ሻ. The constitutive and governi ng equations of the piezoelectric material are expressed as Eq. (1) and Eq. (2), respectively, 𝜎௜௝ൌ𝐶 ௜௝௠௡ 𝜀௠௡൅𝑒 ௡௜௝ሺെ𝐸 ௡ሻ 𝐷௜ൌ𝑒 ௜௠௡𝜀௠௡െ𝜅 ௜௡ሺെ𝐸 ௡ሻ (1) 𝜎௜௝,௜ൌ0 𝐷௜,௜ൌ0 (2) where 𝐶௜௝௠௡ ,𝑒௡௜௝, and 𝜅௜௡ are the elastic constant, piezoelectric constant, and dielectr ic constant, respectively. The repeated index is a dummy index, wh ich implies the summation over all values, 1 to 3. Eq. (1) can be s hortened by adopting t he notation suggested by Barnett and Lothe (Barnett and Lothe, 1975). 𝛴௜௃ൌ𝐿 ௜௃ெ௡ 𝑍ெ௡ (3) By adapting the notation, one can construct a constitutive equa tion having a similar form with linear elastostatics (Hooke’s l aw) or heat conduction (Fourier’ s law). The repeated capital subscripts are summed over 1–4 and 𝐿௜௃ெ௡ is a piezoelectric moduli tensor which can be constructed by 𝐿௜௃ெ௡ ൌ ⎩⎨⎧𝐶௜௝௠௡ 𝐽, 𝑀 ൌ 1,2,3 𝑒௡௜௝ 𝐽 ൌ 1,2,3, 𝑀 ൌ 4 𝑒௜௠௡ 𝐽 ൌ 4, 𝑀 ൌ 1,2,3 െ𝜅 ௜௡ 𝐽, 𝑀 ൌ 4 (4) 𝑍ெ௡ is the extended strain obtained from the extended displacement 𝑈ெ, which is written as 𝑍ெ௡ൌ𝑑𝑈 ெ 𝜕𝑥௡,𝑈 ெൌ൜𝑢ெ 𝑀 ൌ 1,2,3 𝜙𝑀 ൌ 4 (5) where 𝑢ெ and 𝜙 are the displacement and poten tial field, respectively. Simila rly, the extended stress 𝛴௜௃ is expressed as Σ௜௃ൌ൜𝜎௜௃𝐽 ൌ 1,2,3 𝐷௜𝐽ൌ4. (6) Based on the notation, one can express the governing equation f or the Green’s function of a piezoelectric material as 𝐿௜௃௄௟𝜕ଶ𝐺௄ெሺ𝒙െ𝒚 ሻ 𝜕𝑥௜𝜕𝑥௟൅𝛿 ௃ெ𝛿ሺ𝒙െ𝒚 ሻൌ0 . (7) which is mathematically analogou s with the steady-state governi ng equations of other physical phenomena such as elastostatics (Lee et al., 2019b), h eat conduction (Lee et al., 2018; Quang et al., 2011), and th ermoelectricity (Jung et al., 2018). Following previous studies (Dunn and Taya, 1993a), the Green’s function of a piezo electric material is expressed as 𝐺௄ெሺ𝒙െ𝒚 ሻൌ1 8𝜋|𝒙െ𝒚 |න𝐾௄ெିଵ𝛿൭𝒛∙൬𝒙െ𝒚 |𝒙െ𝒚 |൰൱𝑑𝑆ሺ𝒛ሻ |𝒛|ୀଵ (8) where 𝐾௄ெିଵ is the inverse of 𝐾௄ெൌ𝑧 ௜𝑧௡𝐿௜௄ெ௡ (9) For elasticity and the heat conduction problem, the Green’s fun ctions have been derived in closed form even for an anisotropic material (Lee et al., 2018; Mura, 1982). However, it is difficult to derive the closed form of the Green’s function of a piezoelectric material owing to the mathematical complexity arising from the coupling (piezoele ctric) constant (Huang and Yu, 1994). 2.2 Eshelby’s single inclusion pr oblem in the absence of interf acial imperfection The single inclusion problem of a piezoelectric material can be explained by the following four steps (Fig. 1). When the inclusion with a nonzero extended eigenstrain ሺ𝑍஺௕∗ሻ is released after its insertion into the matrix having a same s ize hole, the inclusion deforms less than the (extended) eigenstrain owing to the interaction w ith the matrix. Based on the mathematical analogy between ela sticity and piezoelectricity, o ne can realize that the constrained extended strain field within the inclusion is unifo rm when the shape of the inclusion is ellipsoid. Therefor e, the constrained (extended) d isplacement and strain field are expressed as 𝑈ெൌන 𝐺 ெ௃ሺ𝒙െ𝒚 ሻ𝛴௤௃∗𝑛௤𝑑𝒚 డஐൌන 𝐺 ெ௃ሺ𝒙െ𝒚 ሻ𝐿௤௃஺௕ 𝑍஺௕∗𝑛௤𝑑𝒚 డஐ 𝑍ெ௡ൌ𝑆 ெ௡஺௕ 𝑍஺௕∗,𝑆 ெ௡஺௕ ൌන 𝑃 ெ௡௃௤ 𝑑𝒚 ஐ𝐿௤௃஺௕ (10) where 𝑃ெ௡௃௤ ൌ ⎩⎪⎨⎪⎧1 2ቆ𝜕ଶ𝐺ெ௃ 𝜕𝑥௡𝜕𝑦௤൅𝜕ଶ𝐺௡௃ 𝜕𝑥ெ𝜕𝑦௤ቇ 𝑀 ൌ 1,2,3 𝜕ଶ𝐺ெ௃ 𝜕𝑥௡𝜕𝑦௤𝑀ൌ4 where 𝜕Ω is the surface of the inclusion and 𝑆ெ௡஺௕ is Eshelby tensor of the piezoelectric material. As suggested by Dunn and Taya (Dunn and Taya, 1993a, c), the Eshelby tensor of a piezoelectric ellipsoid having semi-axes as 𝑎ଵ,𝑎ଶ, and 𝑎ଷ can be expressed in alternative form as 𝑆ெ௡஺௕ ൌ ⎩⎪⎨⎪⎧1 8𝜋𝐿௜௃஺௕ නන ൣ 𝛷 ெ௃௜௡ ሺ𝒛ሻ൅𝛷 ௡௃௜ெ ሺ𝒛ሻ൧𝑑𝜃𝑑𝜉 ଷଶగ ଴ଵ ିଵ𝑀 ൌ 1,2,3 1 4𝜋𝐿௜௃஺௕ නන ൣ 𝛷 ெ௃௜௡ ሺ𝒛ሻ൧𝑑𝜃𝑑𝜉 ଷଶగ ଴ଵ ିଵ𝑀ൌ4 (11) where 𝛷ெ௃௜௡ ሺ𝒛ሻൌ𝑧 ௜𝑧௡𝐾ெ௃ିଵሺ𝒛ሻ and 𝑧ଵൌඥ1െ𝜉ଷଶcos 𝜃 𝑎 ଵൗ, 𝑧 ଶൌඥ1െ𝜉ଷଶsin 𝜃 𝑎 ଶൗ, 𝑧 ଷൌ 𝜉ଷ𝑎ଷ⁄. As for the Green’s function, it is also difficult to obtain t he closed-form expression because of the piezoelectric coupling constants (Huang and Yu, 1994; Mikata, 2001). 2.3 Imperfect interface modell ed by linear spring model Eq. (10) assumes perfect bonding at the interface without displ acement or electric potential discontinuity. We adopt the linear spring model to ac count for the interfacial damage for mechanical response for mathematical simplicity, and consid er a dielectric layer at the interface to model the interfacial damage for electrical respon se (Fig. 2). For a piezoelectric material, the constitutive equation at the damaged interface is expressed as Δ𝑇௃ൌΔ Σ ௜௃𝑛௜ൌቀ Σ ௜௃ሺ𝜕ΩାሻെΣ ௜௃ሺ𝜕Ωିሻቁ𝑛 ௜ൌ0 Δ𝑈 ௃ൌ𝑈 ௃ሺ𝜕Ωାሻെ𝑈 ௃ሺ𝜕Ωିሻൌ𝛬 ௃ெΣ௜ெ𝑛௜ (12) where 𝜕Ωା and 𝜕Ωି denote the interface surface on the matrix and inclusion side, respectively (Fig. 2). Λ௃ெ represents the extended interface spring compliance expressed as 𝛬௃ெൌ൝𝛼𝛿௃ெ൅ሺ𝛽െ𝛼 ሻ𝑛௃𝑛ெ 𝐽 ൌ 1,2,3, 𝑀 ൌ 1,2,3 െ𝛾 𝐽 ൌ 𝑀 ൌ 4 0o t h e r s. (13) 𝛼, 𝛽 represent the tangential and normal direction spring complianc e, respectively , and 𝛾 i s the electric compliance coefficient. 𝛼, 𝛽 are equivalent to the small-deformation-regime compliances of the cohesive zone models that are widely adapted in various FEA studies (Bouhala et al., 2013; Fan and Tadmor, 2019; Springer et al., 2 019). Thus, 𝛼, 𝛽 a r e a l w a y s positive. 𝛾 represents the effect of a dielectric layer at the interface a nd is mathematically analogous with Kapitza’s resistance in a heat conduction proble m (Lee et al., 2018; Quang et al., 2011). As expressed in Eq. (12), large interfacial damage results in large displacement (potential) jump across the in terface and zero interfacial dama ge implies perfect bonding. The extended strain field within the single inclusion having th e interfacial imperfection can be expressed as follows: 𝑍ெ௡ሺ𝒙ሻൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗൅න Δ 𝑈 ௃𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜𝑑𝒚 డஐ ൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗൅න Λ ௃௉Σ௤௉𝑛௤𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜𝑑𝒚 డஐ ൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗൅න Λ ௃௉𝐿௤௉ோ௦ ሺ𝑍ோ௦െ𝑍ோ௦∗ሻ𝑛௤𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜𝑑𝒚 డஐ (14) In general, the constrained strai n field within an ellipsoidal inclusion is non-uniform when the interfacial compliances are finite (Appendix A). Thus, it is di fficult to express the strain field in the closed form because Eq. (14) is an implicit equation inv olving 𝒁 in both the left and right sides. Previous studies (Wang et al., 2014a; Wang et al., 2014b; Wang et al., 2014c) have calculated the vo lume averaged modified Eshelby tensor to consider a non-uniform field, as expressed in Eq. (15). 𝑍̅ெ௡ൌ1 Ωන𝑍 ெ௡ሺ𝒙ሻ ஐ𝑑𝒙 ൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗൅1 ΩනቈනΛ ௃௉𝐿௤௉ோ௦ ሺ𝑍ோ௦െ𝑍ோ௦∗ሻ𝑛௤𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜𝑑𝒚 డஐ቉ ஐ𝑑𝒙 (15) The volume averaged modified Eshelby tensor was simplified afte r changing the order of surface integral and volume integral. A similar simplification has also been used for elasticity (Barai and Weng, 2011; Qu, 1993; Yanase and Ju, 2012) and heat conduction problems (Bonfoh et al., 2017). However, according to the Fubini–Tonelli theorem, it is not legitimate to change the order of two integrals because the Green’s functi on contains a singularity within the integral domain (Fo lland, 2013). We present a repres entative example of the violation of the Fubini–Tonelli theorem by considering a simple 2 D i n t e g r a l p r o b l e m i n Appendix B. In this study, we derive the modified Eshelby tensor without an y approximation. In previous papers on elasticity (Lee et al., 2019b; Zhong and Meg uid, 1997), it has been proved that the strain field within t he spherical inclusion is uniform when 𝛼ൌ𝛽 , and for the heat conduction problem it is known that under finite Kapitza’s resi stance the intensity field (gradient of temperature) is uniform when the shape of the incl usion is spherical (Lee et al., 2018; Quang et al., 2011). Therefore, using a mathematical anal ogy (Ryu et al., 2019), we can deduce that the extended strain field within the piezoelectric inclusion is also uniform when the shape of the inclusion is spherical and 𝛼ൌ𝛽 (Appendix A). For the spherical inclusion problem with interfacial compliance 𝛼ൌ𝛽 , Eq. (14) is simplified as follows: 𝑍ெ௡ሺ𝒙ሻൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗൅Λ ௃௉𝐿௤௉ோ௦ 𝐿௜௃௄௟න𝑃 ெ௡௄௟ 𝑛௜𝑛௤𝑑𝒚 డஐሺ𝑍ோ௦െ𝑍ோ௦∗ሻ ൌ𝑆 ெ௡ோ௦ 𝑍ோ௦∗െ𝛤 ெ௡ோ௦ ሺ𝑍ோ௦െ𝑍ோ௦∗ሻ (16) where 𝛤ெ௡ோ௦ ≡െ න Λ ௃௉𝐿௤௉ோ௦ 𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜𝑛௤𝑑𝒚 డஐ (17) By using 𝑛௤ൌ𝑦 ௤/𝑎 and applying the divergence theorem on 𝛤ெ௡ோ௦ , Eq. (17) is further simplified as 𝛤ெ௡ோ௦ ൌെ න Λ ௃௉𝐿௤௉ோ௦ 𝐿௜௃௄௟𝑃ெ௡௄௟ 𝑛௜ሺ𝒚ሻ𝑛௤ሺ𝒚ሻ𝑑𝒚 డஐ ൌെ Λ ௃௉𝐿௤௉ோ௦ ቈන𝐿 ௜௃௄௟𝑃ெ௡௄௟𝛿௜௤ 𝑎𝑑𝒚 ஐ൅න 𝐸 ௜௃௄௟𝜕𝑃ெ௡௄௟ 𝜕𝑦௜𝑦௤ 𝑎𝑑𝒚 ஐ቉ ൌെ Λ ௃௉𝐿௤௉ோ௦1 𝑎ቈන𝐿 ௟௄௃௤𝑃ெ௡௄௟ 𝑑𝒚 ஐ൅න 𝐸 ௜௃௄௟𝜕𝑃ெ௡௄௟ 𝜕𝑦௜𝑦௤𝑑𝒚 ஐ቉ (18) where 𝑎 is the radius of the inclusion. Noting that the first integral term on the right hand side is the definition of the Eshelby tensor shown in Eq. (10), one can write the Eq. (18) as follows: 𝛤ெ௡ோ௦ ൌെ Λ ௃௉𝐸௤௉ோ௦1 𝑎ቈන𝐸 ௟௄௃௤𝑃ெ௡௄௟ 𝑑𝒚 ஐ൅න 𝐸 ௜௃௄௟𝜕𝑃ெ௡௄௟ 𝜕𝑦௜𝑦௤𝑑𝒚 ஐ቉ ൌെ Λ ௃௉𝐸௤௉ோ௦1 𝑎ቈ𝑆ெ௡௃௤ ൅න 𝐸 ௜௃௄௟𝜕𝑃ெ௡௄௟ 𝜕𝑦௜𝑦௤𝑑𝒚 ஐ቉ (19) To further simplify Eq. (19), we consider the following equatio n which is obtained by multiplying 𝑦௤ after differentiating Eq. (7) with respect to 𝑦௡. 𝐿௜௃௄௟𝜕ଷ𝐺௄ெሺ𝒙െ𝒚 ሻ 𝜕𝑦௡𝜕𝑦௜𝜕𝑦௟𝑦௤൅𝜕 𝜕𝑦௡ቀ𝛿௃ெ𝛿ሺ𝒙െ𝒚 ሻቁ𝑦௤ൌ0 (20) Applying the divergence theorem a fter integrating Eq. (20) for the inclusion volume, we obtain න𝐿 ௜௃௄௟𝜕ଷ𝐺௄ெሺ𝒙െ𝒚 ሻ 𝜕𝑦௡𝜕𝑦௜𝜕𝑦௟𝑦௤𝑑𝒚 ஐൌെ න𝜕 𝜕𝑦௡ቀ𝛿௃ெ𝛿ሺ𝒙െ𝒚 ሻቁ𝑦௤𝑑𝒚 ஐ ൌെ න ቀ 𝛿 ௃ெ𝛿ሺ𝒙െ𝒚 ሻቁ𝜕𝑦௤ 𝜕𝑦௡𝑑𝒚 ஐ ൌ𝛿 ௃ெ𝛿௡௤. (21) Therefore, න𝐿 ௜௃௄௟𝜕𝑃ெ௡௄௟ 𝜕𝑦௜𝑦௤𝑑𝒚 ஐൌቐെ1 2൫𝛿௝௠𝛿௡௤൅𝛿 ௝௡𝛿௠௚൯ 𝑀 ൌ 1,2,3 െ𝛿 ௃ெ𝛿௡௤ 𝑀ൌ4 (22) After putting Eq. (22) into Eq. (19), we can express 𝛤ெ௡ோ௦ a s 𝛤ெ௡ோ௦ ൌΛ ௃௉𝐿௤௉ோ௦1 𝑎൫𝐼ெ௡௃௤ െ𝑆 ெ௡௃௤ ൯. (23) where 𝐼ெ௡௃௤ ൌቐ1 2൫𝛿௝௠𝛿௡௤൅𝛿 ௝௡𝛿௠௚൯ 𝑀 ൌ 1,2,3 𝛿௃ெ𝛿௡௤ 𝑀ൌ4 When 𝛼ൌ𝛽 , the extended interfacial compliance Λ௃௉ is simplified as Eq. (24), and thus, Eq. (23) can be rewritten as Eq. (25) Λ௃௉ൌ൜𝛼𝛿௝௣ 𝐽 ൌ 𝑃 ൌ 1,2,3 െ𝛾𝛿 ௃௉𝐽ൌ𝑃ൌ4 (24) 𝛤ெ௡ோ௦ ൌ෍𝜔௃𝛿௃௉𝐿௤௉ோ௦1 𝑎൫𝐼ெ௡௃௤ െ𝑆 ெ௡௃௤ ൯ସ ௃ୀଵ ൌ ෍𝜔௃ 𝑎൫𝐼ெ௡௃௤ െ𝑆 ெ௡௃௤ ൯𝐿௤௃ோ௦ସ ௃ୀଵ (25) where 𝜔௃ൌ൜𝛼 𝐽 ൌ 1,2,3 െ𝛾 𝐽 ൌ 4 Eq. (25) is mathematically analogous with 𝜞 in other physical phenomena such as elasticity (Lee et al., 2019a; Lee et al., 2019b; Othmani et al., 2011), h eat conduction (Lee et al., 2018; Quang et al., 2011), and thermoelectricity (Jung et al., 2018). Combining Eq. (16) and (25), we can predict the modified Eshelby tensor as 𝒁ൌ𝑺ெ:𝒁∗ 𝑺ெൌሺ𝑰൅𝜞 ሻିଵ:ሺ𝑺൅𝜞 ሻ, (26) where : denotes a double contraction which implies 𝑺ெ:𝒁∗ൌ𝑆ெ௡ோ௦ெ𝑍ோௌ∗. 2.4 Mandel notation To evaluate Eq. (26), the inverse and double contraction of 4ൈ3ൈ4ൈ3 o r 3ൈ 4ൈ4ൈ3 tensors is required. Most papers concerning the micromechanics study of a piezoelectric material have adopted the V oigt notation to carry o u t t h e t e n s o r o p e r a t i o n s (Martinez-Ayuso et al., 2017; Odegard, 2004; Wang et al., 2014c ). Here, we adopt the Mandel notation because it is impossible to correctly calculate the mo dified Eshelby tensor with the V oigt notation. When the V oigt notation is adopted, the extende d stress and strain vector expressed in Eq. (27) and the coefficient multiplied on compone nts in the transformed 9ൈ9 matrix depend on the type of the tensor. 〈𝚺〉୚୭୧୥୲ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝜎11 𝜎22 𝜎33 𝜎23 𝜎31 𝜎12 𝐷1 𝐷2 𝐷3⎦⎥⎥⎥⎥⎥⎥⎥⎤ , 〈𝒁〉୚୭୧୥୲ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝜀11 𝜀22 𝜀33 2𝜀23 2𝜀31 2𝜀12 െ𝐸1 െ𝐸2 െ𝐸3⎦⎥⎥⎥⎥⎥⎥⎥⎤ (27) For example, the transformed piezoelectric tensor 〈𝑳〉୚୭୧୥୲, Eshelby tensor 〈𝑺〉୚୭୧୥୲, and inverse of piezoelectric tensor 〈𝑴〉୚୭୧୥୲ൌ〈𝑳ି𝟏〉୚୭୧୥୲ are expressed as 〈𝑳〉୚୭୧୥୲ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝐿1111 𝐿1122 𝐿1133 𝐿2211 𝐿2222 𝐿2233 𝐿3311 𝐿3322 𝐿3333𝐿1123 𝐿1131 𝐿1112 𝐿2223 𝐿2231 𝐿2212 𝐿3323 𝐿3331 𝐿3312𝐿1141 𝐿1142 𝐿1143 𝐿2241 𝐿2242 𝐿2243 𝐿3341 𝐿3342 𝐿3343 𝐿2311 𝐿2322 𝐿2333 𝐿3111 𝐿3122 𝐿3133 𝐿1211 𝐿1222 𝐿1233𝐿2323 𝐿2331 𝐿2312 𝐿3123 𝐿3131 𝐿3112 𝐿1223 𝐿1231 𝐿1212𝐿2341 𝐿2342 𝐿2343 𝐿3141 𝐿3142 𝐿3143 𝐿1241 𝐿1242 𝐿1243 𝐿1411 𝐿1422 𝐿1433 𝐿2411 𝐿2422 𝐿2433 𝐿3411 𝐿3422 𝐿3433𝐿1423 𝐿1431 𝐿1412 𝐿2423 𝐿2431 𝐿2412 𝐿3423 𝐿3431 𝐿3412𝐿1441 𝐿1442 𝐿1443 𝐿2441 𝐿2442 𝐿2443 𝐿3441 𝐿3442 𝐿3443⎦⎥⎥⎥⎥⎥⎥⎥⎤ (28) 〈𝑺〉୚୭୧୥୲ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝑆1111 𝑆1122 𝑆1133 𝑆2211 𝑆2222 𝑆2233 𝑆3311 𝑆3322 𝑆3333𝑆1123 𝑆1131 𝑆1112 𝑆2223 𝑆2231 𝑆2212 𝑆3323 𝑆3331 𝑆3312𝑆1141 𝑆1142 𝑆1143 𝑆2241 𝑆2242 𝑆2243 𝑆3341 𝑆3342 𝑆3343 2𝑆2311 2𝑆2322 2𝑆2333 2𝑆3111 2𝑆3122 2𝑆3133 2𝑆1211 2𝑆1222 2𝑆12332𝑆2323 2𝑆2331 2𝑆2312 2𝑆3123 2𝑆3131 2𝑆3112 2𝑆1223 2𝑆1231 2𝑆12122𝑆2341 2𝑆2342 2𝑆2343 2𝑆3141 2𝑆3142 2𝑆3143 2𝑆1241 2𝑆1242 2𝑆1243 𝑆4111 𝑆4122 𝑆4133 𝑆4211 𝑆4222 𝑆4233 𝑆4311 𝑆4322 𝑆4333𝑆4123 𝑆4131 𝑆4112 𝑆4223 𝑆4231 𝑆4212 𝑆4323 𝑆4331 𝑆4312𝑆4141 𝑆4142 𝑆4143 𝑆4241 𝑆4242 𝑆4243 𝑆4341 𝑆4342 𝑆4343 ⎦⎥⎥⎥⎥⎥⎥⎥⎤ 〈𝑴〉୚୭୧୥୲ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝑀ଵଵଵଵ 𝑀ଵଵଶଶ 𝑀ଵଵଷଷ 𝑀ଶଶଵଵ 𝑀ଶଶଶଶ 𝑀ଶଶଷଷ 𝑀ଷଷଵଵ 𝑀ଷଷଶଶ 𝑀ଷଷଷଷ2𝑀 ଵଵଶଷ 2𝑀 ଵଵଷଵ 2𝑀 ଵଵଵଶ 2𝑀 ଶଶଶଷ 2𝑀 ଶଶଷଵ 2𝑀 ଶଶଵଶ 2𝑀 ଷଷଶଷ 2𝑀 ଷଷଷଵ 2𝑀 ଷଷଵଶ𝑀ଵଵସଵ 𝑀ଵଵସଶ 𝑀ଵଵସଷ 𝑀ଶଶସଵ 𝑀ଶଶସଶ 𝑀ଶଶସଷ 𝑀ଷଷସଵ 𝑀ଷଷସଶ 𝑀ଷଷସଷ 2𝑀 ଶଷଵଵ 2𝑀 ଶଷଶଶ 2𝑀 ଶଷଷଷ 2𝑀 ଷଵଵଵ 2𝑀 ଷଵଶଶ 2𝑀 ଷଵଷଷ 2𝑀 ଵଶଵଵ 2𝑀 ଵଶଶଶ 2𝑀 ଵଶଷଷ4𝑀 ଶଷଶଷ 4𝑀 ଶଷଷଵ 4𝑀 ଶଷଵଶ 4𝑀 ଷଵଶଷ 4𝑀 ଷଵଷଵ 4𝑀 ଷଵଵଶ 4𝑀 ଵଶଶଷ 4𝑀 ଵଶଷଵ 4𝑀 ଵଶଵଶ2𝑀 ଶଷସଵ 2𝑀 ଶଷସଶ 2𝑀 ଶଷସଷ 2𝑀 ଷଵସଵ 2𝑀 ଷଵସଶ 2𝑀 ଷଵସଷ 2𝑀 ଵଶସଵ 2𝑀 ଵଶସଶ 2𝑀 ଵଶସଷ 𝑀ଵସଵଵ 𝑀ଵସଶଶ 𝑀ଵସଷଷ 𝑀ଶସଵଵ 𝑀ଶସଶଶ 𝑀ଶସଷଷ 𝑀ଷସଵଵ 𝑀ଷସଶଶ 𝑀ଷସଷଷ2𝑀 ଵସଶଷ 2𝑀 ଵସଷଵ 2𝑀 ଵସଵଶ 2𝑀 ଶସଶଷ 2𝑀 ଶସଷଵ 2𝑀 ଶସଵଶ 2𝑀 ଷସଶଷ 2𝑀 ଷସଷଵ 2𝑀 ଷସଵଶ𝑀ଵସସଵ 𝑀ଵସସଶ 𝑀ଵସସଷ 𝑀ଶସସଵ 𝑀ଶସସଶ 𝑀ଶସସଷ 𝑀ଷସସଵ 𝑀ଷସସଶ 𝑀ଷସସଷ ⎦⎥⎥⎥⎥⎥⎥⎥⎤ If the Mandel notation is used, t he extended stress and strain vector is written as 〈𝚺〉୑ୟ୬ୢୣ୪ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡𝜎11 𝜎22 𝜎33 √2𝜎23 √2𝜎31 √2𝜎12 𝐷1 𝐷2 𝐷3⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤ , 〈𝒁〉୑ୟ୬ୢୣ୪ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡𝜀11 𝜀22 𝜀33 √2𝜀23 √2𝜀31 √2𝜀12 െ𝐸1 െ𝐸2 െ𝐸3⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤ (29) and 〈𝑳〉୑ୟ୬ୢୣ୪,〈𝑺〉୑ୟ୬ୢୣ୪,〈𝑴〉୑ୟ୬ୢୣ୪ is expressed as 〈𝑳〉୑ୟ୬ୢୣ୪ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ 𝐿ଵଵଵଵ 𝐿ଵଵଶଶ 𝐿ଵଵଷଷ 𝐿ଶଶଵଵ 𝐿ଶଶଶଶ 𝐿ଶଶଷଷ 𝐿ଷଷଵଵ 𝐿ଷଷଶଶ 𝐿ଷଷଷଷ√2𝐿ଵଵଶଷ √2𝐿ଵଵଷଵ √2𝐿ଵଵଵଶ √2𝐿ଶଶଶଷ √2𝐿ଶଶଷଵ √2𝐿ଶଶଵଶ √2𝐿ଷଷଶଷ √2𝐿ଷଷଷଵ √2𝐿ଷଷଵଶ𝐿ଵଵସଵ 𝐿ଵଵସଶ 𝐿ଵଵସଷ 𝐿ଶଶସଵ 𝐿ଶଶସଶ 𝐿ଶଶସଷ 𝐿ଷଷସଵ 𝐿ଷଷସଶ 𝐿ଷଷସଷ √2𝐿ଶଷଵଵ √2𝐿ଶଷଶଶ √2𝐿ଶଷଷଷ √2𝐿ଷଵଵଵ √2𝐿ଷଵଶଶ √2𝐿ଷଵଷଷ √2𝐿ଵଶଵଵ √2𝐿ଵଶଶଶ √2𝐿ଵଶଷଷ2𝐿ଶଷଶଷ 2𝐿ଶଷଷଵ 2𝐿ଶଷଵଶ 2𝐿ଷଵଶଷ 2𝐿ଷଵଷଵ 2𝐿ଷଵଵଶ 2𝐿ଵଶଶଷ 2𝐿ଵଶଷଵ 2𝐿ଵଶଵଶ√2𝐿ଶଷସଵ √2𝐿ଶଷସଶ √2𝐿ଶଷସଷ √2𝐿ଷଵସଵ √2𝐿ଷଵସଶ √2𝐿ଷଵସଷ √2𝐿ଵଶସଵ √2𝐿ଵଶସଶ √2𝐿ଵଶସଷ 𝐿ଵସଵଵ 𝐿ଵସଶଶ 𝐿ଵସଷଷ 𝐿ଶସଵଵ 𝐿ଶସଶଶ 𝐿ଶସଷଷ 𝐿ଷସଵଵ 𝐿ଷସଶଶ 𝐿ଷସଷଷ√2𝐿ଵସଶଷ √2𝐿ଵସଷଵ √2𝐿ଵସଵଶ √2𝐿ଶସଶଷ √2𝐿ଶସଷଵ √2𝐿ଶସଵଶ √2𝐿ଷସଶଷ √2𝐿ଷସଷଵ √2𝐿ଷସଵଶ𝐿ଵସସଵ 𝐿ଵସସଶ 𝐿ଵସସଷ 𝐿ଶସସଵ 𝐿ଶସସଶ 𝐿ଶସସଷ 𝐿ଷସସଵ 𝐿ଷସସଶ 𝐿ଷସସଷ ⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ (30) 〈𝑺〉୑ୟ୬ୢୣ୪ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ 𝑆ଵଵଵଵ 𝑆ଵଵଶଶ 𝑆ଵଵଷଷ 𝑆ଶଶଵଵ 𝑆ଶଶଶଶ 𝑆ଶଶଷଷ 𝑆ଷଷଵଵ 𝑆ଷଷଶଶ 𝑆ଷଷଷଷ√2𝑆ଵଵଶଷ √2𝑆ଵଵଷଵ √2𝑆ଵଵଵଶ √2𝑆ଶଶଶଷ √2𝑆ଶଶଷଵ √2𝑆ଶଶଵଶ √2𝑆ଷଷଶଷ √2𝑆ଷଷଷଵ √2𝑆ଷଷଵଶ𝑆ଵଵସଵ 𝑆ଵଵସଶ 𝑆ଵଵସଷ 𝑆ଶଶସଵ 𝑆ଶଶସଶ 𝑆ଶଶସଷ 𝑆ଷଷସଵ 𝑆ଷଷସଶ 𝑆ଷଷସଷ √2𝑆ଶଷଵଵ √2𝑆ଶଷଶଶ √2𝑆ଶଷଷଷ √2𝑆ଷଵଵଵ √2𝑆ଷଵଶଶ √2𝑆ଷଵଷଷ √2𝑆ଵଶଵଵ √2𝑆ଵଶଶଶ √2𝑆ଵଶଷଷ2𝑆ଶଷଶଷ 2𝑆ଶଷଷଵ 2𝑆ଶଷଵଶ 2𝑆ଷଵଶଷ 2𝑆ଷଵଷଵ 2𝑆ଷଵଵଶ 2𝑆ଵଶଶଷ 2𝑆ଵଶଷଵ 2𝑆ଵଶଵଶ√2𝑆ଶଷସଵ √2𝑆ଶଷସଶ √2𝑆ଶଷସଷ √2𝑆ଷଵସଵ √2𝑆ଷଵସଶ √2𝑆ଷଵସଷ √2𝑆ଵଶସଵ √2𝑆ଵଶସଶ √2𝑆ଵଶସଷ 𝑆ସଵଵଵ 𝑆ସଵଶଶ 𝑆ସଵଷଷ 𝑆ସଶଵଵ 𝑆ସଶଶଶ 𝑆ସଶଷଷ 𝑆ସଷଵଵ 𝑆ସଷଶଶ 𝑆ସଷଷଷ√2𝑆ସଵଶଷ √2𝑆ସଵଷଵ √2𝑆ସଵଵଶ √2𝑆ସଶଶଷ √2𝑆ସଶଷଵ √2𝑆ସଶଵଶ √2𝑆ସଷଶଷ √2𝑆ସଷଷଵ √2𝑆ସଷଵଶ𝑆ସଵସଵ 𝑆ସଵସଶ 𝑆ସଵସଷ 𝑆ସଶସଵ 𝑆ସଶସଶ 𝑆ସଶସଷ 𝑆ସଷସଵ 𝑆ସଷସଶ 𝑆ସଷସଷ ⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ 〈𝑴〉୑ୟ୬ୢୣ୪ൌ ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡ 𝑀ଵଵଵଵ 𝑀ଵଵଶଶ 𝑀ଵଵଷଷ 𝑀ଶଶଵଵ 𝑀ଶଶଶଶ 𝑀ଶଶଷଷ 𝑀ଷଷଵଵ 𝑀ଷଷଶଶ 𝑀ଷଷଷଷ√2𝑀ଵଵଶଷ √2𝑀ଵଵଷଵ √2𝑀ଵଵଵଶ √2𝑀ଶଶଶଷ √2𝑀ଶଶଷଵ √2𝑀ଶଶଵଶ √2𝑀ଷଷଶଷ √2𝑀ଷଷଷଵ √2𝑀ଷଷଵଶ𝑀ଵଵସଵ 𝑀ଵଵସଶ 𝑀ଵଵସଷ 𝑀ଶଶସଵ 𝑀ଶଶସଶ 𝑀ଶଶସଷ 𝑀ଷଷସଵ 𝑀ଷଷସଶ 𝑀ଷଷସଷ √2𝑀ଶଷଵଵ √2𝑀ଶଷଶଶ √2𝑀ଶଷଷଷ √2𝑀ଷଵଵଵ √2𝑀ଷଵଶଶ √2𝑀ଷଵଷଷ √2𝑀ଵଶଵଵ √2𝑀ଵଶଶଶ √2𝑀ଵଶଷଷ2𝑀 ଶଷଶଷ 2𝑀 ଶଷଷଵ 2𝑀 ଶଷଵଶ 2𝑀 ଷଵଶଷ 2𝑀 ଷଵଷଵ 2𝑀 ଷଵଵଶ 2𝑀 ଵଶଶଷ 2𝑀 ଵଶଷଵ 2𝑀 ଵଶଵଶ√2𝑀ଶଷସଵ √2𝑀ଶଷସଶ √2𝑀ଶଷସଷ √2𝑀ଷଵସଵ √2𝑀ଷଵସଶ √2𝑀ଷଵସଷ √2𝑀ଵଶସଵ √2𝑀ଵଶସଶ √2𝑀ଵଶସଷ 𝑀ଵସଵଵ 𝑀ଵସଶଶ 𝑀ଵସଷଷ 𝑀ଶସଵଵ 𝑀ଶସଶଶ 𝑀ଶସଷଷ 𝑀ଷସଵଵ 𝑀ଷସଶଶ 𝑀ଷସଷଷ√2𝑀ଵସଶଷ √2𝑀ଵସଷଵ √2𝑀ଵସଵଶ √2𝑀ଶସଶଷ √2𝑀ଶସଷଵ √2𝑀ଶସଵଶ √2𝑀ଷସଶଷ √2𝑀ଷସଷଵ √2𝑀ଷସଵଶ𝑀ଵସସଵ 𝑀ଵସସଶ 𝑀ଵସସଷ 𝑀ଶସସଵ 𝑀ଶସସଶ 𝑀ଶସସଷ 𝑀ଷସସଵ 𝑀ଷସସଶ 𝑀ଷସସଷ ⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ , As shown in Eq. (30), all the transformed matrices have identic al coefficients, √2 o r 2, regardless of the tensor type. Us ing the Mandel notation, one c an evaluate Eq. (26) by carrying out 9ൈ9 matrix calculations. Thus, the modified Eshelby tensor can be obtained as 〈𝑺ெ〉୑ୟ୬ୢୣ୪ൌ൫〈𝑰〉୑ୟ୬ୢୣ୪൅〈𝜞〉୑ୟ୬ୢୣ୪൯ିଵ൫〈𝑺〉୑ୟ୬ୢୣ୪൅〈𝜞〉୑ୟ୬ୢୣ୪൯ (31) It is notable that the Eq. (31) can not be obtained by using th e V oigt notation. For example, let 𝑸ൌ𝑺 : 𝑳 ൫𝑄ெ௡ோ௦ ൌ𝑆 ெ௡௃௤ 𝐿௤௃ோ௦൯, then 𝑄ଶଷଶଷ is expressed as follows: 𝑄ଶଷଶଷ ൌ𝑆 ଶଷଵଵ 𝐿ଵଵଶଷ ൅𝑆 ଶଷଶଶ 𝐿ଶଶଶଷ ൅𝑆 ଶଷଷଷ 𝐿ଷଷଶଷ ൅2 𝑆 ଶଷଶଷ 𝐿ଶଷଶଷ ൅2 𝑆 ଶଷଷଵ 𝐿ଷଵଶଷ ൅2 𝑆 ଶଷଵଶ 𝐿ଵଶଶଷ ൅𝑆 ଶଷସଵ 𝐿ଵସଶଷ ൅𝑆 ଶଷସଶ 𝐿ଶସଶଷ ൅𝑆 ଶଷସଷ 𝐿ଷସଶଷ (32) Eq. (32) is rewritten by using the V oigt and Mandel notation an d Eq. (34) satisfies 〈𝑸〉ସସ୑ୟ୬ୢୣ୪ൌ2 𝑄 ଶଷଶଷ ൌ〈𝑺〉ସ௄୑ୟ୬ୢୣ୪〈𝑳〉௄ସ୑ୟ୬ୢୣ୪ whereas Eq. (33) does not satisfy 𝑄ଶଷଶଷ ൌ 𝑘〈𝑺〉ସ௄୚୭୧୥୲〈𝑳〉௄ସ୚୭୧୥୲, where 𝑘 is constant. 𝑄ଶଷଶଷ ൌ1 2〈𝑺〉ସଵ୚୭୧୥୲〈𝑳〉ଵସ୚୭୧୥୲൅1 2〈𝑺〉ସଶ୚୭୧୥୲〈𝑳〉ଶସ୚୭୧୥୲൅1 2〈𝑺〉ସଷ୚୭୧୥୲〈𝑳〉ଷସ୚୭୧୥୲ ൅〈𝑺〉ସସ୚୭୧୥୲〈𝑳〉ସସ୚୭୧୥୲൅〈𝑺〉ସହ୚୭୧୥୲〈𝑳〉ହସ୚୭୧୥୲൅〈𝑺〉ସ଺୚୭୧୥୲〈𝑳〉଺ସ୚୭୧୥୲ ൅1 2〈𝑳〉ସ଻୚୭୧୥୲〈𝑺〉଻ସ୚୭୧୥୲൅1 2〈𝑳〉ସ଼୚୭୧୥୲〈𝑺〉଼ସ୚୭୧୥୲൅1 2〈𝑳〉ସଽ୚୭୧୥୲〈𝑺〉ଽସ୚୭୧୥୲ (33) 𝑄ଶଷଶଷ ൌ1 2〈𝑺〉ସଵ୑ୟ୬ୢୣ୪〈𝑳〉ଵସ୑ୟ୬ୢୣ୪൅1 2〈𝑺〉ସଶ୑ୟ୬ୢୣ୪〈𝑳〉ଶସ୑ୟ୬ୢୣ୪൅1 2〈𝑺〉ସଷ୑ୟ୬ୢୣ୪〈𝑳〉ଷସ୑ୟ୬ୢୣ୪ ൅1 2〈𝑺〉ସସ୑ୟ୬ୢୣ୪〈𝑳〉ସସ୑ୟ୬ୢୣ୪൅1 2〈𝑺〉ସହ୑ୟ୬ୢୣ୪〈𝑳〉ହସ୑ୟ୬ୢୣ୪ ൅1 2〈𝑺〉ସ଺୑ୟ୬ୢୣ୪〈𝑳〉଺ସ୑ୟ୬ୢୣ୪൅1 2〈𝑳〉ସ଻୑ୟ୬ୢୣ୪〈𝑺〉଻ସ୑ୟ୬ୢୣ୪ ൅1 2〈𝑳〉ସ଼୑ୟ୬ୢୣ୪〈𝑺〉଼ସ୑ୟ୬ୢୣ୪൅1 2〈𝑳〉ସଽ୑ୟ୬ୢୣ୪〈𝑺〉ଽସ୑ୟ୬ୢୣ୪ ൌ1 2〈𝑺〉ସ௄୑ୟ୬ୢୣ୪〈𝑳〉௄ସ୑ୟ୬ୢୣ୪ (34) W e c a l c u l a t e t h e m o d i f i e d E s h e l b y t e n s o r f o r P Z T - 5 A w i t h t h e m a terial constants listed in Table 1. When conventional units are used, (the unit of the elastic constants is GPa = 109 Pa and that of the piezoelectric constants is C/m2), numerical errors are likely to accumulate during the matrix inverse operation due to the signi ficant difference in the order of magnitudes between the two constants. To avoid the numerical error, we set the unit of charge as nC = 10-9 C (Ryu et al., 2019). 2.5 Numerical validation T o v a l i d a t e t h e m o d i f i e d E s h e l b y t e n s o r r e s u l t s , w e c a r r y o u t a series of FEA simulations using COMSOL (Comsol, 2015). We use a single inclus ion with diameter of 2 mm and a matrix in cubic shape with edge length of 20 mm (the e dge length is ten times longer than the diameter of the i nclusion to reduce the effect of the finite size as seen as Fig. 3). We use approximately 400,000 and 8,000 3D linear tetrahedro n elements to construct the mesh inside the matrix and inclusion, respectively. The unit ei genstrain is assigned on the inclusion under the fixed outer s urfaces of the matrix, and we then predict the Eshelby tensor using the constrained str ain and electric field. We predict the modified Eshelby tensor for a wide range of inte rfacial ranges and compare the result with FEA calculations as depicted in Fig. 4. F i g . 4 ( a ) – ( c ) s h o w s t h e modified Eshelby tensor when only elastic interfacial damage is considered ሺ𝛼ൌ𝛽്0 , 𝛾ൌ 0ሻ. As 𝛼 goes to infinite, the traction converges to zero, so the inclu sion deforms by eigenstrain ሺ𝒁ൌ𝒁∗ሻ, which implies that 𝑆ଵଵଵଵெൌ1 ,𝑆ଶଶଵଵெൌ𝑆ସଷଵଵெൌ0. Further, when the electric interfacial damage incr eases in the absence of an elas tic interfacial damage, 𝑆ସଷଵଵெ decreases and converges to a finite value at an infinite electr ic interfacial damage (Fig. 4(d)). Even in the infinite electric interfacial damage limit, the str ain field around the inclusion is finite because there is no elastic interfacial damage. Hence, t he finite electric field within the inclusion is induced by the piezoelectric effect, i.e., it is d ue to the charge accumulated at the matrix side of the interface. With the increase in the elastic interfacial damage (the α଴ൌ 1.21 case in Fig. 4(d)), the electric field within the inclusion re duces because of the reduction in the strain field around the inclusion. We note tha t, in the limit of the infinite interfacial damage for both (i.e., 𝛼ൌ𝛽→∞ , 𝛾→∞ ), all diagonal components converge to 1 and the off-diagonal terms go to zero because of the complete mechanical and electric decoupling between the inclusion and the matrix. The result can be evidently explained with Eq. (31). When the interfacial damage increases, 𝜞 goes to infinite and 〈𝑺ெ〉୑ୟ୬ୢୣ୪ converges to the 9ൈ9 identity matrix. 3. Modified strain concentration tensor 3.1 Single inhomogeneity problem When the extended strain ሺ𝒁଴ሻ is applied to the infinite matrix ሺ𝑳଴ሻ having single inhomogeneity ሺ𝑳ଵሻ, the extended strain field within the inhomogeneity ሺ𝒁ଵሻ can be obtained by using the superposition principle. When the interfa cial damage is absent, the problem is decomposed into tw o homogeneous problems (Dunn and T aya, 1993a, c). The first problem is the homogeneous material ሺ𝑳଴ሻ under the applied extended strain ሺ𝒁଴ሻ, and its solution is given as a uniform strain field, equal to the appli ed extended strain, over the entire region. The second problem is the single inclusion problem havi ng an equivalent eigenstrain expressed as 𝒁ா௤ൌെ ሾሺ𝑳ଵെ𝑳 ଴ሻ:𝑺൅𝑳 ଴ሿିଵ:ሺ𝑳ଵെ𝑳 ଴ሻ:𝒁଴. (35) After superposing the solutions, the extended strain field with in the inhomogeneity is obtained as follows: 𝒁ଵൌ𝒁𝑰൅𝒁𝑰𝑰ൌ𝒁𝟎൅𝑺 :𝒁𝑬𝒒 ൌሾ𝑰൅𝑺 :𝑳଴ିଵ:ሺ𝑳ଵെ𝑳 ଴ሻሿିଵ:𝒁଴ (36) In the presence of the interfacial damage, an additional interf acial contribution due to the extended displacement should be inserted as the 3rd contribution (Fig. 5), 𝒁𝑰𝑰𝑰ൌെ 𝜞 : ሺ𝒁𝟏െ𝒁𝑬𝒒ሻ. (37) We note that the previ ous studies obtained the modified strain concentration tensor, simply by replacing the Eshelby tensor in Eq. (36) with the modified Eshe lby tensor, which results in a mismatching boundary condition at t he interface (Wang et al., 2 014a; Wang et al., 2014c). We have pointed out and corrected similar errors in the micromecha nics studies on elasticity (Lee et al., 2019a; Lee et al., 2019b) and heat conduction (Lee et a l., 2018). Finally, the extended strain field within the singl e inhomogeneity is expressed as 𝒁ଵൌ𝒁𝑰൅𝒁𝑰𝑰൅𝒁𝑰𝑰𝑰ൌ𝒁𝟎൅𝑺 :𝒁𝑬𝒒െ𝜞 : ሺ𝒁𝟏െ𝒁𝑬𝒒ሻ ൌ൫ 𝑰൅𝑺 :𝑳𝟎ି𝟏:ሺ𝑳𝟏െ𝑳 𝟎ሻ൅𝜞 :𝑳𝟎ି𝟏:𝑳𝟏൯: 𝒁଴ ൌ𝑨ெ:𝒁𝟎 (38) In the limit of zero interfacial damage, Eq. (38) converges to the strain concentration tensor for perfect interface in Eq. (36). 3.2 Numerical validation To validate Eq. (38), we compare the theoretical predictions wi th the FEA results by considering a composite having SiC reinforcement and PVDF matri x. The material properties of each phase are listed in Table 1, and the mesh configuration is the same as that the single inclusion problem (Fig. 3). The Mandel notation is used to carr y out tensor operations. As shown in Fig. 6, our predictions m atch well with the FEA result s for a wide range of interfacial damages. For a comparison, we also plot the modifie d concentration tensor obtained by replacing the Eshel by tensor in Eq. (36) with the m odified Eshelby tensor in Eq. (26), which deviates significantl y from the correct value in th e entire range of interfacial damages. In the absence of electric interfacial imperfection ( γ଴ൌ0), 𝐴ଷଷଷଷெ, which refers to the ratio between the 33 strain components inside the inhomogen eity and the far field 33 strains, converges to zero in the very large interfacial elasti c compliance limit ( α଴→∞ ) as shown in Fig. 6(a), because more strain energy is stored at the interfacial spring with increasing interfacial compliance. The same behavior occurs for o t h e r 𝐴௜௝௞௟ெ components concerning the relation between elastic strains (i.e. 𝑖, 𝑗, 𝑘, 𝑙 ൑ 3 ). Besides, 𝐴ସଷଷଷெ ( t h e r a t i o between the interior electric field along the 3rd axis and the applied 33 strains at the far field) remains finite in the limit of infinite interfacial compliance α଴→∞ (with γ଴ൌ0) as depicted in Fig. 6(b). Because t he strain field distribution in the piezoelectric PVDF matrix does not converge to zero but becomes identical to the strain f ield of a matrix having a spherical cavity, the finite surface charge accumulated at the piezoelectric matrix side of the interface produces a uniform electric field inside the SiC rein forcement, even in the α଴→∞ limit. When the interfacial electric imperfection γ଴ is present, 𝐴ସଷଷଷெ reduces and eventually approaches to zero in the γ଴→∞ limit. The same trend is observed for the all components of the interior electric field regardless of applied loading ty pe (mechanical strain or electric field), i.e., 𝐴ସ௝௄௟ெ decreases monotonically with increasing γ଴ and converges to zero as γ଴→∞ (Fig. 6(b)–(d)). 4. Effective properties of p iezoelectric composite 4.1 The Mori–Tanaka method and FEA condition Based on the Mori-Tanaka method, we predict the effective modul us of the SiC particle reinforced PVDF composite with corrected modified Eshe lby tensor and strain concentration tensor. The effective piezoelectric moduli tensor (Wang et al., 2014a; Wang et al., 2014c) of the composite havi ng spring-type interfacial dam age is expressed as 𝑳𝐄𝐟𝐟ൌሺ𝑐଴𝑳𝟎൅𝑐 ଵ𝑳𝟏:𝑨𝑴ሻ:ሺ𝑐଴𝑰൅𝑐 ଵ𝑨𝑴൅𝑐 ଵ𝑹: 𝑳 𝟏:𝑨𝑴ሻି𝟏 (39) where 𝑐଴ and 𝑐ଵ denote the volume frac tion of the matrix a nd particle, respect ively so 𝑐଴൅ 𝑐ଵൌ1. 𝑹 represents the intensity of the i nterfacial damage, which is e xpressed as 𝑅ெ௡௃௤ ൌ𝜔௃ 𝑎𝐼ெ௡௃௤ , 〈𝑹〉୑ୟ୬ୢୣ୪ൌ1 𝑎 ⎣⎢⎢⎢⎢⎢⎢⎢⎡𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 െ𝛾 െ𝛾 െ𝛾⎦⎥⎥⎥⎥⎥⎥⎥⎤ (40) Using the effective modulus, we can predict the extended strain /stress field within each phase. The applied extended strain ሺ𝒁ഥሻ and stress ሺ𝚺ഥሻ for a composite can be written in terms of the volume averaged extended strain 𝒁ഥ𝒊 and stress field 𝚺ഥ𝒊 within each phase as 𝒁ഥൌ𝑐 ଴𝒁ഥ𝟎൅𝑐 ଵ𝒁ഥ𝟏൅𝑐 ଵ𝑹: 𝑳 𝟏:𝒁ഥ𝟏 𝚺ഥൌ𝑐 ଴𝚺ഥ𝟎൅𝑐 ଵ𝚺ഥ𝟏 (41) We can rewrite the stres s equation in Eq. (41) as 𝑳𝐞𝐟𝐟𝒁ഥൌ𝑐 ଴𝑳𝟎:𝒁ഥ𝟎൅𝑐 ଵ𝑳𝟏:𝒁ഥ𝟏. (42) Hence, the strain and stress field within each phase are expres sed as Eq. (43) and (44), respectively. 𝒁ഥ𝟎ൌ1 𝑐଴ሺ𝑳𝟎െ𝑳 𝟏:ሺ𝑰൅𝑹 :𝑳 𝟏ሻି𝟏ሻି𝟏:൫ 𝑳𝐞𝐟𝐟െ𝑳 𝟏:ሺ𝑰൅𝑹 :𝑳 𝟏ሻି𝟏൯ 𝒁ഥ𝟏ൌ1 𝑐ଵሺ𝑳𝟏െ𝑳 𝟎:ሺ𝑰൅𝑹 :𝑳 𝟏ሻି𝟏ሻି𝟏:൫ 𝑳𝐞𝐟𝐟െ𝑳 𝟎൯ (43) 𝚺ഥ𝟎ൌ1 𝑐଴ሺ𝑴𝟎െ𝑴 𝟏െ𝑹 ሻି𝟏:൫ 𝑴𝐞𝐟𝐟െ𝑴 𝟏െ𝑹 ൯ 𝚺ഥ𝟏ൌ1 𝑐ଵሺ𝑴𝟏െ𝑴 𝟎൅𝑹 ሻି𝟏:൫ 𝑴𝐞𝐟𝐟െ𝑴 𝟎൯ (44) For the validation of the theor etical predictions, we calculate the effective modulus of the particle reinforced represent ative volume element (RVE) usi ng FEA. We construct RVEs having about 20 particles and calculate the average effective m odulus of 10 independent RVEs with a different particle distribution to obtain a statist ically meaningful value (Fig. 7). To obtain the effective piezoelectric modulus, we calculate the ratio between the volume averaged stress/electric displ acement field and the applied uni t strain/electric field. 4.2 Effective properties of compo site having one type of interf acial damages F i r s t , w e s h o w t h e e f f e c t o f e l a stic interfacial damage on the effective modulus, dielectric constant, and piezoel ectric constant of SiC-PVDF com posites with two different volume fractions in Fig. 8. The effective modulus predicted wit h our modified Eshelby tensor and strain concentration tensor matches well with the FEA resul ts up to 10% volume fraction, showing a relative error less than 2% (Fig. 8(a)) . In the small interfacial damage regime, the elastic stiffness of a composite with high volume fraction is h igher than that of a composite having low volume fraction because the composite is stiffened w ith additional SiC particles. However, as the interfacial damage increases, the stiffness of a composite with higher SiC volume fraction decreases compared with that of a composite wit h lower SiC volume fraction, because SiC effectively serves as a soft inhomogeneity when the load transfer between the SiC particle and the PVDF matrix is excessively weakened. As th e elastic interfacial damage approaches infinity, the predicted modulus converges to the por ous matrix result (whereas the prediction from the previous st udy converges to zero) (Fig. 8(f )). The effective piezoelectric constant ሺെ𝑒 ଷଷଷሻ also decreases as the elastic interfacial damage increases whe reas the effective dielectric constant was n o t sig n if ic an t l y af f ected b y the elastic interfacial damage (Fig. 8(b) and (c)). To understand the origin of decrease in th e effective piezoelectric constant, we plot the volume averaged strain of each phase in Fig. 8(d). The volume averaged strain within both matrix and particles decrease with the increase in the interfacial damage, because the interfacial spring e xtends more with inc reasing spring comp liance, which leads to reduced charge induction from the piezoel ectric matrix (Fig. 8(e)). In the infinitely large elastic interfacial damage without elect ric interfacial damage, the com posite has higher effective piezoelectric constant than that of the porous material because t h e S i C p a r t i c l e i s a l s o polarized owing to the finite in terior electric field (Fig. 8(g )). Second, we present the effective properties of the composite in the presence of the electric interface damage (Fig. 9). For a small electric interf acial damage, the effective dielectric constant κଷଷ increases with the SiC volume fraction because the dielectric constant of SiC is larger than that of the PVDF matrix. At a la rge electric interfacial damage, a composite with higher volume fraction of SiC particles has lo wer effective dielectric constant because the SiC particles are not polarized effectivel y. When the electric interfacial damage approaches infinity, the effective dielectric constant d ecreases below the dielectric constant of the porous matrix ( Fig. 9(a) and (f)). This is beca use while the infinite electric interfacial electric damage indi cates zero electric field withi n the SiC particle, the porous material with vacuum or air void s can have a finite electric fi eld within the void. This is different from the elasticity or heat conduction proble ms reported in the previous studies (Lee et al., 2019a; Lee et al., 2018) where th e effective property of a composite in the infinite interfa cial damage limit is equal to the effective property of a porous matrix. At infinite interfacial damage, the Mori–Tanaka method with the interface spring model converges to the results obtained with the material prope rties of a particle set to be zero. In elasticity and heat conduc tion, the elastic stiffness and thermal conductivity of the void (vacuum) are zero, so the two results (porous and infinite interfacial damage) become identical. However, when addressi ng the dielectric response, gi ven that the vacuum void has unit relative permittivity (i.e., 𝜅௜௝ൌ𝛿 ௜௝), the effective dielectric constant at infinite electric interfacial damage is lower t han that of a porous material. The effective stiffness is almost constant regardless of the el ectric interfacial damage and the effective piezoelectric constant ሺെ𝑒 ଷଷଷሻ decreases as the interfacial electric damage increases (Fig. 9(b) and (c)). When the electric field is appli ed to the composite, the electric field within the matrix decrease s with the interfacial electric compliance 𝛾 ( S e e F i g . 9 ( d ) and (e)). Therefore, less mechanical strain is induced from the p i e z o e l e c t r i c m a t r i x , w h i c h implies a decreased െ𝑒 ଷଷଷ, and owing to the mechanically “tie” condition at the interfac e(i.e. 𝛼ൌ0 ), the െ𝑒 ଷଷଷ at infinitely large 𝛾 is larger than the lower limit (Fig. 9(g)). 4.3 Effective properties of compo site having two types of inter facial damages. We also predict the effective properties of the composite havin g both types of interfacial damages. Fig. 10 shows that the effective stiffness mostly depends on the elastic interfacial damage, and it decrea ses with an increase in elasti c interfacial damage, as expected. Our prediction had the maximum error of 0.78% compare d with the FEA results. Similarly, the effective dielectric constant is affected mostly by the electric in terfacial damage. However, the effective െ𝑒 ଷଷଷ depends on both interfacial damages because it considers the coupling between mechanical and electrical responses. As the tw o interfacial damages simultaneously approach infinity, the effective piezoelectric a nd dielectric constants become lower than the porous results, whereas the effective stiffness converges to almost the same modulus as that with a porous material (Fig. 10(d)–(f)). Here, we use the normalized interfacial damage parameter ሺα଴ൌ𝛼 𝐶ଵଵଵଵ଴/𝑎, 𝛾 ଴ൌ𝛾 𝜅ଵଵ଴/𝑎ሻ; thus, the small size of the particle results in a large normalized interfacial damage. Henc e, the effective properties depend on the size of SiC particles at fixed 𝛼 and 𝛾 (Appendix C). Then, we obtain the figure of merit in 33 directions ሺ𝐹𝑂𝑀 ଷଷሻ by predicting the effective 𝑑ଷଷଷ and 𝜅ଷଷఙ in the strain-charge form constitutive equation, 𝜺ൌ𝑪ିଵ:𝝈൅𝒅்∙𝑬 𝑫ൌ𝒅 : 𝝈൅𝜿𝝈∙𝑬 . (45) We predict the effective െ𝑑 ଷଷଷ in terms of both electric and elastic interfacial damage (Fig. 1 1 ) . T h e e f f e c t i v e െ𝑑 ଷଷଷ decreases with an increase in the electrical damage, whereas i t increases as the elastic interfacial damage increases. As expre ssed in Eq. (37), when an external stress is applied, the stress field within the matrix increases with the elastic interfacial damage because the average stress field within the particles approaches zero (Fig. 11(b)). Using the 𝜅ଷଷఙ shown in Fig. 11, we predict the effective 𝐹𝑂𝑀 ଷଷ by using Eq. (42). 𝐹𝑂𝑀 ଷଷൌሺ𝑑ଷଷଷሻଶ 𝜅ଷଷఙ (46) As shown in Fig. 11, the effective 𝐹𝑂𝑀 ଷଷ increases with both interf acial damages. At a large interfacial damage, 𝐹𝑂𝑀 ଷଷ is even higher than that of the pure matrix. This is similar t o the porous material results studied by other researchers (Dunn and Taya, 1993b; Zhang et al., 2017). For a porous material, the effective 𝐹𝑂𝑀 ଷଷ increases with porosity because 𝜅ଷଷఙ decreases more rapidly than ሺ𝑑ଷଷଷሻଶ with respect to the porosity. Similar to a porous material, a composite with a piezoelectric matrix and elastic p articles has high 𝐹𝑂𝑀 ଷଷ w h e n the interfacial damage is large. 5. Conclusion We predict the effective material properties of a piezoelectric composite, considering the elastic and electric interfacial damage. We derive the modi fied Eshelby tensor after correcting the equations used in previous studies and validate it by comparison with the FEA results for a wide range of two interfacial damages. We also co rrectly predict the modified extended strain concentration tensor by decomposing the single inhomogeneity problem into three independent homogeneous problems. Using the modified Eshe lby tensor and extended strain concentration tensor, we predict the effective propertie s of the reinforced composite, and note that our model agrees well with the FEA results up to a 10% volume fraction. The effective moduli are well bounded between two limit cases, perf ect bonding and porous results, except the effective dielectric constant. We suggest a 2D map of effective properties with respect to the two inter facial damage parameters. Acknowledgements This research was supported by the Basic Science Program(NRF-20 19R1A2C4070690), the Creative Materials Discovery P rogram(NRF-2016M3D1A1900038), and the International Research & Development Program (NRF-2018K2A9A2A12000223) of the National Research Foundation of Korea(NRF) funded by the Ministry of Science and ICT of Korea. Appendix A. Non-uniform strain fiel d within the inclusion. The extended strain field within the inclusion in non-uniform w hen the elastic interfacial spring compliance in normal and tangential directio n are different. As shown in Fig. A1, the strain field of s pherical inclusion is non-uniform whereas the strain field when 𝛼ൌ𝛽 is uniform. For ellipsoidal case, the strain field is nonunifo rm even if the spring compliances are same. Appendix B. Violation of Fubini-Tonelli theorem In this section, we explain the violation of Fubini-Tonelli the orem by suggesting a simple double integral problem. When the integral function has singularity within the integral domain, the order of integral can not be exchanged. For example , for the function 𝑥ଶെ𝑦ଶ ሺ𝑥ଶ൅𝑦ଶሻଶൌെ𝜕ଶ 𝜕𝑥𝜕𝑦arctan ሺ𝑦/𝑥ሻ, (B.1) the results of two integrals which has different integral order are different, as shown in Eq.(B.2) නቆ න𝑥ଶെ𝑦ଶ ሺ𝑥ଶ൅𝑦ଶሻଶଵ ଴𝑑𝑦ቇ 𝑑𝑥ଵ ଴ൌെ𝜋 4 නቆ න𝑥ଶെ𝑦ଶ ሺ𝑥ଶ൅𝑦ଶሻଶଵ ଴𝑑𝑥ቇ 𝑑𝑦ଵ ଴ൌ𝜋 4 (B.2) Appendix C. Effect of particle size on effective properties In the main text, we investigate the effect of interfacial dama ge on the effective properties using normalized interfacial damage, 𝛼଴ൌ𝛼 𝐶ଵଵଵଵ଴/𝑎, 𝛾଴ൌ𝛾 𝜅ଵଵ଴/𝑎, where 𝑎 i s radius of the particle in the composite. If we predict the effe ctive properties using interfacial damage, not normalized parameter, the effective properties decr ease with the size of the particle (See Fig. C1). As the size of the particle decreases a t a fixed volume fraction, the total interface area increases, which results in large interfac ial damage, so the effective properties decrease. Table 1 Material constants of the considered materials. The constants are obtained from the previous studies (Odegard, 2004; Wang et al., 2014c). Property PZT-5A PVDF SiC 𝐶ଵଵଵଵ (GPa) 121 3.8 483.7 𝐶ଵଵଶଶ (GPa) 75.4 1.9 99.1 𝐶ଵଵଷଷ (GPa) 75.2 1.0 99.1 𝐶ଶଶଶଶ (GPa) 121 3.2 483.7 𝐶ଶଶଷଷ (GPa) 75.2 0.9 99.1 𝐶ଷଷଷଷ (GPa) 111 1.2 483.7 𝐶ଶଷଶଷ (GPa) 21.1 0.7 192.3 𝐶ଷଵଷଵ ( G P a ) 21.1 0.9 192.3 𝐶ଵଶଵଶ ( G P a ) 2.28 0.9 192.3 𝑒ଵଷଵ ( C / m2) 12.3 0.0 0.0 𝑒ଷଵଵ ( C / m2) -5.4 0.024 0.0 𝑒ଷଶଶ ( C / m2) -5.4 0.001 0.0 𝑒ଷଷଷ ( C / m2) 15.8 -0.027 0.0 𝜅ଵ/𝜅଴ 916 7.4 10.0 𝜅ଶ/𝜅଴ 916 9.3 10.0 𝜅ଷ/𝜅଴ 830 7.6 10.0 𝜅଴ൌ 8.854187817 ൈ 10ିଵଶ C2/Nm2 Figures and captions Fig. 1. Schematic of the single inclusion problem Fig. 2. (a) Schematic of the interface spring model. The right figure shows the deformed state having a finite (extended) displacement jump across the interfa ce. Schematic representation of (b) displacement and (c) poten tial field of the single inclu sion problem having interfacial damage. Fig. 3. Mesh configuration of the sin gle inclusion/inhomogeneity probl em. Fig. 4. (a) 𝑆ଵଵଵଵெ, (b) 𝑆ଶଶଵଵெ, and (c) 𝑆ସଷଵଵெ as a function of the norma lized elastic interfacial damage. (d) 𝑆ସଷଵଵெ with respect to normalized electric interfacial damage. We nor malized the interfacial damage by 𝛼଴ൌ𝐶 ଵଵଵଵ 𝛼/𝑎, 𝛾଴ൌ𝜅 ଵଵ𝛾/𝑎. We also present the results obtained from equations used in the previous study (Wang et al., 2014c). Fig. 5. Schematic of the single inhom ogeneity problem having interface spring. Fig. 6. (a) 𝐴ଷଷଷଷெ, (b) 𝐴ସଷଷଷெ, and (c) 𝐴ସଷସଷெ as a function of the normalized elastic interfacial damage. (d) 𝐴ସଷଷଷெ with respect to the normalized electric interfacial damage. We normalized the interfacial damage by 𝛼଴ൌ𝐶ଵଵଵଵ଴𝛼/𝑎, 𝛾଴ൌ𝜅ଵଵ଴𝛾/𝑎 w h e r e 𝐶ଵଵଵଵ଴ and 𝜅ଵଵ଴ a r e t h e elastic stiffness and dielectric constant of the matrix, respec tively. We also present the results obtained from equations used in t he previous study (Wang et al. , 2014c). Fig. 7. Mesh configuration of the indep endent RVEs at two volume fract ions. Fig. 8. (a) Effective elastic stiffness, (b) piezoelectric constant, a nd (c) dielectric constant of SiC–PVDF composite with respect to the normalized elastic inter facial damage. (d) V olume averaged strain ሺ𝜀̅௜ሻଷଷ within each phase under the applied strain ሺ𝜀̅ଷଷሻ. (e) FEA results of 𝜀ଷଷ field under the applied 𝜀ଷଷ଴ for two interfacial damages (sm all elastic interfacial damage of 0.038, large elastic interfacial damage of 3.8). (f) Effecti ve stiffness and (g) െ𝑒 ଷଷଷ with respect to the SiC particle vol ume fraction changing the elasti c interfacial damage. Fig. 9. (a) Effective dielectric constant, (b) piezoelectric constant, and (c) elastic stiffness of the SiC–PVDF composite with respect to the normalized electric interfacial damage. (d) V olume averaged electric field ሺ𝐸ത௜ሻଷ within each phase under the applied electric field ሺ𝐸തଷሻ. (e) FEA results of the 𝐸ଷ field under the applied 𝐸ଷ଴ for two interfacial damages (small electric interfacial damage of 0.0074, large electric interfaci al damage of 0.2059). (f) Effective stiffness and (g) െ𝑒 ଷଷଷ with respect to the SiC partic le volume fraction changing the electric interfacial damage. The upper and lower dashed lin es represent the perfect bonding and porous results ൫𝜅௣/𝜅଴ൌ1 ൯ , respectively. Fig. 10. (a) Effective stiffness, (b) piezoelectric constant, and (c) d ielectric constant for two interfacial damages under fixed volume fraction of 5%. (d) Effe ctive stiffness, (e) piezoelectric constant, and (f) dielectric constant when the tw o normalized interfacial damages are the same ሺ𝛼଴ൌ𝛾 ଴ሻ. Fig. 11. ( a ) N o r m a l i z e d e f f e c t i v e 𝐹𝑂𝑀 ଷଷ of the composite, where the 𝐹𝑂𝑀ଷଷ଴ i s t h e 𝐹𝑂𝑀 ଷଷ of the pure matrix (PVDF). (b) V olume averaged stress field wi thin each phase under the applied stress. (c) Effective െ𝑑 ଷଷଷ and (d) 𝜅ଷଷఙ,୉୤୤/κ଴ for two interfacial damages. The volume fraction is 5%. Fig. A1. (a) Strain field within the inclusion (PZT-5A) for different s hape and interfacial damage parameters. Non-uniform strain field within (b) sphere a nd (c) spheroid with aspect ratio of 2. The minor axis is used for the normalized interfaci al damage. Fig. C1. Effective (a) elastic stiffness and (c) dielectric constant of the SiC–PVDF composite with respect to the size of the particle at fixed interfacial d amages. The volume fraction of the particle is 5%.