PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
 
Proc Math Phys Eng Sci. 2015 July 8; 471(2179): 20150174.
PMCID: PMC4528662

Boundary effect on the elastic field of a semi-infinite solid containing inhomogeneities

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.

Keywords: elastic boundary-value problem, Mindlin's problem, Green's function, elastic materials, Eshelby's equivalent inclusion method, unit cell

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 [35], electroelastic moduli [6], thermal steady-state heat conductivity [79] and thermal expansion coefficient [1013]. 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) [1720]. 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 [2123], 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 [3337] 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 (x3≥0) containing a subdomain Ω centred at xc with a surface S (x3=0). Following Mura's definition [17], the subdomain Ω is called an ‘inclusion’, which exhibits the same stiffness as the matrix but as an eigenstrain ϵij(x) applied as a fictitious non-mechanical strain. On the other hand, an ‘inhomogeneity’ is a subdomain Ω 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.

Figure 1.
A semi-infinite domain D (x3≥0) containing a subdomain Ω with a surface S. (Online version in colour.)

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 x1x2 plane at x3=0. A concentrated force F is applied in the interior of the upper half at x=(x1,x2,x3) with x3>0. The stiffness tensors for the two phases of materials are denoted by C0 and C¯0 for x3≥0 and x3<0, respectively. Specifically, the shear modulus and Poisson's ratio are written as μ0 and v0 for x3≥0 and μ¯0 and v¯0 for x3<0, respectively.

Figure 2.
Two joined semi-infinite domains (x3≥0) containing a concentrated force F at x with the image force F¯ at x¯.

The image point of x is written as x¯ with the coordinate (x1,x2,x3). The corresponding image force is F¯=(F1,F2,F3). Notice that the image of a vector, say F and x, can be written by a mirror projection as

(¯)=Q()withQ=100010001,
2.1

where (•) and (¯) mean any vector in the original and mirror projected coordinate, such as F¯i=QjiFj and Fi=QijF¯j.

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

ui(x)=Gij(x,x)Fj(x),
2.2

where the form of Gij(x,x) depends on the location of x. For x3≥0,

4πμ0Gij(x,x)=ϕδij14(1v0)ψ,ij+Aϕ¯δijDQimQjnψ¯,mnCx3[Qjmψ¯,im3+4(1v0)δj3ϕ¯,i+2(12v0)δi3Qjmϕ¯,mQjmx3ϕ¯,im]+B(δi3δjkδikδj3)[ln(x3+x3+ψ¯)],k(G+B)Qjm[(x3+x3)ln(x3+x3+ψ¯)ψ¯],im
2.3

and for x3<0,

4πμ0Gij(x,x)=AϕδijDψ,ijx3F[ln(x3+x3+ψ)],ij+B(δi3δjkδikδj3)[ln(x3+x3+ψ)],k(G+B)Qim[(x3+x3)ln(x3+x3+ψ)ψ],jm,
2.4

where A=(μ0μ¯0)/(μ0+μ¯0), B=2μ0(12v0)(μ0μ¯0)/(μ0+μ¯0)[μ0+μ¯0(34v0)], C=(μ0μ¯0)/2(1v0)(μ0+(34v0)μ¯0), D=(3−4v0)/2C, F=2μ0[μ0(12v¯0)μ¯0(12v0)]/[μ0+μ¯0(34v0)][μ¯0+μ0(34v¯0)], G=μ0[μ¯0(12v¯0)(34v0)μ0(12v0)(34v¯0)]/[μ0+μ¯0(34v0)][μ¯0+μ0(34v¯0)], and ψ=|xx|, ϕ=1/|xx|, ψ¯=|xx¯|, ϕ¯=1/|xx¯|.

When μ¯0=0, the above Green's function is reduced to Mindlin's solution [33]

4πμ0Gij(x,x)=ϕδij14(1v0)ψ,ij+ϕ¯δij(34v0)4(1v0)QimQjnψ¯,mnx32(1v0)×[Qjmψ¯,im3+4(1v0)δj3ϕ¯,i+2(12v0)δi3Qjmϕ¯,mQjmx3ϕ¯,im]+2(12v0)(δi3δjkδikδj3)[ln(x3+x3+ψ¯)],k(12v0)Qjm[(x3+x3)ln(x3+x3+ψ¯)ψ¯],im.
2.5

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 [3337]. For μ¯0, the above Green's function is reduced to Rongved's solution [35,38].

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 xxc with the reference point at the particle's centre, such as

