## Abstract

The boundary effect of one inhomogeneity embedded in a semi-infinite solid at different depths has firstly been investigated using the fundamental solution for Mindlin's problem. Expanding the eigenstrain in a polynomial form and using the Eshelby's equivalent inclusion method, one can calculate the eigenstrain and thus obtain the elastic field. When the inhomogeneity is far from the boundary, the solution recovers Eshelby's solution. The method has been extended to a many-particle system in a semi-infinite solid, which is first demonstrated by the cases of two spheres. The comparison of the asymptotic form solution with the finite-element results shows the accuracy and capability of this method. The solution has been used to illustrate the boundary effects on its effective material behaviour of a semi-infinite simple cubic lattice particulate composite. The local field of a semi-infinite composite has been calculated at different volume fractions. A representative unit cell has been taken with different depths to the surface. The average stress and strain of the unit cell have been calculated under uniform loading conditions of normal or shear force on the surface, respectively. The effective elastic moduli of the unit cell not only depend on the material proportion, but also on its distance to the surface. The present model can be extended to other types of particle distribution and ellipsoidal particles.

## 1. Introduction

For one ellipsoidal inhomogeneity or particle embedded in an infinite solid under a uniform far field, Eshelby originally proposed the elastic solution through replacing the inhomogeneity with an inclusion with the same material properties as the rest solid but introducing an eigenstrain on the inclusion to meet its stress equivalence [1,2], which has been called the equivalent inclusion method (EIM). Eshelby's EIM has been widely used to predict the effective material properties of composite materials, such as homogenized elastic moduli [3–5], electroelastic moduli [6], thermal steady-state heat conductivity [7–9] and thermal expansion coefficient [10–13]. Recently, Eshelby's EIM is further applied to the investigation of the mechanical properties of nanomaterial [14,15], metallic polycrystals [16] and other advanced materials. EIM has become a foundation for micromechanical modelling of composites through the homogenization of the local mechanical fields for the effective material behaviour, in which the composite material size is assumed to be much larger than the particles' size [17,18]. Unfortunately, no infinitely large composite actually exists. Simplifications and assumptions have been commonly used in micromechanics-based modelling.

When a composite is uniformly fabricated, in a statistical sense at a larger length scale, the composite can be treated as a statistically homogeneous material, and the stress and strain at a material point can be evaluated by the averages of stress and strain on a representative volume element (RVE) [17–20]. Here an RVE in a continuum body is a material volume that statistically represents the neighbourhood of a material point. The microstructure can be periodic [21–23], random [24,25] or even functionally graded materials [12,26,27]. From the relation between averaged stress and strain, we can derive an effective mechanical constitutive law of the RVE [28]. Inevitably, the overall elastic properties depend on such factors as the volume fraction of the particles, the elastic moduli of the matrix and particles, and the microstructure of the composites.

Whether an RVE can provide an accurate prediction of effective material behaviour has been an interesting problem. Drugan & Wills [29] used the ensemble average to obtain the effective elastic properties and found a maximum error of 5% when the minimum RVE size is no less than twice the particle size for any volume fraction. Their formulation has been derived using Green's function in the infinite domain, which has been used by Eshelby [1]. However, because there commonly exist particles close to the boundary of a composite, Eshelby's assumption for one particle embedded in an infinite matrix cannot be exactly satisfied, and thus the boundary effect cannot be included in the modelling [17].

Numerical methods have been commonly used in stress analysis of an RVE or a unit cell containing a finite number of particles, which can represent the microstructure of the composite. In this way, the boundary effect can be considered through periodic boundary conditions or other treatments [30,31]. However, owing to the high computational cost, the material system is limited to a certain scale. For example, Gusev [32] used Monte Carlo realizations to generate the random microstructure of a unit cell including 8, 27 and 64 spherical particles and studied the effective moduli of the RVE under the periodic boundary condition. Kanit *et al.* [24] adopted the Voronoi mosaic model to generate the random microstructure, and used the finite-element method (FEM) to compute the effective elastic moduli and thermal conductivity considering three types of boundary conditions for different sizes of RVE. They found that a large number of realizations are needed when a small RVE is used and different boundary conditions may produce different results for stress and displacement, respectively.

In summary, when the RVE size is finite, the effective material properties not only depend on the microstructure and volume fraction of particles, but also change with the load types (force versus displacement) or boundary effects on the RVE. Given a large composite but a small RVE, when an RVE is selected at a different distance to the boundary, different average stress–strain relations may be obtained. Durgan & Wills [29] suggested ‘using the half-space Green's function instead of that for the whole space’ to analytically investigate the boundary effect.

This paper introduces Mindlin's [33] fundamental solution of ‘a concentrated force in the interior of a semi-infinite solid’ into the EIM and studies the boundary effect on the material behaviour of a composite. When a uniform normal or shear stress is applied on the surface of a half-space containing particles, the local stress field can be calculated by applying the fundamental solutions for Mindlin's problem [33–37] along with Eshelby's EIM [1].

Note that Rongved [34] has also developed a fundamental solution for the point force acting on the semi-infinite domain with fixed boundary. The Green function of a semi-infinite domain with a fixed boundary has been used with Eshelby's EIM to investigate the stress concentration of a microvoid embedded in an adhesive layer under the far-field load [38]. The mismatch of material constants is represented by the eigenstrain. Extending this method from solid mechanics to fluid mechanics, we successfully simulated the sedimentation process of multiple particles [39].

