Interphases and interfaces are essential microstructural components of heterogeneous solids that play important, sometimes dominant, roles in transport processes in composite media. Even though the concept of an interphase (a perfectly bonded finite-thickness layer) is somewhat more physical than that of an interface (a zero-thickness surface, across which the fields are discontinuous), both concepts represent idealized mechanical models introduced to describe complex mechanical processes at adjacent boundaries of dissimilar solids. Thus, formulations of the models that accurately and efficiently describe the effects of interphases and interfaces on thermomechanical behavior of micro-and nanostructured solids are of utmost importance in the mechanics of materials. In view of this, there exists a large body of literature in which the interphase and interface models have been proposed and extensively studied in the context of heat conduction, elasticity, and thermoelasticity; see comprehensive reviews in [1][2][3][4].
Interphase models are primarily used to model composites with coated reinforcements. The coatings are typically designed to mitigate the effects of stress concentrations and enhance the toughness of the composites. In such cases, the layer properties are known. Another application of interphase models is in the modeling of transition regions that appear as the result of damage, diffusion, or chemical reactions (such regions are characterized by reduced rather than increased stiffness). In the latter case, the properties of the layers are average values obtained by, for example, some homogenization procedure. For a recent comprehensive review and comparison of analytical and numerical models for the thermoelastic behavior of composites reinforced by coated spheres, see [5] and the references therein. In the framework of computational micromechanics, a periodic homogenization approach using the finite-element method has been applied ( [5][6][7], among others). The available analytical models [8][9][10][11] are limited to a single coated inhomogeneity and uniform far-field loading. It is noteworthy that most publications on the problem deal with an isotropic elastic interphase. We are aware of only a few papers in which the thermoelastic behavior of a particulate composite with an anisotropic interphase has been addressed; see [12][13][14] and the references therein.
It is also possible to model a thin interphase as a zero-thickness (also called imperfect) interface. For coated reinforcements with thin and ultra-thin coating layers, this model is mostly used to reduce computational cost and bypass associated problems (e.g. ill-conditioning). In the case of a transition zone, this model is used because of the transition zone's typically unknown material properties and small thickness, which are difficult to estimate. The available interface models can be divided into two groups-phenomenological and asymptoticbased.
Phenomenological models, used mostly for the description of transition zones, endow an interface with its own energetic structure and require additional data (interface constitutive laws, material parameters, jump conditions across the interface). For example, cohesive models imply continuity of the traction at the interface, but allow for a jump in the displacements, while so-called coherent elastic interfaces imply continuity of the displacements, but allow for a jump in the traction. More elaborately, the so-called general interface models allow for jumps in both displacement and traction fields. Development of phenomenological models started as early as in 1940s, see, for example, [15]; the history of development and a long list of references on the topic can be found in [3,4,[16][17][18][19].
Asymptotic-based interface models are derived analytically from the fully resolved interphase problem using various types of asymptotic analysis, for example, a Taylor series expansion [20] or a perturbation method [21]. As such, they do not require any additional assumptions or data. Typically, the asymptotic-based imperfect interface models are general, in the sense that both the displacement and normal traction vectors undergo jumps across the interface. However, it has been shown that, for limiting behavior, the asymptotic analysis-based and some phenomenological models concur. Important contributions to development of asymptotic-based models were made by Bövik [20], Hashin [22], and Benveniste [23], among others. The analysis of the literature relevant to particulate composites with imperfect interfaces leads to conclusions similar to those already listed for the interphase models, namely, (i) most analytical solutions deal with the case of a single particle and uniform far-field loading, while (ii) most interface models involve isotropic interfaces. Only a few papers (e.g. [23,24]) deal with curved anisotropic interfaces between two anisotropic media. In [25], the self-consistent and Mori-Tanaka homogenization schemes are extended to the case of an elastic particulate composite with a general imperfect interface by taking three (two in-plane and one orthogonal) elastic moduli of interfaces into account. Also, asymptotic interface models have been applied to continua with microstructure and multiphysics problems [26][27][28].
In the last 20 years, attention to the topic of interfaces and interphases has increased significantly in connection with developments in modern nanotechnologies and nanomaterials. For nanostructured materials, the interphase or interface effects on the local thermomechanical fields and macroscopic properties are even more substantial, owing to higher interface area-to-volume ratios than in traditional materials. In search of tools that can adequately model nanoscale phenomena, researchers turned their attention to the theory of material surfaces developed in the 1970s by Gurtin and Murdoch [29] and generalized in the 1990s by Steigmann and Ogden [30]. Both Gurtin-Murdoch (G-M) and Steigmann-Ogden (S-O) theories became very popular and were extensively used to study composite materials with nanosized reinforcements, see the reviews in [3,31,32]. In the context of particulate composites with spherical reinforcements, analytical solutions for a single spherical particle with the G-M interface were obtained and used to model the elastic fields [33][34][35][36][37] and effective properties [10,38,39] of particulate nanocomposites. Similar solutions for the S-O interfaces are reported in [40][41][42][43][44]. More advanced, finite-cluster [45] and representative unit cell [46] models of spherical particle composite with the G-M interface have been developed.
While theoretical interest in the G-M and S-O models remains high, enthusiasm about their applicability to nanostructured solids has recently subsided. There are several reasons for this. First, both theories are phenomenological, and, as such, require additional data, for example, on surface elastic properties, which are currently lacking or can even be inconsistent. For example, the results of molecular dynamics calculations revealed that the two models predicted dramatically different elastic moduli of nanostructures under bending and under tension [47]. For some combinations of material parameters and loading, the G-M model is inconsistent [48]. Second, these theories were proposed for free-surface problems and model the surface as a two-dimensional prestressed membrane or a shell of vanishing thickness that adheres to a three-dimensional bulk solid without slipping. In application to the interface problems, this means that the interface is treated as a coherent (elastic) interface, that is, the displacement vector is continuous across it while the traction undergoes a jump. However, the validity of such an assumption for nanomaterials has never been established. Third, these theories predict the elastic contact to be perfect in micro-and macro-heterogeneous materials, which is far from being always true.
The interfaces in real heterogeneous solids are, as a rule, incoherent, owing to inconsistency of the atomic lattices of contacting solids, dislocations, vacancies, and so on [49]; this gives a sound reason to consider the interface as a zone of reduced (rather than increased) stiffness. The material interface model allowing for the displacement discontinuity is considered in a few papers. Gurtin et al. [50] proposed a general theory of curved deformable solid interfaces in a polycrystalline solid. Another generalized continuum framework for modeling the elastic coherent and incoherent interfaces under general loading conditions was proposed in [51]. A common feature of these theories is an enlarged number (five and four, respectively) of interface elastic constants. This means that they cannot, after all, be derived in the isotropic elasticity framework.
The pertinent question is: What could be regarded as a "proper" model of an interface? In our opinion, the requirements for such a model must include its applicability at various structural levels or length scales (from the nano-to the macroscale) and an ability to catch essential features of a real interface (incoherency, size effect, surface stress, etc.). An important point is also the physical significance of the model parameters and the possibility of their theoretical or laboratory assessment. Among those known in the literature, the theory proposed in [50] seems to meet these requirements to the largest (although not full) extent. However, this advanced and promising model has not yet found application in the mechanics of materials.
The aim of this work is twofold and consists of (i) developing a model of an anisotropic imperfect interface in an elastic particulate composite with initial stress and (ii) applying the model to thermoelastic matrix type composites with incoherent interfaces. The paper is structured as follows. In Section 2, the thermoelastic problem for a particulate composite with a spherically anisotropic interphase and initial stress is formulated. In Section 3, the formal solution to this problem is derived for the case of a hydrostatic far-field load and the first-order accurate asymptotic procedure is applied to reduce the problem to that of an anisotropic imperfect interface. These solutions are used in the Maxwell-type estimates for the effective bulk modulus and the thermal expansion coefficient of the composite under study. In Section 4, the imperfect interface model [23] is generalized to account for the uniform eigenstress. The relevance of this model to the theory of material interfaces and its applicability in the nanomechanics context are discussed. In Section 5, a rigorous analytical solution for the unit cell model of a thermoelastic spherical particle composite with a transversely isotropic interphase and an anisotropic imperfect interface is obtained using the multipole expansion method. The accurate numerical data given in Section 6 reveal a profound effect of the interphase or interface anisotropy and eigenstress on the stress concentrations and effective thermoelastic properties of a spherical particle composite. In Section 7, we present a discussion of our results and conclusions. The background theory is provided in Appendices A to C.
Consider an elastic particulate composite comprising a homogeneous isotropic matrix solid and spherical inhomogeneities of equal radii R 1 . Each inhomogeneity consists of an isotropic core of radius R h = R 1 -h and an anisotropic interphase layer of thickness h. The composite medium is subjected to the uniform far-field stress σ far and the initial stress σ 0 . To be specific, we assume that the latter is the thermal stress owing to the uniform temperature change T and the difference in the coefficient of thermal expansion (CTE) of constituents.
In the linearly elastic solid, the small strain tensor ε = ε ij i i ⊗ i j relates the displacement vector u = u i i i (u = u (0) in the matrix, u = u (1) in the core inhomogeneity and u = u (c) in the interphase) as ε = (∇ ⊗u+∇ ⊗u T )/2. The stress tensor σ = σ ij i i ⊗ i j relates ε as σ = C : ε + σ 0 , where C is the fourth-rank elastic stiffness tensor. The matrix (i = 0) and core (i = 1) materials are isotropic, with the Poisson ratio ν = ν i and shear modulus µ = µ i . For them, the Duhamel-Neumann law is σ (i) = 2µ i ε (i) + λ i trε (i)β i T I, where I is the second-rank unit tensor and λ i = 2µ i ν i / (1 -2ν i ) is the Lamé constant. Also, β i = 3k i α i , where k i = (2µ i + 3λ i )/3 is the bulk modulus and α i is the CTE of the ith material. The thermal stress tensor is σ (i) 0 = -β i TI. The elastic equilibrium requires that div σ = div C : ∇u (i) = 0 (i = 0, 1).
The elastic moduli of the interphase layer possess spherical transverse isotropy. Introduced by Saint-Venant [52], spherical anisotropy implies spatial variation of the components of the tensor C ij = C ij (x) in a way that ∂C ij /∂α ≡ 0 (α = r, θ, ϕ). Hereinafter, two-index notation is adopted for the components of the elastic stiffness tensor C. In the spherical coordinate system Orθϕ, with the Or axis aligned with the anisotropy axis of the transversely isotropic material, Hooke's law reads
where Both the core-to-interphase and interphase-to-matrix interfaces are assumed to be perfect, which means that the displacement u and normal traction t = σ • n vectors are continuous across these interfaces:
Here, n = e r is the outward unit vector normal to the spherical surface, [[w]] = w + -w -denotes the jump of the w field across the interface, and the superscript "+" ("-") indicates the fields in the domains with the outward (inward) normal.
To make our presentation more clear, we start with the simple one-dimensional problem. Specifically, we consider an infinite solid with a single inhomogeneity loaded by the external hydrostatic pressure p and temperature step T. Spherical symmetry of the geometry and loading determines the spherical symmetry of the elastic fields. This implies that only the radial component of the displacement vector is nonzero and that all the fields are functions of the radial coordinate r. Analytical expressions for the displacement, strain, and stress fields in the matrix and core inhomogeneity for the spherically symmetrical problem are well-known (see e.g. [53]). The matrix fields are
where C is the equiaxial far-field strain and p = 3k 0 C is the corresponding far-field hydrostatic pressure. In equation (3.1), A is the unknown constant, whereas C is regarded as the known loading parameter. In the case of unconstrained thermal expansion C = α T, the total (elastic plus thermal) stress σ (0) rr = pβ 0 T vanishes at infinity. The elastic fields in the core inhomogeneity are
where D is the unknown constant.
The solution for the interphase layer is somewhat more involved. We take the radial displacement in the form
where E, F, and G are the constants to be found. It is noteworthy that the term Gr in equation ( 5) represents the particular solution to the non-homogeneous equilibrium equation aiming to counterbalance the body forces caused by the initial (thermal, in our case) stress in an anisotropic solid. This displacement generates the strains
and stresses
This stress field must obey the equilibrium equation div σ = 0, written in spherical basis as
Substitution of equation ( 7) into equation ( 8) yields
where υ = 2 (C 11 + C 12 -C 13 ) /C 33 . These results are consistent with the general theory provided in Appendix B.
The interface conditions of equation ( 2) are fulfilled by taking the appropriate constants A, D, E, and F. The first two of these (interphase-to-core continuity) are written, in our case, as
Substitution of the explicit expressions for u r and σ rr of equations ( 4) and ( 7) into equation ( 9) gives the following two linear equations:
Another two equations are obtained from the matrix-to-interphase continuity conditions of equation ( 2), written explicitly as
From equations ( 10) and ( 12), the constants A, D, E, and F are uniquely determined.
Now, we approximate the anisotropic interphase layer as the imperfect interface located in the midpoint R c = Rh/2 of the interphase layer (Figure 1). Specifically, our task is to find the appropriate functions F 1 and F 2 for the interface conditions
To this end, we apply the procedure developed by Benveniste [23], which provides the first-order approximation of the elastic fields owing to the interphase layer. Next, we outline this procedure in a slightly modified form.
As the first step, we expand u (c) r (R c ) -u (1) r (R c ) into the Taylor series in the vicinity of the point r = R h :
We are looking for the O(h) approximation, so all the O(h m ) (m ≥ 2) terms are neglected. By taking equations ( 4), (6), and ( 9) into account, we find that
It is instructive to compare equation (15) with expansion of the equation
This comparison says that the difference between hε rr (R c ) and hε rr (R h ) is of the order of O(h 2 ) and can be neglected in the first-order approximation. The same applies equally to all first-order terms of the Taylor series expansion that we consider next. It also follows from equations ( 4) and ( 7) that
By combining equations ( 15) and ( 16) with equation ( 9) and taking this remark into account, we find that
where
Analogous operations with equation ( 11) lead to
Now, summation of equations ( 17) and ( 19) yields the interface displacement jump condition in the form
where the interface average operator • is defined as w = (w (0) + w (1) )/2. Here, the constants c 1 , c 2 , and c 6 are the short notation for the c 1i , c 2i , and c 6i , respectively (i = 0, 1), defined by equation (18). Derivation of the second, normal traction jump condition follows a similar way. We expand σ (c) rr (R c )σ (1) rr (R c ) into the Taylor series in a vicinity of the point r = R h and employ the continuity conditions of equation ( 9) to get
Next, we express ∂σ rr /∂r in terms of σ rr and u r using equation (8). We have
Similarly,
Now, we combine equations ( 22) and (23) with equations ( 21) and ( 9) to obtain
where
Together with the expression
equation ( 24) results in
Equations ( 20) and ( 27) provide the explicit expression of F i in equation ( 13).
Solution to the single spherical inhomogeneity problem with the interface defined by equations ( 20) and ( 27) is straightforward. It follows from equations (3.1) and (4) that
Substitution of these expressions into equations (3.1) and ( 4) yields, after simple algebra, the following algebraic equations:
from which A and D are determined.
The obtained solutions enable an estimate of the effective bulk modulus k * and CTE α * of a spherical particle composite with transversely isotropic interphase and imperfect interface in the framework of the Maxwell homogenization scheme [54]. To be specific, we consider a composite with the volume fraction c of inhomogeneities imperfectly bonded to the matrix. In this case, the equivalent inhomogeneity is a sphere of radius R * (R c /R * ) 3 = c with the unknown effective properties k * and α * , perfectly bonded to the matrix. Continuity of the radial displacement and stress fields
The Maxwell scheme reads A * /C * = A/C, where A is found from equation (28). To find the effective bulk modulus k * of the composite with imperfect interface defined by equations ( 20) and (27), one has to find A from equation (28) for C = 1 and T = 0 and substitute the obtained value in place of A * in equation (29). The resulting formula is
Then, by solving equation (28) for C = 0 and T = 1, we get
The effective bulk modulus and CTE of the composite with interphase are also given by equations ( 30) and (31), provided that A is found from equations ( 10) and ( 12).
To complete this section, we note the following. Levin [55] has derived the formula relating the effective CTE to the effective bulk modulus of a two-phase heterogeneous solid with a perfect interface. In [10], the analogous formula is derived for the two-phase spherical particle composite with semi-imperfect (either with displacement or stress discontinuity) interface. Equations ( 28) and ( 29) can be regarded as an extension of Levin's formula to a composite with an anisotropic interphase or interface. To get it in explicit form, one has to derive the analytical expression of A from equations (10) and (12) or equation (28), respectively, and substitute it into equation (29). The resulting formula is cumbersome and we do not report it here.
In [23], the imperfect interface model is derived as the first-order approximation of a thin anisotropic interphase layer. Adding the initial stress to this model is straightforward and analogous to that done in the previous section for the particular problem. By analogy with [56], we rewrite the Benveniste model of a transversely isotropic layer with the thermal stress in compact form. Specifically, the displacement jump condition is
where
) c 1i and a 2i are defined by equation ( 18), and
Also, N = nn and P = I -N, where I is the second-rank unit tensor. The normal traction jump condition reads
where
and c 4i , c 5i , and c 7i are defined by equation (25). In equation (36), ε S = grad S u+ grad S u T /2 is the surface strain tensor, grad S is the surface gradient, and div S is the surface divergence, see [23] for their definition and properties. Equations ( 32) and ( 35) are reduced to equations (20) and (27) in the particular case of spherical symmetry.
The model of an incoherent interface of equations ( 32) and ( 35) involves eight independent parameters, all with clear physical meaning and conventional dimensionality. They are five elastic constants (C 11 , C 12 , C 13 , C 33 , and C 44 ), two CTEs (α 11 and α 33 ) and the length parameter h. We mention a few degenerate cases of this model for a very thin (h → 0) anisotropic interphase. Assuming that C 13 /C 33 = O (1), we have c 1 → 0 and hence
with the normal γ n = h/C 33 and tangential γ t = h/C 44 spring stiffness and two surface elastic constants,
Then, in the limit µ s , λ s → 0 (soft interface) we come to the thermoelastic spring layer model consistent with [10]:
] = 0. The opposite limit γ n , γ t → 0 (hard interface) yields the G-M model of coherent interface with zero bulk stress and interface stress σ 0 = -d t T:
In the case of an isotropic interface, these results are consistent with those reported in [10,57,58]. In the trivial case γ n = γ t = µ s = λ s = 0, we arrive at the conventional perfect interface conditions
Gurtin et al. [50] developed the deformation theory of a solid microstructure with curved incoherent interfaces.
To this end, the interface between amorphous phases with equal elastic constants was considered and the reference state (u = 0) was chosen such that the bulk stress at the interface is hydrostatic pressure p. If the interfacial free energy is independent of the relative displacement gradient, the force balance across the interface and the internal stress relation can be rewritten in our notation as
where f is the interface stress in the reference state and a k (k = 1, 2, . . . , 5) are the elastic moduli of the interface. The units of f , a 1 , and a 2 are [N/m], the units of p and a 3 are N/m 2 , and the units of a 4 and a 5 are N/m 3 . It is argued in [50] that the proposed theory applies to solids with a nanometer-scale microstructure. Importantly, equations ( 32) and (35) extend this theory to an incoherent interface between dissimilar elastic materials and provide a certain insight into the interface elastic moduli. For the sake of comparison, we rewrite these equations as
In the particular case c i0 = c i1 = c i (where the matrix and inclusion are made of the same material), the compared models coincide, provided we take
Recall that our approach is not confined to the thermal stress problem. The theory that we have developed is valid for eigenstresses of any kind, hydrostatic in the isotropic constituents, and hydrostatic plus hoop stress at the interface.
The following is also worth mentioning here. It is asserted in [50] that the moduli a 4 and a 5 are positive (negative) if the interface is more (less) compliant than the bulk material, while the opposite applies to a 1 and a 2 . Equations ( 18) and ( 25) are consistent with this assertion and specifically express a i in terms of the bulk and interface elastic constants. In particular, these equations give a simple explanation to the often debated issue of possible negative values of the interface elastic constants (being, in fact, a difference between two positive numbers).
To complete this discussion on the incoherent interface model, we note that h enters it as a free length (skin) parameter. Note that equation ( 39) is written in a such way that its right-hand-side terms describe perturbation owing to an interface. The small dimensionless number δ = h/R 1 is a parameter governing the interface contribution to the elastic fields and effective moduli. For a fixed h and the interface constants c i , this contribution can be interpreted as the inhomogeneity size effect, inherent in nanostructured solids.
Consider a periodic composite comprising a homogeneous matrix and spherical inhomogeneities composed of core particles of radii R h = R 1 -h and interphase layers of thickness h. To keep things simple, we assume that the inhomogeneities are arranged in a simple cubic (SC) array. The unit cell of this structure is a cube with side length a, containing a single inhomogeneity. The volume fraction of inhomogeneities is c = V 1 /V , where V = a 3 and V 1 = 4πR 3 1 /3 are the volumes of the unit cell and the inhomogeneity, respectively. Consideration of the many-particle representative unit cell model follows the same pattern [59].
The macroscopically uniform stress field in the composite bulk is assumed. This implies uniformity of the macroscopic strain ε and stress σ tensors, defined as
where S 0 is the outer surface of the unit cell. Importantly, this definition holds true for composites with interphases and imperfect interfaces. Periodicity of the composite microstructure results in quasi-periodicity of the displacement vector u (x
and periodicity of the corresponding strain and stress fields. In equation ( 41), E = ε is the uniform macroscopic strain tensor. Owing to the periodicity of the local fields, the unit cell serves as the representative volume element of the composite. The displacement vector u and the normal traction vector t = σ • n are continuous across the core-to-interphase and interphase-to-matrix interfaces (equation ( 2)). It is noteworthy that the periodic displacement boundary conditions of equation ( 41) automatically ensure the balance of angular momentum at the macroscale and guarantee symmetry of the macroscopic stress tensor defined by equation ( 40) (e.g., [60]).
To solve the model boundary value problem, we use the multipole expansion method. Its application to composites with isotropic constituents is discussed in detail elsewhere (e.g., [46,61]). Here, we outline the idea of the method and provide the necessary formulas. An appropriate formalism for the anisotropic interphase has been recently developed [56].
The periodicity conditions of equation ( 41) are fulfilled by taking the displacement vector in the form
where u far (x) = E • x is the linear displacement field corresponding to the uniform strain field ε = E and u dis is the spatially periodic perturbation field. In turn, u dis is expressed in terms of the periodic vector functions U (i) ts [46,61] as
Hereinafter, the short notation i,t,s = 3 i=1 ∞ t=0 t s=-t is used. The series expansion coefficients a (i)
ts are the complex constants to be found from the interface conditions, e.g., equation (2).
Obtaining an infinite set of the linear algebraic equations for the coefficients a (i) ts involves: (i) local series expansion of u (0) in the local spherical coordinate basis; (ii) substitution of the transformed u (0) , together with u (i) , into the interface conditions; and (iii) decomposition of the functional equalities using orthogonality of S (i) ts . The obtained infinite linear system is then appropriately truncated to t ≤ t max and solved numerically.
The local series expansion of u far (x) in a vicinity of inhomogeneity is
where the regular vector functions u
ts The final formula is [46,61]
where
ts are the irregular vector functions defined in Appendix A and
The series expansion coefficients η
ktls in equation ( 47) are the lattice sums providing periodicity of u dis . For their explicit expression, see [46,61]. The regular part of equation ( 46) is the disturbance field induced by all other inhomogeneities surrounding the selected one.
The displacement u (1) inside the core inhomogeneity is represented by the series over the regular vector functions u
where
ts are the unknown complex constants. The series expansion of the displacement vector u (c) in the spherical interphase layer is analogous to equation ( 46) but uses the set of vector functions v (i) ts and V (i) ts (x), (equation ( 70)):
where e
ts and f
ts are the series expansion coefficients. By analogy with equation ( 5), equation ( 49) additionally involves the linear term Gr aiming to counterbalance the body forces owing to thermal stress in the anisotropic interphase layer.
Obtaining the resolving linear system for the coefficients of these series expansions is straightforward. Let us consider the first equation of equation ( 2), namely, [[u]] = u (c)u (1) r=R h = 0. We substitute equations ( 46) and ( 49) into equation ( 2) and take the orthogonality property, equation ( 63), of S (j) ts into account to get the infinite set of linear algebraic equations for t > 0. We write it in the matrix-vector form as
where
ts } T . The matrices UM t , VM t , and VG t are defined by equations ( 66) and (71), respectively.
Decomposition of the second of equation ( 2), namely, [[t]] = (t (c)t (1) ) r=R h = 0 follows the same pattern and gives us another set of equations:
where the matrices TM t , WM t , and WG t are defined by equations ( 69) and (77), respectively. Fulfilling the matrix-to-coating (r = R 1 ) interface conditions is analogous and yields
where a kl = {a
ts } T , and c ts = {c
Equations ( 50) to ( 53) constitute a closed linear system, from which all the unknowns can be found with any desirable accuracy using the truncation method. Numerical solution of the truncated linear system enables accurate evaluation of the local displacement, strain, and stress fields at every point of the model composite. Also, the unit cell model perfectly matches the Rayleigh homogenization scheme for the effective stiffness of the composite.
To obtain the resolving linear system for a unit cell model with an anisotropic interface, we substitute equations ( 42) and ( 48) and corresponding normal traction vectors into the right-hand side of equations ( 32) and ( 35) and expand them in terms of S (j) ts . The derivation procedure is discussed in detail elsewhere [56]. Here, we give only the final formulas. The infinite set of linear equations resulting from equations ( 32) and ( 35) is as follows. For t ≥ 1, they are
For t = 0, they resemble equation ( 28) of Section 3:
For the explicit formulas of the matrices in equations ( 54) and (55), see [56] and Appendix C.
The macroscopic Duhamel-Neumann law reads
where the macroscopic strain ε and stress σ are defined by equation ( 40), C * is the effective elastic stiffness tensor, and α * is the tensor of the effective CTEs. Recall that in our model ε = E (equation ( 41)), whereas σ is expressed in terms of the induced elastic dipole moment t = t ij i i i j of inhomogeneity by the formula
For the spherical inhomogeneity,
,
,
where κ = 8πµ 0 (1ν 0 ) [62]. Equation ( 59) is valid for the inhomogeneities of arbitrary structure and interface bonding type and hence applies to both the interphase and interface problems under study. The Rayleigh scheme for elasticity (see e.g. [61]) is conveniently formulated in terms of the stiffness contribution tensor N of inhomogeneity related to the dipole moment t by the formula t = V 1 N : ε [63]. The components of the N tensor are found as N ijkl = t ij , where t is induced by the macroscopic strain
Provided that T = 0, the effective stiffness tensor is given by the simple exact formula
where C 0 is the elastic stiffness tensor of matrix solid. Equations ( 59) and ( 60) apply equally to composites with an anisotropic interphase or an anisotropic imperfect interface. Then, the effective CTEs are found from equation ( 57) as
where t is the dipole moment induced by the constrained thermal stress for ε = 0.
In this section, we give a few numerical examples to illustrate (i) the accuracy and validity limits of the firstorder approximation of the transversely isotropic thermoelastic interphase layer, (ii) an effect of the interphase anisotropy on the thermal stress and effective CTE of a composite, and (iii) its applicability as a model of the nanolevel incoherent interface. In what follows, we will conveniently use the technical elastic constants in parallel with the components of stiffness tensor. The Young moduli E i , shear moduli G ij , and Poisson ratios ν ij of transversely isotropic solid are related to C ij by
where = (C 11 + C 12 )C 33 -2(C 13 ) 2 and only five of these constants are independent.
The following fixed set of elastic constants is used in all numerical tests: ν 0 = 0.3, ν 1 = 0.2, µ 1 = 10µ 0 , ν 12 = ν 13 = 0.1, G 12 = G 13 = 5µ 0 , and E 1 = 4E 3 = 30µ 0 . By taking µ 0 = 1, we thereby assume that all the stress and elastic constants reported in this section are normalized by µ 0 . To minimize the number of model parameters and highlight the contribution of the interphase or interface, we assign α 0 = α 1 = 0. The variable parameters in our study are the interphase layer thickness h, interphase CTEs α 11 and α 33 , volume fraction c, and radius R 1 of inhomogeneities. To obtain the fully convergent series solution, the equations and unknowns with t ≤ t max were retained in the truncated linear system, where t max = 2 and t max = 25 for the single inhomogeneity and unit cell model, respectively. rr (R 1 ) as a function of interphase thickness h.
In the first numerical test, we estimate the accuracy of the first-order approximation of equations ( 20) and ( 27) in the case of spherically symmetric loading of an infinite solid with a single inhomogeneity. We consider the unconstrained thermal expansion C = α T, in this case, σ rr (r) → 0 with r → ∞. In Figure 2, the matrix thermal stress σ (0) rr (R 1 ) owing to the temperature change T = 1 is shown as a function of the interphase thickness normalized by the inhomogeneity radius: 0 ≤ h/R 1 ≤ 0.2. The solid curves represent the interphase model; in this case, σ (0) rr is calculated at the matrix side of the interface between the matrix and interphase layer. The dash-dotted curves represent its approximation by the interface model and the stress is evaluated in the matrix bulk r = R 1 > R c . Three CTE combinations are considered: α 11 = 0, α 33 = 2 (curve 1), α 11 = α 33 = 1 (curve 2), and α 11 = 2, α 33 = 0 (curve 3). It can be seen from the figure that (a) the thermal stress is greatly affected by the interphase CTE anisotropy and (b) the developed interface model works quite well. The maximum relative error of approximation reaches 4% for α 11 = 2, α 33 = 0, and h/R 1 = 0.2. For h/R 1 < 0.2 (recall that we assume the interphase layer to be thin), a closer agreement between the compared models is expected. Now, we proceed to the periodic composite containing a simple cubic array of spherical inhomogeneities. This composite possesses cubic symmetry of the effective elastic moduli and is isotropic with respect to the effective CTE. Again, we assume the unconstrained thermal expansion of the composite defined in this case by the condition σ = 0 or, alternatively, ε = α * T, see equation (57). An accurate solution to the unit call model problem involves the higher-order spherical harmonics. This narrows the validity range of equations (32) and (35) [56,59]. In Figure 3, the stress σ rr variation in the matrix along the arc r = R 1 , 0 ≤ θ ≤ π/2, ϕ = π/2 owing to the thermal load T = 1 of the composite with c = 0.45 and h/R 1 = 0.1 is shown. For this volume content and packing type, the minimal separation between the particle surfaces is as small as 0.1R 1 , which assumes significant local stress concentration. By analogy with Figure 2, the solid and dashdotted curves represent the interphase and interface model, respectively. Also, curves 1 to 3 correspond to the previously defined combinations of the interphase CTE. As seen, the compared models predict very close stress σ (0) rr values. However, for other components of the stress tensor, the agreement is not as good. In Figure 4, the stress σ θθ variation along the arc r = R 1 , 0 ≤ θ ≤ π/2, ϕ = π/2 is shown for h/R 1 = 0.02 (curve 1), h/R 1 = 0.05 (curve 2), and h/R 1 = 0.10 (curve 3). The first-order approximation is, expectedly, quite accurate for a thin (h = 0.02R 1 ) interphase but becomes gradually worse with the increase in h/R 1 , especially in the narrow gap between the inhomogeneities. Moreover, the first-order model of the interface that we consider becomes unphysical for h > 0.1R 1 in the high-filled particulate composite [56,59].
It has already been discussed that the dimensionless number δ = h/R 1 is a parameter governing the interface contribution to the elastic fields and effective moduli of the composite. For a fixed skin parameter h and the interface constants c i , this contribution can be regarded as the inhomogeneity size effect. Also, the surface tension σ 0 is conveniently modeled by the thermal stresses β ii T. In Figure 5, the stress σ (0) rr at the pole point r = R 1 , θ = 0 is shown as a function of R 1 /h for c = 0.15, 0.30, 0.40, and 0. on the volume fraction and size of inhomogeneities, as well as on the other structure parameters of the composite and deserves a comprehensive study. This task, however, is beyond the scope of this paper.
In Table 1, the effective CTE α * of a periodic composite with a simple cubic array of spherical particles predicted by the two models is given for h/R 1 = 0.1, α 11 = 2, α 33 = 0, 0.1 ≤ c ≤ 0.5. The data in the two first columns are obtained using the Maxwell scheme (equation (31)). The next two columns give the data corresponding to the accurate multipole expansion solution (equation ( 61)). The data shown in the last column are found from numerical, finite-element analysis of the unit cell model of the composite with a transversely isotropic interphase. Our accurate data for the composite with the interphases practically coincide with the finite-element solution. The interface model overestimates α * , but the deviation is small; even for c = 0.5 (which is close to dense packing), it is of order 0.2% . The Maxwell scheme also provides a reasonably good estimate, with a relative error of 5% for c = 0.5.
In Table 2, the effective CTE of a periodic composite with an imperfect interface is given as a function of h/R 1 , c, and the interface anisotropy type. As expected, α * is nearly proportional to the skin parameter h and volume fraction of inhomogeneities c. For fixed values of h and c, α * varies about three times for the combinations of α 11 and α 33 that we consider. In Figure 6, the effective CTE of the periodic composite with anisotropic imperfect interface is shown as a function of R 1 /h. By analogy with Figure 5, these data can be interpreted in the context of nanomechanics as the inhomogeneity size-dependent macroscopic elastic bulk strain of nanostructured solid owing to interface stress.
The findings of this study can be summarized as follows.
The model of an anisotropic imperfect interface in a heterogeneous solid with an initial stress is derived as the first-order approximation of a thin transversely isotropic interphase layer. The developed model assumes discontinuity of both the displacement and normal traction vector fields across the interface and involves eight independent parameters. This set involves five elastic constants, two eigenstress-related parameters (CTEs, in the case of thermal stress), and one length parameter. This ensures applicability of the model at various structural levels and flexibility in modeling the interfaces in composites under the presence of interface residual stress and discontinuity of the displacement and stress fields. The proposed theory covers most of the known interface models. In application to polycrystalline solids, our model is consistent with the theory of curved deformable interfaces [50] and provides a certain insight into the interface elastic moduli. It is believed that taking the incoherency of the interface into account makes the model more realistic and thus increases the reliability of predicting the elastic properties of nanostructured solids. The developed model is not limited to the spherical geometry considered here and can be applied to a variety of heterogeneous solids with incoherent interfaces. Such an extension might require computational efforts that could be a subject of future investigations. Another promising direction of future work is an extension of the developed approach to composites with a graded interphase.
To illustrate the derivation procedure and essential features of the anisotropic imperfect interface model, the spherically symmetrical problem for an infinite elastic solid with a single inhomogeneity is considered first. The analytical expressions for the displacement, strain, and stress fields in the matrix, anisotropic interphase, and core inhomogeneity are derived. Retaining the O(h) terms in the asymptotic expansion of the solution yields the displacement and stress jump conditions. The developed model is applied to estimate the effective bulk modulus and CTE of a spherical particle composite with the transversely isotropic interphase and imperfect interface in the framework of a Maxwell homogenization scheme. These results can be regarded as an extension of Levin's formula to composites with an anisotropic interphase or an imperfect interface.
The rigorous solution for the unit cell model of a spherical particle composite with a transversely isotropic interphase and thermal stress has been obtained using the multipole expansion method. Accurate fulfillment of the matrix-to-interphase and interphase-to-core inhomogeneity contact conditions reduces the model boundary value problem of thermoelasticity to the linear algebraic system for multipole strengths and provides a highly efficient algorithm for numerical study. This solution is readily incorporated in the many-particle representative unit cell model and thus enables consideration of a composite comprising spherical particles of diverse sizes and properties, with adequate account of their arrangement, interactions, and anisotropy of interphases. The analogous approach is used to find a rigorous solution for the unit cell model of a spherical particle composite with the newly developed model of an anisotropic incoherent interface. Numerical comparison of the interphase and interface models supports our theory and provides an estimate of its accuracy and applicability limits in a multiparticle environment. Both models are applied in the framework of a Rayleigh homogenization scheme for evaluation of the effective CTE of a spherical particle composite with imperfect interfaces.
The reported accurate numerical data reveal a significant effect of the interphase or interface anisotropy and eigenstress on the stress concentration and effective thermoelastic properties of particulate composite. This is in sharp contrast to coherent interface models that predict the size effect to be weak and are confined to composites with soft nanoinclusions and nanoporous solids. The extended set of material parameters used in the proposed interface model may significantly improve the reliability or accuracy of the micromechanical simulation of heterogeneous media.
ts = e r χ s t
For t < 0, they are defined as S
where the over bar indicates a complex conjugate, α
A complete set of the vector functions u
ts obeying Lamé equation is as follows [61]. The regular (bounded everywhere but for r → ∞) functions u
1s , which represent translation and rotation, respectively, of a rigid solid, are excluded from consideration.
The roots υ i and m i are arranged in such a way that 0 < m 1 < m 2 < m 3 for the regular functions and m 3 < m 2 < m 1 < 0 for the irregular ones. Given m i ,
where WM t (r) = r m r1 -1 p r1t 0 r m r3 -1 p r3t 0 r mr 2 -1 (m r2 -1) 0 r m r1 -1 q r1t 0 r m r3 -1 q r3t Again, the WG t matrix for t(V (i) ts ) has the same form as WM t , with m i and k i for the irregular functions. For more details, see [56].
The matrices entering equation ( 54) are defined as [56] UG t (R, ν 0 ) = UG t (R, ν 0 ) -FG t (R, ν 0 ), UM t (R, ν 0 ) = UM t (R, ν 0 ) -FM t (R, ν 0 ), UM t (R, ν 1 ) = UM t (R, ν 1 ) + FM t (R, ν 1 ), where FG t (R, ν 0 ) = h 2 [a 10 FG t (R, ν 0 ) + 2µ 0 C 0 TG t (R, ν 0 )] ,
Here,
The matrices entering equation ( 55) are TG t (R, ν 0 ) = TG t (R, ν 0 ) -TG t (R, ν 0 ),