ϵij(x)=ϵij0+ϵijk1(xkxkc)+ϵijkl2(xkxkc)(xlxlct)+,xΩ.
2.6

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

σij=Cijkl(uk,lϵkl),
2.7

where ϵkl=0forxDΩ. The substitution of the above stress into the equilibrium equation yields

Cijkluk,li=Cijklϵkl,i.
2.8

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:

ui(x)=D[Gij(x,x)Cmjklϵkl(x)xm]dx.
2.9

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

ui(x)=ΩGijxmCmjklϵkldx.
2.10

Defining Γijmu=Gijxm and substituting it into equation (2.10), one can obtain

ui(x)=ΩΓijmCmjklϵkldx,
2.11

where

Γimn=12(Γimnu+Γinmu),
2.12

in which Cnmkl=Cmnkl is considered to make Γimn=Γinm.

Using the kinematic equation, the strain is written as

ϵij=12(ui,j+uj,i)=12Ω(Γimn,j+Γjmn,i)Cnmklϵkldx=ΩΓijmneCmnklϵkldx,
2.13

where

Γijmne=14(Γimn,ju+Γjmn,iu+Γinm,ju+Γjnm,iu),
2.14

which satisfies both minor and major symmetries of fourth-rank tensors.

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:

ϵij(x)=Eij0+Eijk1(xkxkc)+Eijkl2(xkxkc)(xlxlc),xΩ0,xDΩ.
2.15

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

ϵmn=Emn0+Emnp1(xpxpc)+Emnst2(xsxsc)(xtxtc).
2.16

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

ϵij=ΩΓijkleCklmn[Emn0+Emnp1(xpxpc)+Emnst2(xsxsc)(xtxtc)]dx=14πμ0Dijkl0CklmnEmn0+Dpijkl1CklmnEmnp1+Dstijkl2CklmnEmnst2,
2.17

where Dijkl0, Dpijkl1 and Dstijkl2 are provided in the electronic supplementary material.

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 C1. 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-xc with the origin at the centre of the ellipsoid, such as

ϵij(x)=Eij0+Eijk1(xkxkc)+Eijkl2(xkxkc)(xlxlc)+,xΩ,
3.1

where the eigenstrain shares the same physical meaning of polarization strain as Eshelby's problem.

Substituting equation (2.7) into equilibrium equation yields

Cijkl(uk,liϵkl,i)=Cijkl1uk,li.
3.2

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

ϵij(x)=14πμ0[Dijkl0(xc)CklmnEmn0+Dpijkl1(xc)CklmnEmnp1+Dstijkl2(xc)CklmnEmnst2]14πμ0Dijkl,r0(xc)CklmnEmn0+Dpijkl,r1(xc)CklmnEmnp1+Dstijkl,r2(xc)CklmnEmnst2(xrxrc)18πμ0Dijkl,yz0(xc)CklmnEmn0+Dpijkl,yz1(xc)CklmnEmnp1+Dstijkl,yz2(xc)CklmnEmnst2(xyxyc)(xzxzc).
3.3

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 ϵij0=Cijkl1σkl0 will be induced, the equivalent inclusion condition is written as

Cijkl1(ϵkl0+ϵkl)=Cijkl(ϵkl0+ϵklϵkl),
3.4

which can be rewritten as

ΔC:(ϵ0+ϵ)=C:ϵ,
3.5

where ΔC=C1C. 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:

[ΔCijklDklst(xc)Cstmn4πμ0Cijmn]Emn0+ΔCijklDpklst(xc)CstmnEmnp1+ΔCijklDpqklst(xc)CstmnEmnpq2=4πμ0ΔCijklϵkl0[ΔCijklDpklst,r(xc)Cstmn4πμ0Cijmnδpr]Emnp1+ΔCijklDklst,r(xc)CstmnEmn0+ΔCijklDpqklst,r(xc)CstmnEmnpq2=0andΔCijklDklst,yz(xc)CstmnEmn0+ΔCijklDpklst,yz(xc)CstmnEmnp1+[ΔCijklDpqklst,yz(xc)Cstmn8πμ0Cijmnδpyδqz]Emnpq2=0.
3.6