To demonstrate and quantify the boundary effect on material behaviour, this paper will first investigate one inhomogeneity embedded in a semi-infinite solid located at different depths. Expanding the eigenstrain in a polynomial form and using the EIM, we calculate the eigenstrain and thus derive the elastic field. Eshelby stated that the eigenstrain is uniformly distributed inside the inclusion when the RVE including one ellipsoidal particle is infinitely large [1,2]. This statement is adopted by most of micromechanics models. However, the distribution of eigenstrain is unknown and can be very complex for a bounded RVE. The present study makes an asymptotic analysis and reveals that quadratic expansion of eigenstrain can provide high accuracy even for particles that are close to the boundary. When the inhomogeneity is far from the boundary, the solution recovers Eshelby's existing solution. Then the formulation can be extended to a many-particle system in a semi-infinite solid, which is demonstrated by the cases of two spheres. The comparison of the asymptotic form solution with that of FEM shows the accuracy and capability of this method. Finally, the solution has been used to illustrate the boundary effects on effective material behaviour for a semi-infinite simple cubic lattice particulate composite under a uniform force on the boundary. The elastic field is periodic in the horizontal plane parallel to the boundary but changes in the depth direction. A representative unit cell, which includes one inhomogeneity of the simple lattice, has been selected to describe the average elastic fields. The average strain of the unit cell changes in the depth direction, while the average stress remains the same. The effective elastic moduli of the unit cell has been calculated under uniform loading conditions of normal or shear force on the surface. Although for simplicity the present method has only been used in periodic particle distribution, it can be extended to other microstructures with higher order computational costs [24,29].

Particle-reinforced composite materials have been widely used in engineering applications. The effective material properties are the key features for material design and quality control. The conventional micromechanics models state that the all the RVE have the same effective properties regardless of their positions. This study reveals that the effective elastic moduli not only depend on the material proportion, but also depend on the location of the unit cell. The boundary effect may soften the particle-reinforced composites and the effective material properties will change with the size of samples to be tested.

In summary, given a microstructure of a composite sample, the proposed method can calculate the local field in the sample with boundary effect and predict the effective material properties. This work also provides a promising approach to conducting virtual experiments to predict the material behaviour for the novel material design and development.

In what follows, §2 reviews fundamental solutions for a concentrated force in a semi-infinite domain [37] and introduces an explicit form Green's function for Mindlin's problem [33]. Using the Green function technique [25], we derive the elastic field for a semi-infinite solid containing inclusions. A polynomial form of the eigenstrain field has been adopted in the derivation. Section 3 introduces Eshelby's EIM to solve the inhomogeneity problem, in which the inhomogeneities are replaced by inclusions with an eigenstrain field on each inclusion. Based on the stress equivalence, the eigenstrains are calculated and the elastic field is obtained. Section 4 presents some numerical results for a semi-infinite domain containing one and two particles as well as a composite with periodically distributed particles in a simple lattice. The present solution is compared with high fidelity finite-element models. Asymptotic analysis of the polynomial form of the eigenstrain using uniform, linear and quadratic distributions demonstrates the convergence and accuracy of the proposed method. The effects of boundary and particle interaction on the local elastic field has been discussed. For the semi-infinite periodic composite with a simple lattice microstructure of particles, using unit cells at different depths from the boundary, we compute the effective elastic moduli of the unit cells and investigate the boundary effect. Finally, some concluding remarks are provided in §5.

## 2. Formulation of one inclusion in a semi-infinite domain

Figure 1 illustrates a semi-infinite domain *D* (*x*_{3}≥0) containing a subdomain *Ω* centred at **x**^{c} with a surface *S* (*x*_{3}=0). Following Mura's definition [17], the subdomain *Ω* is called an ‘inclusion’, which exhibits the same stiffness as the matrix but as an eigenstrain *Ω* in a homogeneous solid with different material properties from the rest (*D*−*Ω*) but on which an existing eigenstrain is unnecessary. Therefore, a particle in a matrix is an inhomogeneity instead of an inclusion. However, when an inhomogeneity is subjected to a body force or a far-field stress, the mismatch between inhomogeneity and the matrix can be simulated by an appropriately chosen eigenstrain on the subdomain so that the elastic field can be obtained through solving the inclusion problem [25,40]. Notice that the inclusion problems for a semi-infinite domain or two half-infinite domains have been studied [36,37,41,42]. However, the inhomogeneity problem for a semi-infinite domain has not been solved in the literature yet. Using the fundamental solution, this section will review and formulate the inclusion problem, which will be used for the inhomogeneity problem in §3.

The fundamental solution for a semi-infinite domain with a concentrated force can be traced back to Boussinesq's solution [43], in which a concentrated force on the surface of the domain is considered. Mindlin [33,44] provided the solution for a concentrated force in the interior of the semi-infinite domain with traction-free boundary. Using a similar procedure, Rongved [34] obtained the solution for a concentrated force in the interior of the half-space with fixed boundary. He also provided a solution for two joint semi-infinite solids with an interior force [35]. Yu & Sanday [36] derived the elastic fields of two joint semi-infinite solids in which frictionless and fully bonded interfaces are analysed separately. In their derivation, Galerkin vectors [33] due to infinitesimal inclusion are first presented, with which the displacement is derived. Then the displacement caused by a finite inclusion is obtained by direct integration over the inclusion *Ω*. Fully bonded and frictionless interfaces are considered separately and different kinds of eigenstrains are analysed. Introducing the imaging displacement, Walpole [37] derived the general solutions for the two joined half-spaces with bonded/smooth interfaces. The Green function of the two fully bonded semi-infinite domains including a concentrated force will be introduced below. Assuming the second half domain is rigid or of zero stiffness, one can obtain the fundamental solution for a semi-infinite domain with fixed or traction-free boundaries, respectively [33,34].