Solving this linear equation system, one can obtain the eigentrain terms of E0, E1and E2, with which the displacement is written as

ui(x)=ΩΓimnCmnkl(Ekl0+Eklp1(xpxpc)+Eklpq2(xpxpc)(xqxqc))dx=14πμ0(Mimn0CmnklEkl0+Npimn1CmnklEklp1+Tpqimn2CmnklEklpq2),
3.7

where the explicit forms of Mimn0, Npimn1 and Tpqimn2 are available in the electronic supplementary material. The strain is calculated with equation (2.17), and the stress can be obtained by the constitutive law.

(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 xI=(x1I,x2I,x3I). All particles being treated as including the same material properties as the matrix, the eigenstrain can be introduced with a polynomial form as

ϵij,I(x)=Eij0,I+Eijk1,I(xkxkI)+Eijkl2,I(xkxkI)(xlxlI)+,xΩI,
3.8

where no summation is applied to the superscript, while repeated subscripts still mean summation. Then the equivalent condition equation (3.4) for one of the particles, say ΩJ, reads

Cijkl1(ϵkl0+ϵkl)=Cijklϵkl0+ϵklϵkl,J,xΩJ.
3.9

Figure 3.
A semi-infinite domain D (x3≥0) containing many particles. (Online version in colour.)

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:

ui(x)=14πμ0I=1N(Mimn0,ICmnklEkl0,I+Npimn1,ICmnklEklp1,I+Tpqimn2,ICmnklEklpq2,I).
3.10

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 (xrxrJ) terms:

I=1NΔCijklDklst0(xc)CstmnEmn0,I4πμ0CijmnEmn0,J+I=1NΔCijklDpklst1,I(xc)CstmnEmnp1,I+I=1NΔCijklDpqklst2,I(xc)CstmnEmnpq2,I=4πμ0ΔCijklϵkl0I=1NΔCijklDpklst,r1,I(xc)CstmnEmnp1,I4πμ0CijmnδprEmnp1,J+I=1NΔCijklDklst,r0,I(xc)CstmnEmn0,I+I=1NΔCijklDpqklst,r2,I(xc)CstmnEmnpq2,I=0andI=1NΔCijklDklst,yz0,I(xc)CstmnEmn0,I+I=1NΔCijklDpklst,yz1,I(xc)CstmnEmnp1,I+I=1NΔCijklDpqklst,yz2,I(xc)CstmnEmnpq2,I8πμ0CijmnδpyδqzEmnpq2,J=0.
3.11

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,

[A]60×60nJ[E]60n=[F]60J,
3.12

where the eigenstrain terms of all the particles are represented by a vector as [E] whose length is 60n. [A]J is the local coefficient matrix with dimensions 60×60n 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

[A]60n×60n[E]60n=[F]60n.
3.13

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 60N 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]60×60nJ in equation (3.12) is calculated, which will make the global coefficient matrix [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 60N 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 E1=69 GPa, E0=0.8 GPa, ν1=0.334 and ν0=0.38, respectively [38,46,47]. The strain components corresponding to the uniform load are ϵ330=0.001, ϵ110=ϵ220=ν0ϵ330 is applied to the system to represent a far-field tension σ330=0.8 MPa. Under this far-field tension, the elastic field is investigated.

Figure 4 illustrates the distribution of eigenstrain component ϵ33 inside the inhomogeneity at the x1x3 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 4a), linearly (cf. figure 4b) and quadratically distributed (cf. figure 4c). 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 x3=0. In total, about 92 000 unstructured 6-node triangle elements are used to calculate the solutions.

Figure 4.
Distribution of ϵ33 in the inclusion: (a) uniform distribution; (b) linear distribution; and (c) quadratic distribution. (Online version in colour.)

Figure 5a illustrates the variation of stress σ33 along the axis x3 (x1=x2=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 5b), 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.

Figure 5.
Variation of stress σ33 with x3: (a) h = 1.1a and (b) h =2a. (Online version in colour.)

In figure 5a, σ33 is equal to the far-field tension at the boundary x3=0, and decreases with x3 to the bottom of the particle. Inside the particle, it increases with x3 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 2σ330. Then it decreases with x3 and finally converges to the far-field tension σ330. Figure 5b illustrates the variation of σ33 with x3, 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 σ330 at the boundary, and then increases with x3 and reaches the peak value near the bottom of the particle. Then it decreases with x3 to the bottom of the particle. Inside the particle, σ33 increases with x3. 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 x3 and converges to the far-field load. In figure 5b, 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 σ310=0.58MPa is applied to the above system, in which h=11 μm. Figure 6a illustrates the variation of the shear stress with x3. The shear stress σ31 equals to the far-field shear stress σ310 at the surface, and then increases rapidly with x3 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 x3 fast and reaches the minimum value, i.e. ~0.75 σ310. It increases again and then converges to the far-field shear load.

Figure 6.
Variation of shear stress σ31 with x3: (a) h=11 μm for different forms of eigenstrain distribution and (b) quadratic form of eigenstrain for different h values. (Online version in colour.)

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 x3 curves for three cases of h=11, 20 and 30 μm are illustrated in figure 6b. The shear stresses at the surface are equal to the far-field load. When h=11 μm, the stress increases with x3 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 x3 axis and then increases to the inflection point at the bottom of the particle. Inside the particle, σ31 increases with x3. 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 2a to avoid the particles' overlap. The same uniaxial load of σ330=0.8MPa is applied to the system corresponding to ϵ330=0.001, ϵ110=ϵ220=ν0ϵ330.

First, an extreme case is considered first with h=2.0a and s=2.5a. The results from EIM using uniform, linear and quadratic forms of eigenstrains are compared with the FEM results in figure 7a. 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 x3 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 x3 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%.

Figure 7.
Variation of stress σ33 with x3 when two top-down particles are embedded into a semi-infinite domain: (a) s=2.5a and (b) s=3.0a. (Online version in colour.)

When the distance between two particles is longer, the eigenstrain is simpler owing to the decay of particle–particle interaction. Figure 7b illustrates the variation of σ33 with x3 when s=3.0a. 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 7a), the maximum stress locates in the middle of the two particles, whereas when their distance is longer (cf. figure 7b), 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 σ330=0.8 MPa, the stress profile is investigated. Figure 8a illustrates the variation of σ11 with x1, for x3=h. In this case, s=2.5a while h=1.1a. One can observe that σ11 is negative, which indicates that the particle interaction produces compression between particles in the x1 direction. The curve is symmetric. For x10, |σ11| is zero at far field, and then increases with x1 and reaches the local peak value near the left particle at about 0.45σ330. Then |σ11| decreases to a local minimum value, say 0.34σ330 in the particle. After that, |σ11| increases fast and reaches the maximum value 1.47σ330 at the middle point between two particles. One can conclude that the centre of the two particles is under high compression owing to the particle interaction. The EIM results with uniform, linear and quadratic forms of eigenstrain are compared. The maximum values from uniform and linear distributed eigenstrain are much smaller than that for quadratic ones. Because 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.
Variation of stress (a) σ11 and (b) σ33 versus x1/a (a=10 μm). (Online version in colour.)

Figure 8b illustrates the symmetric distribution of σ33 with x1. For x10, the stress σ33 decrease from σ330 at far field to a small negative value at the boundary of the left particle. Then a large jump to 1.71σ330 is observed at the particle boundary. Then it decreases with x1 inside the particle. A large drop is also observed at the other boundary of the particle, where σ33 jumps from approximately 1.43σ330 to approximately 0.69σ330. Finally, it increases to a small positive value at the centre of the two particles. The three forms of eigenstrain distribution exhibit notable differences on the particle domain but the difference between the linear and quadratic forms is smaller than that between the uniform and linear forms, which indicates the quadratic form solution should be close to the convergent solution.