Consider an infinite two-phase elastic solid in figure 2, which consists of two different semi-infinite, homogeneous and isotropic halves, with a fully bonded interface. The Cartesian coordinate has been set up with the origin on the interface of *x*_{1}−*x*_{2} plane at *x*_{3}=0. A concentrated force **F** is applied in the interior of the upper half at **C**^{0} and *x*_{3}≥0 and *x*_{3}<0, respectively. Specifically, the shear modulus and Poisson's ratio are written as *μ*^{0} and *v*^{0} for *x*_{3}≥0 and *x*_{3}<0, respectively.

The image point of **x**^{′} is written as **F** and **x**^{′}, can be written by a mirror projection as

Based on Walpole's formulation [45], the displacement at **x** in the infinite two-phase solid can be rewritten in terms of Green's function as
*G*_{ij}(**x**,**x**^{′}) depends on the location of **x**. For *x*_{3}≥0,
*x*_{3}<0,
*D*=(3−4*v*^{0})/2*C*, *ψ*=|**x**−**x**^{′}|, *ϕ*=1/|**x**−**x**^{′}|,

When

Note that the above mathematical form of Mindlin's solution is more concise than the existing form solution although they are consistent with each other [33–37]. For

Eigenstrain means a non-mechanical strain, which was first introduced by Eshelby [1] and then established by Mura [17] as a formal term in micromechanics to describe local inelastic strain, such as thermal expansion, phase transformation, initial strains, plastic strains and misfit strains. Using this concept, Eshelby [1,2] proposed that the stress disturbance caused by the inhomogeneity could be simulated by an inclusion with the same material properties as the matrix but with an appropriately chosen eigenstrain. Using the stress equivalent condition, the inhomogeneity problem can be reduced to an inclusion problem.

The distribution of the eigenstrain ** ϵ***(

**x**) is typically continuous over

*Ω*, which can be written in terms of the Taylor expansion of

**x**−

**x**

^{c}with the reference point at the particle's centre, such as

The constitutive law over the semi-infinite domain can be rewritten as

The Green function technique has been widely used to solve the above type of problem. By using the Green function which satisfies the boundary conditions, the displacement can be written in terms of eigenstrain as follows:

Using Gauss's theorem and considering no eigenstrain at the boundary, equation (2.9) can be rewritten as

Defining *C*_{nmkl}=*C*_{mnkl} is considered to make *Γ*_{imn}=*Γ*_{inm}.

Using the kinematic equation, the strain is written as

Considering that the eigenstrain is a continuous tensor function over the inhomogeneity, one can write it in a polynomial form of the coordinate as follows:

For simplicity, only constant, linear and quadratic terms are considered in this paper, i.e.

Following Mura's book [17], one can integrate equation (2.13) as

## 3. Formulation of inhomogeneities in a semi-infinite domain

The EIM is first used to solve for the elastic field of one inhomogeneity embedded in a semi-infinite domain and is then generalized to the multiple-inhomogeneity case.

### (a) One inhomogeneity

In the last section, the semi-infinite domain *D* with stiffness **C** contains an inclusion *Ω* with an eigenstrain field *ϵ**(**x**). The elastic field has been solved in equation (2.17). In this section, the subdomain *Ω* in figure 1 is replaced by an inhomogeneity with a different stiffness **C**^{1}. A uniform stress *σ*^{0} is applied in the far field. Owing to the material mismatch, the local field in the neighbourhood of *Ω* will be disturbed. To solve this problem the concept of Eshelby's EIM [1,2] will be used in this section. The particle *Ω* is firstly filled with the same matrix material of **C** but an eigenstrain *ϵ**(**x**) is applied in *Ω* to make the stress *Ω* in the two cases equivalent. Then the boundary-value problems can be equivalently solved.

For one inhomogeneity in an infinite solid subject to a uniform far-field stress, the eigenstrain is uniform [1,2,17]. However, owing to the uniform boundary on the semi-infinite solid, the eigenstrain is no longer constant. The variation of ** ϵ***(

**x**) should be continuous over

*Ω*, which can be written in terms of polynomial of

**x-x**

^{c}with the origin at the centre of the ellipsoid, such as

Substituting equation (2.7) into equilibrium equation yields

Therefore, if the eigenstrain can be derived from the above equation, one can obtain the elastic field. Using Taylor's expansion, the disturbed strain inside the particle can be written as

When a uniform force is applied on the boundary of the semi-infinite domain, a uniform average stress *σ*^{0} will be caused. For the equivalent inclusion problem, a uniform strain *C*=*C*^{1}−*C*. Substituting equation (3.3) into equation (3.5) and comparing the coefficients of the constant, linear and quadratic terms at left and right sides, one can obtain the following equations system:

Solving this linear equation system, one can obtain the eigentrain terms of **E**^{0}, **E**^{1}and **E**^{2}, with which the displacement is written as

### (b) Multiple inhomogeneities

The above formulation can be extended from a single particle to multiple particles in a very straightforward way. When many particles are inserted in the semi-infinite domain, their interactions will significantly affect the elastic fields. Figure 3 illustrates the whole system with many particles. To formulate this problem, a global coordinate system has been set up in figure 3, and each particle, say *Ω*^{I}, is centred at *Ω*^{J}, reads

The displacement of any point in the domain is caused by the eigenstrain of all particles, which can be obtained by the superposition as follows:

Expanding the disturbed strain on particle *Ω*^{J}, substituting it into equation (3.9)) one can obtain the following equation system by comparing the coefficients for different orders of

Owing to symmetry, there are total 60 independent components of the eigenstrain for each particle when quadratic terms are considered. Therefore, equation (3.11) can be rewritten in the matrix notation,
*E*] whose length is 60*n*. [*A*]^{J} is the local coefficient matrix with dimensions 60×60*n* whose elements are the coefficients in equation (3.11). [*F*]^{J} is the local force vector with length 60, whose elements are the summation of the body force multiplying the corresponding coefficients.