To investigate the particle interactions and boundary effects, different configurations of h/a and s/a are considered. Figure 9a illustrates the stress σ33 distribution along the x1 in the case of x3=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.1a, h=2.0a and h=3.0a, the stress σ33 inside the particle caused by σ330 increases with h. The difference of σ33 between h=2.0a and 3.0a is much less than that between h=1.1a and h=2.0a, which means that the boundary effect on the local field in the neighbourhood of particles is small for h≥2.0a. Moreover, σ33 outside the particle are almost the same for the three cases. Figure 9b illustrates the stress σ33 distribution for different centre–centre distances, i.e. s=2.5a, s=3.0a, s=4.0a while h=1.1a. Stress component σ33 is greater than zero at the midpoint of two particles, i.e., x1=0. And the value increases with s. When s=3.0a, σ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.

Figure 9.
Variation of σ33 with x1: (a) comparison of different h and (b) comparison of different s. (Online version in colour.)

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.1a, h=2.0a and h=3.0a. 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=2a, although only a tension stress σ330 is applied, the induced stress σ11 is near 4.0σ330 and while σ33 is near 3.0σ330, which indicates the centre of two particles is under compression. With the increase in s, σ11 and σ33 will increase and finally approach 0 and σ330, respectively. It is interesting that σ33 is negative for s=2a−2.4a and then turns to positive for s>2.4a, which indicates the stress direction can change with the particle–particle distance.

Figure 10.
Variation of stress with s at midpoint of particles for (a) σ11 and (b) σ33. (Online version in colour.)