For each particle, one can write one set of equations as the above. Rewriting all the equations in the matrix notation, assembling all the [*A*]^{J} into the global coefficient matrix and all the [*F*]^{J} into global force vector, one could obtain

Solving this global linear equation of system, one can obtain the eigenstrain terms of all particles. Substituting them into equation (3.10), one can also obtain the displacement field and thus the strain and stress fields.

In the numerical implementation of this algorithm, for one particle with uniform, linear and quadratic forms of eigenstrain distributions, only 6, 18 and 60 d.f. are, respectively, needed. Using quadratic form of eigenstrain as an example, to investigate the elastic fields of a sample with *N* particles, we only need to solve 60*N* unknowns in the linear equation system. In comparison, FEMs need to mesh the entire domain and the domains near the particles need very fine mesh to obtain accurate results. Generally, hundreds of thousands of degrees of freedom are needed for three-dimensional FEM simulation of a large domain with several particles. It is formidable to simulate a large number of spherical particles. Therefore, the computational cost of our approach is fairly low compared with the FEM modelling.

Note that the boundary-value problem is solved by a continuum mechanics based approach, and the particle interactions have been implicitly included. The pairwise interaction exists between any pair of particles and each particle is subject to the effects of all other particles, which generally decay with the particle–particle distance rapidly. Therefore, one can set a cut-off distance when *A*]_{60n×60n} a sparse matrix for a large-scale particle system.

In summary, the well-known Mindlin solution, reorganized in an elegant one, is used to study the elastic field for one particle in a semi-infinite domain. Although previous work has considered a concentrated force on a point or a prescribed strain in an inclusion, this is the first work to use it to study heterogeneous materials. The solution has been extended to many-particle systems for actual composite material samples. Because each particle is represented by an eigenstrain function, this analytical treatment of particles enables the simulation of hundreds of particles without the necessity of meshing the particles [39]. To our knowledge, this can be among or beyond the highest number of particles that FEM has simulated in the literature [32].

## 4. Results and discussion

Section 3 provides the formulation of the stress fields of multiple spherical particles embedded into a semi-infinite domain under the Neumann (stress) boundary condition. This section will investigate the accuracy of the present formulation for a semi-infinite domain containing (i) one particle, (ii) two particles, and (iii) many periodically distributed particles in a simple cubic lattice. The accuracy of the formulation will be verified with the FEM for two axisymmetric cases of one particle and two top-down particles embedded in a semi-infinite domain. The boundary effects on the local field and particle interactions are also discussed. Using equation (2.16), the eigenstrains on particles can be expanded in a polynomial form. Basically, higher accuracy is expected when higher order terms of eigenstrains are used. To compare the accuracy of the proposed EIM, constant, linear and quadratic forms of the eigenstrain distribution will be used. The asymptotic analysis demonstrates that the quadratic form of eigenstrain provides sufficient accuracy for the cases in this study. Finally, the boundary effects on the local field and effective elastic behaviour of a composite containing many particles, which are periodically distributed in a simple cubic lattice, are investigated under a uniaxial load and a simple shear load on the boundary, respectively.

In the rest of this section, FEMs will be used to test the accuracy of the present method. However, the present method actually has many unique features that FEM cannot compete with. Firstly, the present method has much higher computational efficiency than the FEM. To investigate the elastic field of *N* particles near the boundary, this method only needs 60*N* unknowns of the linear equation system for the quadratic eigenstrain components while the FEM needs to mesh the entire domain, and the domain near the particles needs fine mesh to obtain the results with high accuracy. For the spherical domain, much more elements are required to simulate the smooth surface unless B- or T-spline elements are used. Secondly, with the present analytic-approximated solution, one can apply asymptotic analysis to check the accuracy of the solution easily even with some singularity effect (e.g. penny-shape crack), whereas the FEM model may not catch the physics if the mesh is not constructed appropriately. Finally, the complexity order of the present method does not change with microstructure, but the FEM generally needs some treatments for specific microstructure configurations.

### (a) One particle

An aluminium particle with radius *a*=10 μm is embedded in a semi-infinite high-density polyethylene (HDPE) at (0,0,*h*), where *h*=11, 20 and 30 μm to the boundary. The Young moduli and Poisson ratios for aluminium and HDPE are *E*_{1}=69 GPa, *E*_{0}=0.8 GPa, *ν*_{1}=0.334 and *ν*_{0}=0.38, respectively [38,46,47]. The strain components corresponding to the uniform load are

Figure 4 illustrates the distribution of eigenstrain component *x*_{1}–*x*_{3} plane, when the aluminium particle is located at *h*=11 μm which means 1 μm spacing to the boundary of a particle with a diameter of 20 μm. A significant difference is observed when the eigenstrain is assumed to be uniformly (cf. figure 4*a*), linearly (cf. figure 4*b*) and quadratically distributed (cf. figure 4*c*). As Eshelby [1,2] pointed out, the eigenstrain inside an inhomogeneity should be uniform when the inhomogeneity is embedded into an infinite domain subject to a uniform far-field stress. However, owing to the boundary effect, obviously, the eigenstrain is not uniform any more. A polynomial form of eigenstrain may approach the exact solution when higher order terms are used. We will compare the above three solutions with the FEM result to verify their accuracy. The FEM models built with the commercial software ABAQUS are used to verify our analytical solution. An axisymmetric domain with 0.5×1 mm is used to represent the semi-infinite domain with a free boundary at *x*_{3}=0. In total, about 92 000 unstructured 6-node triangle elements are used to calculate the solutions.

Figure 5*a* illustrates the variation of stress *σ*_{33} along the axis *x*_{3} (*x*_{1}=*x*_{2}=0), when *h*=11 μm. One can observe a considerable difference in the results predicted by EIM with a constant eigenstrain from the FEM result. When linear terms are introduced, much better results are obtained. When quadratic terms are considered, the results almost overlap with the FEM results, which indicates that the quadratic form of eigenstrains are sufficient for the present elastic analysis when the spacing to the boundary is higher than 5% of a particle's size. When the particle is embedded at *h*=20 μm (figure 5*b*), the results obtained by uniformly and linearly distributed eigenstrain are much closer to the FEM results than those in the previous case because the boundary effect decays with *h*. With the distance to the boundary *h* increasing, the difference of the results from uniform, linear and quadratic distributed eigenstrain decreases. When *h* is large enough, the model reduces to one particle embedded into an infinite domain [1,2], in which the eigenstrain is uniform inside the inhomogeneity and higher order terms disappear.

In figure 5*a*, *σ*_{33} is equal to the far-field tension at the boundary *x*_{3}=0, and decreases with *x*_{3} to the bottom of the particle. Inside the particle, it increases with *x*_{3} to the top of the particle. Then a sudden change of the slope is observed and *σ*_{33} increases quickly to its maximum value, which is about *x*_{3} and finally converges to the far-field tension *b* illustrates the variation of *σ*_{33} with *x*_{3}, when the aluminium particle is embedded at *h*=20 μm, whose trend is quite different from the one for *h*=11 μm. Stress *σ*_{33} is equal to the far-field tension *x*_{3} and reaches the peak value near the bottom of the particle. Then it decreases with *x*_{3} to the bottom of the particle. Inside the particle, *σ*_{33} increases with *x*_{3}. One can observe here the change of *σ*_{33} inside the particle is much smaller than that for *h*=11 μm. Then it goes outside the particle and increases to the maximum value quickly near the top of the particle. Finally, it decreases with *x*_{3} and converges to the far-field load. In figure 5*b*, the stress distribution on the particle is not symmetric either, and the peak values near the bottom of the particle are smaller than that near the top. Actually, when the distance between the particle and boundary is large enough, the boundary effect is negligible so that the result will reduce to Eshelby's original problem for one particle embedded in an infinite domain, with which the curve will become absolutely symmetric and the stress inside the particle is uniform.

To investigate the shear deformation, a uniform shear load *h*=11 μm. Figure 6*a* illustrates the variation of the shear stress with *x*_{3}. The shear stress *σ*_{31} equals to the far-field shear stress *x*_{3} and reaches the inflection point at the bottom of the aluminium particle. It continuously increases inside the particle and reaches the maximum value at the top of the particle. It then decreases with *x*_{3} fast and reaches the minimum value, i.e. ∼0.75

The results from EIM with different eigenstrain forms are compared. Observations show that the stress obtained from the uniformly distributed eigenstrain is much lower than the ones calculated by linear and quadratic forms of eigenstrains. Because the shear stress distribution in the particle is fairly close to a linear trend, both linear and quadratic forms of eigenstrains provide a similar prediction. Therefore, the quadratic distributed eigenstrain can be considered as the convergent solution from the asymptotic analysis point of view.

To investigate the boundary effect, *σ*_{31} versus *x*_{3} curves for three cases of *h*=11, 20 and 30 μm are illustrated in figure 6*b*. The shear stresses at the surface are equal to the far-field load. When *h*=11 μm, the stress increases with *x*_{3} to the inflection point at the bottom of the particle; whereas for the cases of *h*=20 and 30 μm, the stress first decreases along *x*_{3} axis and then increases to the inflection point at the bottom of the particle. Inside the particle, *σ*_{31} increases with *x*_{3}. The variation of stress in the particle decreases for higher distance values of *h*/*a*. One can predict that the stress distribution in the particle will approach a constant when the distance between the particle and boundary is high enough, which makes the case reduce to the classical Eshelby problem of one particle embedded into an infinite domain. Note that the simple shear problem is not axisymmetric; FEM modelling has not been conducted for the high-computational cost because it is not the focus of this paper.

### (b) Two top-down particles

When two aluminium particles are embedded into the same semi-infinite HDPE domain, the problem becomes more complicated due to the coupling of the particle–boundary interaction and the particle–particle interaction [48]. For simplicity, we consider two relative positions of particles: top-down and side-by-side. For the top-down case, the two particles are located at (0,0,*h*) and (0,0,*h*+*s*), where *s* is the centre–centre distance of two particles. Obviously, *s* should be higher than 2*a* to avoid the particles' overlap. The same uniaxial load of

First, an extreme case is considered first with *h*=2.0*a* and *s*=2.5*a*. The results from EIM using uniform, linear and quadratic forms of eigenstrains are compared with the FEM results in figure 7*a*. FEM results are obtained by using an axisymmetric model of the commercial software ABAQUS, in which approximately 135 000 unstructured 6-node triangle elements are used. One can find that *σ*_{33} increases with *x*_{3} from boundary to specific peak value and then decreases a little bit when it reaches the bottom of the first particle. Then it increases again until it reaches the midpoint of two particles, i.e. (0,0,*h*+*s*/2). Then it decreases with *x*_{3} until the top of the second particle. After a slight increment, it decreases again and finally converges to the far-field tension load. The four curves are fairly different. EIM with uniform and linear distributed eigenstrain cannot address the stress distribution accurately. The maximum stress from uniform-distributed eigenstrain has a relative error of 30%. However, the EIM with the quadratic form of eigenstrain agrees with the FEM results very well, in which the difference of the maximum stress is less than 3%.