(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 11a. 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.

Figure 11.
(a) Configuration of a simple cubic lattice in a semi-infinite domain. The dashed lines are just used to distinguish the unit cells, not boundaries. (b) The configuration of neighbouring particles whose effects are considered. Here only the first three ...

Given a simple lattice composite in figure 11a, 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 E1=69 GPa, E0=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=2a. The centre–centre distance of neighbouring particle is s=4a. 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 23a.

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 11b. 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 ϵij embedded in a large domain D with the free stress boundary condition. The stress and displacement field in the domain are continuous. The average stress can be written as

σijD=1VDDσij(x)dx=1VDDσik(x)nkxjdx=σik0VDDnkxjdx=σij0,
4.1

where the symmetric and periodic boundary conditions along the opposite lateral surfaces and equilibrium in the vertical direction are used.

Similarly, the average strain can be written as

ϵijD=1VDDϵij(x)dx=1VDD(ϵijϵij)(x)dx+1VDDϵij(x)dx=Cijkl1VDDσkl(x)dx+1VDDϵijdx=Cijkl1σkl0+ηϵij,
4.2

where η 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

ϵijDI=Cijkl1σkl0+1VDΩ[Eij0,I+Eijk1,I(xkxkI)+Eijkl2,I(xkxkI)(xlxlI)]dx=Cijkl1σkl0+ηEij0,I+aI25Eijkk2,I,
4.3

where aI 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:

σ33D=σ330;ϵijDI=Cij331σ330+ηEij0,I+aI25Eijkk2,I.
4.4

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

E=σ330/ϵ33DIandv=ϵ11D/ϵ33D.
4.5

Once the testing load is given, the elastic moduli only depend on the eigenstrain.

Under the far-field tension σ330=0.8MPa, the eigenstrain component ϵ33 distribution at x1x3 planes are illustrated in figure 12. One can observe that due to the boundary effect, the ϵ33 inside unit cell 1 is very different from the others. The distribution along x1 direction is symmetric, but the symmetry about x3 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, ϵ33 is highest at the bottom, and then decreases along the x3 direction. The differences of ϵ33 at the bottom and top are fairly large. For unit cell 2, ϵ33 is also decreasing along the x3 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 x3 axis becomes nearly symmetric. Noted that with the boundary effect decaying, not only the symmetry along x3 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.

Figure 12.
Distribution of ϵ33 in unit cells (a) unit cell 1; (b) unit cell 2; (c) unit cell 3; (d) unit cell 4; (e) unit cell 5; (f) unit cell 6. (Online version in colour.)

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.

Figure 13.
Variation of effective material properties for unit cells: (a) Young's modulus and (b) Possion's ratio. (Online version in colour.)

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 13b 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 σ310 at the boundary, the average stress and average strain can be obtained by rewriting Eqs (4.1) and (4.2), i.e.

σ31D=σ310;ϵijDI=Cij311σ310+ηEij0,I+a25Eijkk2,I.
4.6

The effective shear modulus left angle bracketμright angle bracket can be calculated as

μ=σ312ϵ31.
4.7

Under a uniform shear loading σ310=0.58MPa, the eigenstrain component ϵ31 distributions at x1x3 planes are illustrated in figure 14. The eigenstrain ϵ31 in unit cell 1 is quite different from that of other unit cells. ϵ31 is symmetric about the x1 axis. But the symmetry about the x3 axis does not exist due to the boundary effect. The eigenstrain ϵ31 is highest at the bottom of the particle, and then decreases along the x3 direction. From unit cell 2 to unit cell 6, the symmetry about x3 is gradually recovering. This indicates that the boundary effect is decaying with the increasing distance to the boundary, as for unit cell 6, ϵ31 is almost symmetric along the x3 axis. Eigenstrain distributions in unit cell 5 and unit cell 6 are almost the same, which means that the eigenstrain is convergent.

Figure 14.
Distribution of ϵ31 in unit cells: (a) unit cell 1; (b) unit cell 2; (c) unit cell 3; (d) unit cell 4; (e) unit cell 5; (f) unit cell 6. (Online version in colour.)

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.

Figure 15.
Variation of shear modulus for unit cell. (Online version in colour.)

(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 E1>E0 and G1>G0. Figure 16a 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 x3 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.
Variation of effective material properties for unit cells for different particle size: (a) Young's modulus and (b) shear modulus. Note that different colour inhomogeneities represent the particles with different sizes. (Online version in colour.)

Figure 16b 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 17a. 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 17b 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.

Figure 17.
Variation of effective material properties with volume fraction: (a) Young's modulus and (b) shear modulus. (Online version in colour.)

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.

Supplementary Material

supplementary_material.tex:

Supplementary Material

EIM_Library.cpp:

Acknowledgements

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

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.

References

1. Eshelby JD. 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. A 241, 376–396. (doi:10.1098/rspa.1957.0133)
2. Eshelby JD. 1959. The elastic field outside an ellipsoidal inclusion. Proc. R. Soc. Lond. A 252, 561–569. (doi:10.1098/rspa.1959.0173)
3. Takao Y, Chou T-W, Taya M. 1982. Effective longitudinal Young's modulus of misoriented short fibre composites. J. Appl. Mech. 49, 536–540. (doi:10.1115/1.3162509)
4. Chen C, Cheng C. 1996. Effective elastic moduli of misoriented short-fibre composites. Int. J. Solids Struct. 33, 2519–2539. (doi:10.1016/0020-7683(95)00160-3)
5. Yin HM, Buttlar WG, Paulino GH, Di Benedetto H. 2008. Assessment of existing micromechanical models for asphalt mastics considering viscoelastic effects. Road Mater. Pavement 9, 31–57. (doi:10.1080/14680629.2008.9690106)
6. Dunn ML, Taya M. 1993. Micromechanics predictions of the effective electroelastic moduli of piezoelectric composites. Int. J. Solids Struct. 30, 161–175. (doi:10.1016/0020-7683(93)90058-F)
7. Hatta H, Taya M. 1986. Equivalent inclusion method for steady-state heat conduction in composites. Int. J. Eng. Sci. 24, 1159–1172. (doi:10.1016/0020-7225(86)90011-X)
8. Yin HM, Paulino G, Buttlar WG, Sun LZ. 2005. Effective thermal conductivity of functionally graded composites. J. Appl. Phys. 98, 063704 (doi:10.1063/1.2039998)
9. Yin HM, Paulino G, Buttlar WG, Sun LZ. 2008. Effective thermal conductivity of functionally graded particulate nanocomposites with interfacial thermal resistance. J. Appl. Mech. 75, 051113 (doi:10.1115/1.2936893)
10. Takei T, Hatta H, Taya M. 1991. Thermal expansion behavior of particulate-filled composites II: multi-reinforcing phases (hybrid composites). Mat. Sci. Eng. A 131, 145–152. (10.1016/0921-5093(91)90353-O)
11. Takei T, Hatta H, Taya M. 1991. Thermal expansion behavior of particulate-filled composites I: single reinforcing phase. Mat. Sci. Eng. A 131, 133–143. (doi:10.1016/0921-5093(91)90352-N)
12. Yin HM, Paulino G, Buttlar W, Sun LZ. 2007. Micromechanics-based thermoelastic model for functionally graded particulate materials with particle interactions. J. Mech. Phys. Solids 55, 132–160. (doi:10.1016/j.jmps.2006.05.002)
13. Sakata S, Ashida F, Kojima T. 2010. Stochastic homogenization analysis for thermal expansion coefficients of fiber reinforced composites using the equivalent inclusion method with perturbation-based approach. Comput. Struct. 88, 458–466. (doi:10.1016/j.compstruc.2009.12.007)
14. Duan HL, Wang J, Huang ZP, Karihaloo BL. 2005. Eshelby formalism for nano-inhomogeneities. Proc. R. Soc. A 461, 3335–3353. (doi:10.1098/rspa.2005.1520)
15. Gao X, Huang ZP, Qu JM, Fang DN. 2014. A curvature-dependent interfacial energy-based interface stress theory and its applications to nano-structured materials: (I) general theory. J. Mech. Phys. Solids 66, 59–77. (doi:10.1016/j.jmps.2014.01.010)
16. Xiao XZ, Song DK, Xue JM, Chu HJ, Duan HL. 2015. A self-consistent plasticity theory for modeling the thermo-mechanical properties of irradiated FCC metallic polycrystals. J. Mech. Phys. Solids 78, 1–16. (doi:10.1016/j.jmps.2015.01.011)
17. Mura T. 1987. Micromechanics of defects in solids. Dordrecht, The Netherlands: Kluwer Academic Publishers.
18. Nemat-Nasser S, Hori M. 1999. Micromechanics: overall properties of heterogeneous materials. Amsterdam, The Netherlands: Elsevier.
19. Hill R. 1963. Elastic properties of reinforced solids-some theoretical principles. J. Mech. Phys. Solids 11, 357–372. (doi:10.1016/0022-5096(63)90036-X)
20. Willis J. 1981. Variational and related methods for the overall properties of composites. Adv. Appl. Mech. 21, 1–78. (doi:10.1016/S0065-2156(08)70330-2)
21. Yin HM, Sun LZ, Chen JS. 2002. Micromechanics-based hyperelastic constitutive modeling of magnetostrictive particle-filled elastomers. Mech. Mater. 34, 505–516. (doi:10.1016/S0167-6636(02)00178-3)
22. Xia Z, Zhang Y, Ellyin F. 2003. A unified periodical boundary conditions for representative volume elements of composites and applications. Int. J. Solids Struct. 40, 1907–1921. (doi:10.1016/S0020-7683(03)00024-6)
23. Yin HM, Sun LZ. 2005. Elastic modelling of periodic composites with particle interactions. Phil. Mag. Lett. 85, 163–173. (doi:10.1080/09500830500157413)
24. Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. 2003. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int. J. Solids Struct. 40, 3647–3679. (doi:10.1016/S0020-7683(03)00143-4)
25. Yin HM, Sun LZ, Chen JS. 2006. Magneto-elastic modeling of composites containing chain-structured magnetostrictive particles. J. Mech. Phys. Solids 54, 975–1003. (doi:10.1016/j.jmps.2005.11.007)
26. Reiter T, Dvorak GJ, Tvergaard V. 1997. Micromechanical models for graded composite materials. J. Mech. Phys. Solids 45, 1281–1302. (doi:10.1016/S0022-5096(97)00007-0)
27. Yin HM, Sun LZ, Paulino G. 2004. Micromechanics-based elastic model for functionally graded materials with particle interactions. Acta Mater. 52, 3535–3543. (doi:10.1016/j.actamat.2004.04.007)
28. Weng G. 1984. Some elastic properties of reinforced solids, with special reference to isotropic ones containing spherical inclusions. Int. J. Eng. Sci. 22, 845–856. (doi:10.1016/0020-7225(84)90033-8)
29. Drugan W, Willis J. 1996. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. J. Mech. Phys. Solids 44, 497–524. (doi:10.1016/0022-5096(96)00007-5)
30. Fish J. 2007. Discrete-to-continuum scale bridging. In Multiscaling in molecular and continuum mechanics: interaction of time and size from macro to nano (ed. Sih GC, editor. ), 85–102. Dordrecht, The Netherlands: Springer.
31. Aboudi J, Arnold SM, Bednarcyk BA. 2013. Micromechanics of composite materials: a generalized multiscale analysis approach. Amsterdam, The Netherlands: Elsevier.
32. Gusev AA. 1997. Representative volume element size for elastic composites: a numerical study. J. Mech. Phys. Solids 45, 1449–1459. (doi:10.1016/S0022-5096(97)00016-1)
33. Mindlin RD. 1936. Force at a point in the interior of a semi-infinite solid. J. Appl. Phys. 7, 195 (doi:10.1063/1.1745385)
34. Rongved L. 1955. Point force in the interior of a semi-infinite solid with fixed boundary. J. Appl. Mech. 22, 545–546.
35. Rongved L. 1955. Force interior to one of two jointed semi-infinite solids. In Proc. 2nd Midwestern Conf. on Solid Mechanics, Lafayette, IN, 23 April, pp. 1–13. Lafayette, IN: Purdue University.
36. Yu HY, Sanday SC. 1991. Elastic field in joined semi-infinite solids with an inclusion. Proc. R. Soc. Lond. A 434, 521–530. (doi:10.1098/rspa.1991.0111)
37. Walpole L. 1996. An elastic singularity in joined half-spaces. Int. J. Eng. Sci. 34, 629–638. (doi:10.1016/0020-7225(95)00120-4)
38. Liu YJ, Yin HM. 2014. Stress concentration of a microvoid embedded in an adhesive layer during stress transfer. J. Eng. Mech. 140, 04014075 (doi:10.1061/(ASCE)EM.1943-7889.0000786)
39. Liu YJ, Yin HM. 2014. Equivalent inclusion method-based simulation of particle sedimentation toward functionally graded material manufacturing. Acta Mech. 225, 1429–1445. (doi:10.1007/s00707-013-1058-0)
40. Yin HM, Lee P-H, Liu YJ. 2014. Equivalent inclusion method for the Stokes flow of drops moving in a viscous fluid. J. Appl. Mech. 81, 071010 (doi:10.1115/1.4027312)
41. Mura T, Shodja H. 1996. Inclusion problems. Appl. Mech. Rev. 118–127. (doi:10.1016/0020-7225(95)00120-4)
42. Avazmohammadi R, Hashemi R, Shodja HM, Kargarnovin MH. 2010. Ellipsoidal domain with piecewise nonuniform eigenstrain field in one of joined isotropic half-spaces. J. Elast. 98, 117–140. (doi:10.1007/s10659-009-9220-6)
43. Boussinesq J. 1885. Applications des potentiels à l' étude de l'équilibre et du mouvement des solides élastiques. Paris, France: Gauthier-Villars.
44. Mindlin RD. 1953. Force at a point in the interior of a semi-infinite solid. In Proc. 1st Midwestern Conf. Solid Mech, Urbana, IL, 24–24 April, pp. 56–59. Urbana, IL: Hunter & Wroughton, Inc.
45. Walpole L. 1997. An inclusion in one of two joined isotropic elastic half-spaces. IMA J. Appl. Math. 59, 193–209. (doi:10.1093/imamat/59.2.193)
46. Yin HM, Yang DJ, Kelly G, Garant J. 2013. Design and performance of a novel building integrated PV/thermal system for energy efficiency of buildings. Sol. Energy 87, 184–195. (doi:10.1016/j.solener.2012.10.022)
47. Lee P-H, Yin HM. 2014. Size effect on functionally graded material fabrication by sedimentation. J. Nanomech. Micromech. 5, A4014008 (doi:10.1061/(ASCE)NM.2153-5477.0000087)
48. Moschovidis ZA, Mura T. 1975. Two ellipsoidal inhomogeneities by equivalent inclusion method. J. Appl. Mech. 42, 847–852. (doi:10.1115/1.3423718)
49. Fish J. 2014. Practical multiscaling. Chichester, UK: Wiley.

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society