When the distance between two particles is longer, the eigenstrain is simpler owing to the decay of particle–particle interaction. Figure 7*b* illustrates the variation of *σ*_{33} with *x*_{3} when *s*=3.0*a*. One can see that the difference of the four curves is smaller although the EIM results for uniform eigenstrain are still fairly distinct from the other three. When quadratic terms are introduced in the eigenstrain, much better agreement with the FEM results is obtained, which highlights the advantage of using the quadratic form of eigenstrain.

Note that when two particles are closer (cf. figure 7*a*), the maximum stress locates in the middle of the two particles, whereas when their distance is longer (cf. figure 7*b*), the maximum stress locates close to the boundary of the two particles.

### (c) Two side-by-side particles

When the relative positions of two particles are side-by-side, the problem is no longer axisymmetric. Therefore, FEM modelling will not be conducted. The two aluminium particles are centred at (−*s*/2,0,*h*) and (*s*/2,0,*h*). Under the same far-field load of *a* illustrates the variation of *σ*_{11} with *x*_{1}, for *x*_{3}=*h*. In this case, *s*=2.5*a* while *h*=1.1*a*. One can observe that *σ*_{11} is negative, which indicates that the particle interaction produces compression between particles in the *x*_{1} direction. The curve is symmetric. For *σ*_{11}| is zero at far field, and then increases with *x*_{1} and reaches the local peak value near the left particle at about *σ*_{11}| decreases to a local minimum value, say *σ*_{11}| increases fast and reaches the maximum value *h* and *s* are fairly small, the effects of boundary and particle–particle interactions are very significant, so that the eigenstrain distribution is very complex. Therefore, the EIM with constant and linear terms cannot provide accurate results for this case.

Figure 8*b* illustrates the symmetric distribution of *σ*_{33} with *x*_{1}. For *σ*_{33} decrease from *x*_{1} inside the particle. A large drop is also observed at the other boundary of the particle, where *σ*_{33} jumps from approximately

To investigate the particle interactions and boundary effects, different configurations of *h*/*a* and *s*/*a* are considered. Figure 9*a* illustrates the stress *σ*_{33} distribution along the *x*_{1} in the case of *x*_{3}=*h* and *s*/*a*=2.5. Only the right half of the configuration is illustrated considering the symmetry. For the three cases of *h*=1.1*a*, *h*=2.0*a* and *h*=3.0*a*, the stress *σ*_{33} inside the particle caused by *h*. The difference of *σ*_{33} between *h*=2.0*a* and 3.0*a* is much less than that between *h*=1.1*a* and *h*=2.0*a*, which means that the boundary effect on the local field in the neighbourhood of particles is small for *h*≥2.0*a*. Moreover, *σ*_{33} outside the particle are almost the same for the three cases. Figure 9*b* illustrates the stress *σ*_{33} distribution for different centre–centre distances, i.e. *s*=2.5*a*, *s*=3.0*a*, *s*=4.0*a* while *h*=1.1*a*. Stress component *σ*_{33} is greater than zero at the midpoint of two particles, i.e., *x*_{1}=0. And the value increases with *s*. When *s*=3.0*a*, *σ*_{33} is close to the far-field applied stress. The peak value of *σ*_{33} appears at the left and the right boundary of the particle. When particle–particle distance is larger, the peak value of stress is higher. Moreover, the stress near the boundary of the particle is negative. The absolute value is higher when *s* is smaller.

Notice that the stress components at the middle point between the two particles are affected by the particle–particle distance and the distance to the boundary. Figure 10 illustrates stress components at the middle point between two particle centres varying with *s* for three cases of *h*=1.1*a*, *h*=2.0*a* and *h*=3.0*a*. The difference between the three cases is fairly small, which indicates that the boundary effect on the stress at the middle point of two particles is limited. However, the particle interaction produces significant effect. When the particle–particle distance is very close, say *s*=2*a*, although only a tension stress *σ*_{11} is near *σ*_{33} is near *s*, *σ*_{11} and *σ*_{33} will increase and finally approach 0 and *σ*_{33} is negative for *s*=2*a*−2.4*a* and then turns to positive for *s*>2.4*a*, which indicates the stress direction can change with the particle–particle distance.

### (d) Boundary effect on effective elastic behaviour of a simple cubic lattice composite

The formulation has been extended to a semi-infinite domain containing many particles, which can represent a general composite material. The boundary effect may change the local field in the composite considerably and lead to different effective material behaviour, which has not been studied in the literature yet. Therefore, this is the first piece of work to investigate the boundary effect on the effective material behaviour of RVE after this problem is pointed out [29].

To simplify the problem, we consider a simple cubic lattice composite illustrated in figure 11*a*. Because the microstructure is periodic, using one unit cell containing one particle, one can calculate the average stress and strain under a load on the boundary. The relation between the average stress and strain can represent the effective elasticity of the composite. In general cases, the periodic boundary condition is used to solve for effective elasticity of a periodic composite [18]. However, owing to the boundary effect, the periodic boundary condition in the depth direction is not valid any more. Therefore, the effective elasticity of the unit cells at different depths will be different. In the following, we number the unit cell starting with the traction-free boundary, i.e. the bottom of first unit cell at the boundary.

Given a simple lattice composite in figure 11*a*, an aluminium particle with radius *a*=10 μm is embedded in a semi-infinite HDPE. Young's moduli and Poisson's ratios for aluminium and HDPE are *E*_{1}=69 GPa, *E*_{0}=0.8 GPa, *ν*_{1}=0.334 and *ν*_{0}=0.38, respectively. The distance between the bottom and the centre of first layer particle is *h*=2*a*. The centre–centre distance of neighbouring particle is *s*=4*a*. Owing to the geometry of this microstructure, the elastic field in the plane parallel to the boundary surface will exhibit a periodic pattern for unit cells, so that the particles in the plane will share the same eigenstrain distributions. Therefore, the eigenstrain distribution on a particle will only change along the depth direction. From the above discussion, when the centre–centre distance is large, the particle interaction decays significantly. For simplicity, only 26 particles are taken into account to solve the eigenstrain, i.e. the cut-off of the effective zone is

For each layer, the equivalent stress formula is built as equation (3.11). Using the particle–particle interaction cut-off, the number of neighbouring particles considered in each formula is only 26, as shown in figure 11*b*. It is expected that the boundary effect will decay along the depth direction. Therefore, the eigenstrain will converge to a stable value, when the the distance of the unit cell to the boundary is far enough. Then we can assume the eigenstrain of *n* layer and *n*+1 layer are the same.

Assembling the equations of all the layers into a linear equation system, the eigenstrain can be obtained. With the eigenstrain, the local displacement and strain can be obtained from equations (2.17) and (2.9). Furthermore, using the homogenization technique, the average stress and strain can be calculated and therefore the effective elastic moduli can be derived. Uniform uniaxial tension and simple shear along the surface will be considered to compute the effective Young's modulus, Poisson's ratio and the shear modulus, respectively.

#### (i) Uniaxial tensile loading on the boundary

Consider a unit cell at a certain distance to the surface containing an inclusion *Ω* with an eigenstrain *D* with the free stress boundary condition. The stress and displacement field in the domain are continuous. The average stress can be written as

Similarly, the average strain can be written as
*η* is the volume fraction of the inhomogeneity. Recall that there is only one particle in one unit cell. Therefore, substituting equation (3.8) into equation (4.2), one can obtain
*a*_{I} is the radius of particle *I*.

When the composite material is under uniaxial tension, equations (4.1) and (4.2) can be rewritten as follows:

Therefore, the effective Young's modulus and Poisson's ratio of unit cell *I* for the present loading condition can be written as

Under the far-field tension *x*_{1}–*x*_{3} planes are illustrated in figure 12. One can observe that due to the boundary effect, the *x*_{1} direction is symmetric, but the symmetry about *x*_{3} axis does not exist due to the boundary effect. Noted that in all the unit cells, no eigenstrain exists outside the particle. In unit cell 1, *x*_{3} direction. The differences of *x*_{3} direction, but the difference is much smaller compared with that of unit cell 1. From unit cell 3 to unit cell 6, the boundary effect becomes smaller and smaller, and therefore the stress along the *x*_{3} axis becomes nearly symmetric. Noted that with the boundary effect decaying, not only the symmetry along *x*_{3} is gradually recovered, the distribution is uniform. This is because the linear and quadratic terms of eigenstrain become smaller, i.e. the uniform component plays a more dominant role.

The effective Young's modulus and Poisson's ratio are predicted by equation (4.5). Figure 13 illustrates the effective material properties of the composite materials. The effective Young's modulus of the first unit cell is much lower than those of other unit cells, which indicates that the boundary effect causes the lower stress in the particle of the first unit cell, which leads to a lower effective Young's modulus. One can conclude that the free surface will soften the composites with particulate reinforcements. Then the effective Young's modulus increases along the depth as it approaches a stable value after a distance of 5–6 unit cells, which means that the effect of the boundary disappears and the effective Young's modulus converges to the effective Young's modulus of a composite material with the same volume fraction but no boundary effect.

One can observe that the convergence value is only lower than the value predicted by the FEM homogenization [49]. It is higher than the dilute model, the self-consistent model and the Mori–Tanaka model [18]. The predicted effective Young's modulus by the dilute model and the self-consistent model are fairly close.

Figure 13*b* illustrates the effective Poisson ratio of the composite material. One can observe that the Poisson ratio of the first unit cell is smaller than the others. However, the second unit cell shows a higher Poisson's ratio. Finally, the effective Poisson ratio converges to a stable value, which is between the values of two pure materials. The converged value is higher than the one predicted by FEM homogenization [49], but lower than the Mori-Tanaka model [18], the dilute model and the self-consistent model. One may observe a trend among the five predictions that the model with a higher Young's modulus prediction provides a lower Poisson's ratio prediction.

#### (ii) Uniform simple shear loading on the boundary

For a composite material with simple lattice subject to a uniform shear loading

The effective shear modulus 〈*μ*〉 can be calculated as

Under a uniform shear loading *x*_{1}–*x*_{3} planes are illustrated in figure 14. The eigenstrain *x*_{1} axis. But the symmetry about the *x*_{3} axis does not exist due to the boundary effect. The eigenstrain *x*_{3} direction. From unit cell 2 to unit cell 6, the symmetry about *x*_{3} is gradually recovering. This indicates that the boundary effect is decaying with the increasing distance to the boundary, as for unit cell 6, *x*_{3} axis. Eigenstrain distributions in unit cell 5 and unit cell 6 are almost the same, which means that the eigenstrain is convergent.

Figure 15 illustrates the effective shear modulus of the unit cells. Similar to the effective Young's modulus, the effective shear modulus of the first unit cell is much lower than others, which again highlights that the boundary effect can soften the particle-reinforced composites. Then the effective shear modulus increases and finally converges to a stable value. The convergence value is between the prediction of the FEM homogenization [49] and the Mori–Tanaka model [18]. The predictions of the dilute model and the self-consistent model are also illustrated.

#### (iii) Parametric study

Given the size of a cubic unit cell, say 40μm each side, if the particle size becomes larger, the volume fraction of the particle will be much higher. It is expected that the effective Young's modulus and shear modulus will also increase with the particle size, as *E*_{1}>*E*_{0} and *G*_{1}>*G*_{0}. Figure 16*a* illustrates the variation of effective Young's modulus with unit cells for the composites with different aluminium particle sizes, say *a*=10 μm, 12.5 μm and 15 μm, respectively. As expected, the composite has larger particle size shows higher effective Young's modulus. Moreover, the effective Young's moduli of the three cases all increase along the *x*_{3} direction. Furthermore, the boundary effect of the case *a*=15 μm is much more significant than that of the other two, which causes higher variation of the effective Young's modulus versus unit cell number. Actually, when the centre of particle remains the same at the centre of unit cells but its size increases, the spacing between the particle to the free boundary will decrease, so that the boundary effect is amplified.

Figure 16*b* illustrates the variation of shear modulus with unit cell number along the depth direction for the composites with different particle sizes. Similar to Young's modulus, the composite with a larger particle size not only shows higher effective shear modulus, but also exhibits more significant variation in the depth direction owing to the boundary effect. Note that the effect of particle–particle interaction will also increase when particle size becomes larger. If the particle size is very large, to reach the equivalent accuracy, a longer cut-off is needed and more neighbouring particles should be considered.

Note that the effective material properties also vary with the volume fraction of the inhomogeneity. Figure 17 illustrates the variation of the effective moduli with the volume fraction of the aluminium for the unit cell at the boundary (unit cell 1) and far from the boundary (unit cell 8). The effective Young's modulus increases with the volume fraction for both the unit cell at the boundary or far from the boundary, as illustrated in figure 17*a*. However, the effective Young's modulus of the unit cell far from boundary increases much faster than that of the unit cell at the boundary. The similar trend has also been observed in figure 17*b* for the effective shear modulus changing with volume fractions. Although the same geometry and volume fraction of particle are in the unit cells, the effective elastic moduli of the unit cells is lower when it is close to the boundary owing to the boundary effect.

The effective Young's modulus, Poisson's ratio and shear modulus of the composite material are analysed above with the aid of the EIM. It reveals that the boundary effect can soften the particle-reinforced composites. Note that although only the example of the composites with a simple cubic lattice is analysed, this method can be applied to composites with other periodic microstructures, e.g. FCC (face centred cubic) or BCC (body centred cubic). Moreover, the model can be applied to the composite with different particle shapes, say ellipsoids, and different constitutive properties, say viscoelastic or elastoplastic materials.

## 5. Conclusion

Using the fundamental solution of a semi-infinite domain with a point force inside, this paper derives the elastic fields of a semi-infinite solid containing one, two and a simple cubic lattice of particles based on Eshelby's EIM. The particle–particle interaction and boundary effect are investigated. For the simple cubic lattice distributed particles, the effective Young's modulus, Poisson's ratio and shear modulus are calculated. The dependence of the effective elastic moduli on the distance to the boundary is demonstrated with aluminium particles in HDPE. The conclusions are summarized as follows.

Firstly, for the semi-infinite HDPE domain containing one aluminium particle, the elastic fields are calculated under a uniform far-field normal and shear stress. The stress inside the particle is much higher than that outside the particle, and reaches peak values at the boundary of the particle. The effect of the stress free boundary decreases with the distance to the boundary. The accuracy has been verified by FEM, and the results obtained from the quadratic form of eigenstrain show excellent agreement with the high fidelity FEM results.

Secondly, for the semi-infinite HDPE domain containing two aluminium particles, the stress fields are investigated. The particle–particle interaction and boundary effect both make significant contributions to the elastic fields. When the centre–centre distance of two particles is small, the stress reaches its peak value at the midpoint of two particles. However, when the centre–centre distance is long, the particle–particle interaction is so small that the two particles can be treated as ‘isolated’. Therefore, the peak value of the stresses appears at the boundary of the particles. When the two particles are side-by-side, the stress at the midpoint of two particles increases with the centre–centre distance. When the centre–centre distance is short, the midpoint is under high compression. But when the centre–centre distance is long, the stress at the midpoint converges to the far-field load. Moreover, the variation of the stress with centre–centre distance at the midpoint is almost independent of the boundary effect.

Thirdly, for the semi-infinite HDPE domain containing a simple cubic lattice of aluminium particles, the eigenstrains remain the same in the particles located at the same layer but vary in the depth direction. They converge to stable values when the particles are far from the boundary. The effective Young's modulus and shear modulus vary with the unit cells along the depth direction. The effective Young's modulus and shear modulus are lowest at the first unit cell and then increase along the depth direction, which means that the existence of a free boundary can soften the composite under the stress loading condition. The boundary effect decays fast and the effective Young's modulus and shear modulus quickly converge to the stable values in 5–6 unit cells.

The present formulation can be extended to simulate a large-scale particle system with various microstructures and generalized to ellipsoidal particles embedded in the semi-infinite domain. Similar to the general Eshelby problem, this formulation can be extended to different particle shapes, such as fibres, slits and penny-shape cracks. It can be used to study interactions between particles with different sizes and their interactions with boundary.

## Competing interests

We declare we have no competing interests.

## Funding

This work is sponsored by the National Science Foundation CMMI 0954717, CMMI 1301288 and CMMI 1301160, whose support is gratefully acknowledged. The authors also appreciate the support of the Henry Mitchell Weitzner Research Fund, which has been and will be used in their research of roofing materials for solar energy applications and technologies.

## Acknowledgements

Mr Charles Iselin, Mr Eduardo Lavocat and Ms Gisele Ribeiro have helped with the manuscript preparation.

- Received March 10, 2015.
- Accepted May 14, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.