# Eigenstrain Boundary Integral Equations with Local Eshelby Matrix for Stress Analysis of Ellipsoidal Particles

Authors: Hang Ma; Cheng Yan; Qing-hua Qin
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Aiming at the large scale numerical simulation of particle reinforced materials, the concept of local Eshelby matrix has been introduced into the computational model of the eigenstrain boundary integral equation (BIE) to solve the problem of interactions among particles. The local Eshelby matrix can be considered as an extension of the concepts of Eshelby tensor and the equivalent inclusion in numerical form. Taking the subdomain boundary element method as the control, three-dimensional stress analyses are carried out for some ellipsoidal particles in full space with the proposed computational model. Through the numerical examples, it is verified not only the correctness and feasibility but also the high efficiency of the present model with the corresponding solution procedure, showing the potential of solving the problem of large scale numerical simulation of particle reinforced materials.


---

## Body

## 1. Introduction

The inclusion and inhomogeneity problems have been a focus of solid mechanics for several decades since the pioneering work of Eshelby in the fifties of last century [1]. The elastic behavior of an inclusion embedded in a matrix is of considerable importance in a wide variety of physical and engineering problems [2]. Following Eshelby’s work on the eigenstrain solution and equivalent inclusion, numerous investigations have been carried out and reported in the literature [3–5]. The eigenstrain solution can correspond to thermal mismatch, lattice mismatch, phase transformation, microstructural evolution, quantum dot [6], intrinsic strains in residual stress problems [7], and so forth. With the replacement of equivalent inclusion, the eigenstrain solution can correspond also to various inhomogeneity [5], cavity, and even the crack problems [8], showing the great significance of this research.The analytical methods can afford the fine details of stress distributions within or outside the particles and the basis for further investigation. However, the available solutions apply generally to only simple geometries such as single ellipsoidal, cylindrical, and spherical inclusions in full or half space. Therefore, the FEM based numerical methods [9] as well as the boundary element methods (BEM) [10–12] become the tools for the analysis of particle problems with various shapes and materials. For solids containing a large number of particles with random size and space distributions, the difficulty of numerical methods lies mainly in the fact that the solution scale is too large owing to the description of interfaces in discretization. One of the procedures for the difficulty using the BEM is to introduce the special technique of fast multiple expansions [13, 14], which leads to much complexity of the algorithm.Recently a new computational model with an iteration procedure has been proposed by introducing the concepts of Eshelby’s eigenstrain and the equivalent inclusion into the boundary integral equations (BIE) for the 2D and 3D stress analysis and the overall properties of solids with a large number of particles [15–17]. For the densely populated particles, however, the interaction among particles will affect the convergence of iteration, depending on the distance between particles. The shorter the distance is, the stronger the interaction will be. In order to overcome this difficulty, the local Eshelby matrix for the newly defined group of near field particles has been introduced with the aid of the discrete form of eigenstrain BIE in full space to improve the original algorithm. Taking the subdomain BEM in full space [18] as the control, the 3D stress analysis has been carried out for several ellipsoidal particles with various Young’s modulus ratios and different shapes to verify the feasibility and efficiency of the improved eigenstrain BIE algorithm.


## 2. Computational Models

### 2.1. Eigenstrain Boundary Integral Equations

In the present model, the perfect adhesion between the particle and matrix, both being isotropic materials, is assumed; that is, the displacement continuity and the traction equilibrium hold true along their interfaces. In the solution domain considered, Ω represents the domain of matrix with the outer boundaryΓ, and ΩI is the domain of particle with the boundary or interface ΓI (ΓI=ΩI∩Ω). The displacement and the traction eigenstrain boundary integral equations are, respectively, as follows:
(1)C(p)ui(p)=∫Γτj(q)uij*(p,q)dΓ(q)hhhhhhhhhii-∫Γuj(q)τij*(p,q)dΓ(q)hhhhhhhhhii+∑I=1NI∫ΩIεjk0(q)σijk*(p,q)dΩ(q),(2)C(p)σij(p)=∫Γτk(q)uijk*(p,q)dΓ(q)hhhhhhhhhhh-∫Γuk(q)τijk*(p,q)dΓ(q)hhhhhhhhhhh+∑I=1NI∫ΩIεkl0(q)σijkl*(p,q)dΩ(q)hhhhhhhhhhh+εkl0(p)Oijkl*(p,q),
where
(3)Oijkl*(p,q)=lim⁡Ωε→0∫ΔΓεxlτijk*(p,q)dΓ(q)represents the free term resulted from the domain integral of (2). Ωε with the boundary Γε is the infinitesimal region ε in ΩI. xl=xl(q)-xl(p) stands for the projection of the distance between the field and source points, q and p. In (1) and (2), uij*, τij*, and σij* are the Kelvin fundamental solutions of the displacement, the traction, and the stress, respectively. uijk*, τijk*, and σijk* are the corresponding derivations of the fundamental solutions. NI is the total particle number in Ω. εij0 is the eigenstrain. In fact, the boundary integral equations describe only the displacement and stress fields of homogeneous elastic media. In order to describe the solids with inhomogeneous particles, the replacement by the equivalent inclusion has to be carried out with the aid of the concept of Shelby tensors.


### 2.2. Eshelby Tensor and Equivalent Inclusions

According to the work of Eshelby [1], the tress-free constrained strain εijC and the eigenstrain εij0 of a single particle in full space are correlated with Eshelby tensor, Sijkl; that is,
(4)εijC=Sijklεkl0.The value of Eshelby tensor depends on the shape ofΩI. For simple geometries, Eshelby tensor can be expressed explicitly or found in the literature. In general case, Eshelby tensor can be expressed in integral form [17] and computed numerically:
(5)Sijkl=12(δikδjl+δilδjk)iiiiiiiiii+14μ∫ΓIxl{τijk*+τjik*-2ν1+ντijτmmk*}dΓ,
where μ and ν are the shear modulus and Poisson’s ratio of the matrix, respectively. It needs to be pointed out that (5) is applicable only for the uniform distribution of eigenstrain in ΩI. For the multiple particles, the condition for the uniform distribution is that the volume of particles would be small enough or the distance between particles is kept sufficiently large. In the present work, the uniform distribution of eigenstrain is assumed. Define the Young’s modulus ration =EI/EM, where the subscripts I and M denote the particle and matrix, respectively. According to Hooke’s law, if a particle under the applied strain εij is replaced by an equivalent inclusion without altering its stress state, the following relation should be satisfied:
(6)(1-β1)εijC+β2δijεkkC-εij0-νM1-2νMδijεkk0=-(1-β1)εij-β2δijεkk,
where
(7)β1=1+νM1+νIβ,hhhiiiiiiiiβ2=νM1-νM-νI1-2νIβ.With (4) and (6), the eigenstrain of a single particle can be computed using its applied strains so that the replacement of the equivalent inclusion is realized. In this way, the solids containing particles become apparently homogeneous which can be described by the boundary integral equations (1) and (2). For a single particle, the matrix form of (6) after discretization can be written as
(8)[S0]{ε0}={ε}.


### 2.3. Local Eshelby Matrix

As stated previously, the eigenstrains of a particle need to be determined by using the applied strains of the particle when the replacement is carried out by an equivalent inclusion. However, for the case of multiple particles, there are interactions among particles. In addition to its own state of the particle, the applied strains on it will be affected by the surrounding particles owing to the interactions among them. The gratitude of the interactions depends on the distances of particles, which would be a primary factor interfering the convergence of iteration procedures [15–17]. As shown in Figure 1, the group of the near-field particles having a number NL is defined as those within the circle of dashed line for the current particle I while the group of the far-field particles can be defined correspondingly as those outside the circle. With such definitions, the near-field particles belong to the short distance group with relatively strong effect of interactions while the far-field particles belong to the long distance group with relatively weak effect of interactions to the current particle.Figure 1
The group definitions for multiple particles.In the full space, if only the near-field group is taken into consideration while neglecting tentatively the far-field group, the stresses of the current particleI can be expressed as follows, using the stress BIE (2) by combining the constitutive relation as well as the integral type transformation [17]:
(9)σijI(p)=EijklεklI=Eijkl(εklCI-εkl0I)iiiiiiii=∑J=1NL∫ΓJεkl0Jxlτijk*dΓ,p∈ΩI,
where Eijkl stands for the elastic tensor of matrix. As there is no outer boundary in full space, the two boundary integrals in (2) vanish. By making extension of the concept of Eshelby tensor to different particles and using (9) and (5), the constrained strains of the current particle and the corresponding Eshelby tensor can be written, respectively, as
(10)εijCI=εij0I+∑J=1NLCijklJ∫ΓJεmn0Jxnτklm*dΓ=∑J=1NLSijklIJεkl0J,(11)SijklIJ=12(δikδjl+δilδjk)δIJiiiiiiiii+14μ∫ΓJxl{τijk*+τjik*-2ν1+νδijτmmk*}dΓ.If the superscriptsI≠J, the Eshelby tensor SijklIJ correlates the eigenstrains and the constrained strains of different particles in the near-field group. Similar to the relation of replacement of (6), according to Hooke’s law, the relation of replacement by equivalent inclusions should be followed for the particles in the near-field group:
(12)(1-β1)εijCI+β2δijεkkCI-εij0I-ν1-2νδijεkk0I=-(1-β1)εijI-β2δijεkkI,I=1,2,…,NL.Combining (10) and (11) and after discretization, the matrix form of (12) can be written as
(13)[S110S12⋯S1,NLS21S220S2,NL⋮⋱⋮SNL,1SNL,2⋯SNL,NL0]{ε10ε20⋮εNL0}={ε1ε2⋮εNL}.In this way, with (13), the eigenstrain vector can be computed using the applied strain vector for the near-field particles. Similar to the analysis for the multiple cracks [8], the matrix in (13) is named Eshelby matrix, which can be looked at as an extension of the concept of Eshelby tensor and the equivalent inclusion in the numerical applications. It needs to be pointed out that (9)–(12) are analytical equations, which reflect special stress-strain correlations among particles of the near-field group in full space, describing their state precisely. Nevertheless, (13) is a numerical form after discretization, describing the state and correlations of the near-field particles approximately.For a problem to be solved, once the radius is given for the dashed line circle as shown in Figure1, each particle in the solid will have, in general case, a unique near-field group and a unique Eshelby matrix correspondingly. For the convenience of computing, rewrite (13) as
(14)ε(I)0=SIε(I),I=1,2,…,NI,
where the matrix SI in (14) is obtained by the inversion then reduction of Eshelby matrix, which corresponds to the eigenstrain vector, ε(I)0, of the current particle. ε(I)={ε1,ε2,…,εNL}T is the applied strain vector of the near-field particles, being the same as that of the right hand side of (13), corresponding to the current particle. The subscript I in (14) goes through all of the particles in the solid; the total number of them being denoted as  NI.In this way, based on the proposed algorithm of the eigenstrain BIE with Eshelby matrix, the computation of eigenstrains of each particle in solids consists of three parts: the first is computed from the stresses or the applied strains resulted from the load. The second comes from the effects of particles in the near-field group, having relatively strong interactions, which is computed using (14). The third is from the effects of particles in the far-field group, having relatively weak interactions, which is computed directly using the domain integrals of the BIE. It is obvious that the minimum number of particles in the near-field group is one. If so (NL= 1), the proposed algorithm would be reduced to that in the previous work [15–17], while (14) becomes (8).


## 2.1. Eigenstrain Boundary Integral Equations

In the present model, the perfect adhesion between the particle and matrix, both being isotropic materials, is assumed; that is, the displacement continuity and the traction equilibrium hold true along their interfaces. In the solution domain considered, Ω represents the domain of matrix with the outer boundaryΓ, and ΩI is the domain of particle with the boundary or interface ΓI (ΓI=ΩI∩Ω). The displacement and the traction eigenstrain boundary integral equations are, respectively, as follows:
(1)C(p)ui(p)=∫Γτj(q)uij*(p,q)dΓ(q)hhhhhhhhhii-∫Γuj(q)τij*(p,q)dΓ(q)hhhhhhhhhii+∑I=1NI∫ΩIεjk0(q)σijk*(p,q)dΩ(q),(2)C(p)σij(p)=∫Γτk(q)uijk*(p,q)dΓ(q)hhhhhhhhhhh-∫Γuk(q)τijk*(p,q)dΓ(q)hhhhhhhhhhh+∑I=1NI∫ΩIεkl0(q)σijkl*(p,q)dΩ(q)hhhhhhhhhhh+εkl0(p)Oijkl*(p,q),
where
(3)Oijkl*(p,q)=lim⁡Ωε→0∫ΔΓεxlτijk*(p,q)dΓ(q)represents the free term resulted from the domain integral of (2). Ωε with the boundary Γε is the infinitesimal region ε in ΩI. xl=xl(q)-xl(p) stands for the projection of the distance between the field and source points, q and p. In (1) and (2), uij*, τij*, and σij* are the Kelvin fundamental solutions of the displacement, the traction, and the stress, respectively. uijk*, τijk*, and σijk* are the corresponding derivations of the fundamental solutions. NI is the total particle number in Ω. εij0 is the eigenstrain. In fact, the boundary integral equations describe only the displacement and stress fields of homogeneous elastic media. In order to describe the solids with inhomogeneous particles, the replacement by the equivalent inclusion has to be carried out with the aid of the concept of Shelby tensors.


## 2.2. Eshelby Tensor and Equivalent Inclusions

According to the work of Eshelby [1], the tress-free constrained strain εijC and the eigenstrain εij0 of a single particle in full space are correlated with Eshelby tensor, Sijkl; that is,
(4)εijC=Sijklεkl0.The value of Eshelby tensor depends on the shape ofΩI. For simple geometries, Eshelby tensor can be expressed explicitly or found in the literature. In general case, Eshelby tensor can be expressed in integral form [17] and computed numerically:
(5)Sijkl=12(δikδjl+δilδjk)iiiiiiiiii+14μ∫ΓIxl{τijk*+τjik*-2ν1+ντijτmmk*}dΓ,
where μ and ν are the shear modulus and Poisson’s ratio of the matrix, respectively. It needs to be pointed out that (5) is applicable only for the uniform distribution of eigenstrain in ΩI. For the multiple particles, the condition for the uniform distribution is that the volume of particles would be small enough or the distance between particles is kept sufficiently large. In the present work, the uniform distribution of eigenstrain is assumed. Define the Young’s modulus ration =EI/EM, where the subscripts I and M denote the particle and matrix, respectively. According to Hooke’s law, if a particle under the applied strain εij is replaced by an equivalent inclusion without altering its stress state, the following relation should be satisfied:
(6)(1-β1)εijC+β2δijεkkC-εij0-νM1-2νMδijεkk0=-(1-β1)εij-β2δijεkk,
where
(7)β1=1+νM1+νIβ,hhhiiiiiiiiβ2=νM1-νM-νI1-2νIβ.With (4) and (6), the eigenstrain of a single particle can be computed using its applied strains so that the replacement of the equivalent inclusion is realized. In this way, the solids containing particles become apparently homogeneous which can be described by the boundary integral equations (1) and (2). For a single particle, the matrix form of (6) after discretization can be written as
(8)[S0]{ε0}={ε}.


## 2.3. Local Eshelby Matrix

As stated previously, the eigenstrains of a particle need to be determined by using the applied strains of the particle when the replacement is carried out by an equivalent inclusion. However, for the case of multiple particles, there are interactions among particles. In addition to its own state of the particle, the applied strains on it will be affected by the surrounding particles owing to the interactions among them. The gratitude of the interactions depends on the distances of particles, which would be a primary factor interfering the convergence of iteration procedures [15–17]. As shown in Figure 1, the group of the near-field particles having a number NL is defined as those within the circle of dashed line for the current particle I while the group of the far-field particles can be defined correspondingly as those outside the circle. With such definitions, the near-field particles belong to the short distance group with relatively strong effect of interactions while the far-field particles belong to the long distance group with relatively weak effect of interactions to the current particle.Figure 1
The group definitions for multiple particles.In the full space, if only the near-field group is taken into consideration while neglecting tentatively the far-field group, the stresses of the current particleI can be expressed as follows, using the stress BIE (2) by combining the constitutive relation as well as the integral type transformation [17]:
(9)σijI(p)=EijklεklI=Eijkl(εklCI-εkl0I)iiiiiiii=∑J=1NL∫ΓJεkl0Jxlτijk*dΓ,p∈ΩI,
where Eijkl stands for the elastic tensor of matrix. As there is no outer boundary in full space, the two boundary integrals in (2) vanish. By making extension of the concept of Eshelby tensor to different particles and using (9) and (5), the constrained strains of the current particle and the corresponding Eshelby tensor can be written, respectively, as
(10)εijCI=εij0I+∑J=1NLCijklJ∫ΓJεmn0Jxnτklm*dΓ=∑J=1NLSijklIJεkl0J,(11)SijklIJ=12(δikδjl+δilδjk)δIJiiiiiiiii+14μ∫ΓJxl{τijk*+τjik*-2ν1+νδijτmmk*}dΓ.If the superscriptsI≠J, the Eshelby tensor SijklIJ correlates the eigenstrains and the constrained strains of different particles in the near-field group. Similar to the relation of replacement of (6), according to Hooke’s law, the relation of replacement by equivalent inclusions should be followed for the particles in the near-field group:
(12)(1-β1)εijCI+β2δijεkkCI-εij0I-ν1-2νδijεkk0I=-(1-β1)εijI-β2δijεkkI,I=1,2,…,NL.Combining (10) and (11) and after discretization, the matrix form of (12) can be written as
(13)[S110S12⋯S1,NLS21S220S2,NL⋮⋱⋮SNL,1SNL,2⋯SNL,NL0]{ε10ε20⋮εNL0}={ε1ε2⋮εNL}.In this way, with (13), the eigenstrain vector can be computed using the applied strain vector for the near-field particles. Similar to the analysis for the multiple cracks [8], the matrix in (13) is named Eshelby matrix, which can be looked at as an extension of the concept of Eshelby tensor and the equivalent inclusion in the numerical applications. It needs to be pointed out that (9)–(12) are analytical equations, which reflect special stress-strain correlations among particles of the near-field group in full space, describing their state precisely. Nevertheless, (13) is a numerical form after discretization, describing the state and correlations of the near-field particles approximately.For a problem to be solved, once the radius is given for the dashed line circle as shown in Figure1, each particle in the solid will have, in general case, a unique near-field group and a unique Eshelby matrix correspondingly. For the convenience of computing, rewrite (13) as
(14)ε(I)0=SIε(I),I=1,2,…,NI,
where the matrix SI in (14) is obtained by the inversion then reduction of Eshelby matrix, which corresponds to the eigenstrain vector, ε(I)0, of the current particle. ε(I)={ε1,ε2,…,εNL}T is the applied strain vector of the near-field particles, being the same as that of the right hand side of (13), corresponding to the current particle. The subscript I in (14) goes through all of the particles in the solid; the total number of them being denoted as  NI.In this way, based on the proposed algorithm of the eigenstrain BIE with Eshelby matrix, the computation of eigenstrains of each particle in solids consists of three parts: the first is computed from the stresses or the applied strains resulted from the load. The second comes from the effects of particles in the near-field group, having relatively strong interactions, which is computed using (14). The third is from the effects of particles in the far-field group, having relatively weak interactions, which is computed directly using the domain integrals of the BIE. It is obvious that the minimum number of particles in the near-field group is one. If so (NL= 1), the proposed algorithm would be reduced to that in the previous work [15–17], while (14) becomes (8).


## 3. Subdomain BEM in Full Space

The main purpose of the present work is to verify the correctness of the proposed computational model and the feasibility of the algorithm, so that the results from the subdomain BEM are required as the control for comparison. The two-dimensional analysis has been carried out for two particles in full plate in [18]. As there is no outer boundary in the full space, however, the loading condition needs to be dealt with in a special way. Therefore, a brief introduction is necessary for the subdomain BEM in full space.A single particle denoted asΩI in full space is considered with its boundary or interface ΓI under far-field uniform load as shown in Figure 2(a). Decomposing Figure 2(a) into an interior problem (Figure 2(b)) and an exterior problem (Figure 2(c)) while the interface condition is being kept unaltered, there are the displacement uI and the traction τI on the boundary of the particle in Figure 2(b) without far-field load. In addition to the displacement uI and the traction -τI on the boundary of the hole in Figure 2(c), however, there is far-field load. It is noted that the interface condition should be satisfied between the interior and exterior problems (Figures 2(b) and 2(c)); that is, the continuation in displacement and the equilibrium in traction are held true on the interface so that the signs of tractions are in opposite on the two boundaries. As there is no outer boundary in the full space, in order to express the effect of the far-field load with the boundary integral equations, Figure 2(c) needs to be further decomposed into Figures 2(d) and 2(e), respectively, where the displacement uI0 and the traction -τI0 on the fictitious boundary in dashed line can easily be obtained. The hole in Figure 2(e) is suitable for the description using the boundary integral equations in full space.Schematics of the subdomain BEM in full space.
(a)
(b)
(c)
(d)
(e)For a full space containingNI particles decomposed in such a way, any of the particles (Figure 2(b)) need to be described independently using NI boundary integral equations. All of the holes in full space (there is only one empty hole given in Figure 2(e)) need to be described using another boundary integral equation. In the present work, the so-called subdomain BEM in full space is the discrete form of the equations combined from these above mentioned NI+1 boundary integral equations through the interface conditions. The boundary integral equation is written as
(15)C(p)ui(p)=∫Γτj(q)uij*(p,q)dΓ(q)hhhhhhhhhhh-∫Γuj(q)τij*(p,q)dΓ(q).Notice that there is no domain integral in (15). For any particle shown in Figure 2(b), the matrix form of the boundary integral equation after discretization can be written as
(16)TIuI=UIτI,I=1,2,…,NI,
where UI and TI stand for the coefficient matrices of the displacement and the traction fundamental solutions, respectively. The traction vector on the ΓI can be expressed as follows if the matrix UI is invertible:
(17)τI=UI-1TIuI,I=1,2,…,NI.The matrix form of the boundary integral equation can be written as follows for theNI holes in full space (Figure 2(e)) after discretization:
(18)T(u-u0)=U(τ0-τ),
where  τ={τ1,τ2,…,τNI}T, τ0={τ10,τ20,…,τNI0}T, u={u1,u2,…,uNI}T, u0={u10,u20,…,uNI0}T.Inserting (17) into (18) with arrangement, the following is obtained:(19)[T11+U11U1-1T1T12+U12U2-1T2⋯T1,NI+U1,NIUNI-1TNIT21+U21U1-1T1T22+U22U2-1T2T2,NI+U2,NIUNI-1TNI⋮⋱⋮TNI,1+UNI,1U1-1T1TNI,2+UNI,2U2-1T2⋯TNI,NI+UNI,NIUNI-1TNI]u=Tu0+Uτ0,where UIJ and TIJ are the submatrices in U and T in (18), while the subscripts I and J mean that the source point p and the field point q are placed on ΓI and ΓJ, respectively. When I=J the submatrices are located on the main diagonal of the matrix at the left hand side of (19), describing the same single empty hole itself. In fact, the following relations hold according to the properties of integral kernels:
(20)TI=-TII,(21)UI=UIIβI.The premise of (21) is the assumption that the particle and matrix have the same Poisson’s ratio, where βI=EI/EM is the Young’s modulus ratio between the particle and the matrix. After the displacementsu on the interfaces are obtained using (19), the tractions τI on each interface of the particles can be computed step by step using (17).


## 4. Numerical Examples

As stated above, one of the main purposes of the present work is to check the proposed computational model, that is, the correctness of the Eshelby matrix; only a small number of ellipsoidal particles in full space are taken into consideration in the numerical computations with an assumption that the eigenstrain distributes uniformly in particles. Only one near-field group is chosen; the number of them being the same as that of the total number of particles (NL=NI) in full space. Therefore, there is no iteration required in the numerical examples in the present work. Otherwise, if the previous algorithms are to be employed [15–17], that is, NL= 1, the iteration is indispensable. The definition of the axes of the ellipsoid particle is shown in Figure 3(a). The boundary elements in one octant of the ellipsoid particle are shown in Figure 3(b), which are used for the computation of Eshelby tensors with the boundary type integrations [17] and for the computation of Eshelby matrices in the eigenstrain BIE or for the discrete elements in the subdomain BEM.The definition of the axes of ellipsoid (a) and the boundary elements in one octant (b).
(a)
(b)


### 4.1. Stress Distributions of Two Particles

The two horizontally arranged oblate ellipsoidal particles (a=b=3c; see Figure 3(a)) are shown in Figure 4(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the two centers of particles are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It can be seen from Figure 4(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen also from Figure 4(b) that the stress components σ33 have jumps on the interfaces, no matter how hard the particles are.Two horizontally placed oblate ellipsoidal particles (a) and the stress distributions (b).
(a)
(b)The two horizontally arranged spheroidal particles (a=b=c; see Figure 3(a)) with different modulus are shown in Figure 5(a) under a unit uniform load in z direction in full space, where the soft particle is placed in the left while the hard particle is in the right. The stress distributions along the line connecting the two centers of particles are compared in Figure 5(b), showing also that the results of the two algorithms are in good agreement.Two horizontally placed spheroidal particles with different modulus (a) and the stress distributions (b).
(a)
(b)The two vertically arranged spheroidal particles (a=b=c; see Figure 3(a)) are shown in Figure 6(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the two centers of particles are compared in Figure 6(b), showing also that the results of the two algorithms are in good agreement. It can be seen from Figure 6(b) that the stress components σ11 have jumps on the interfaces in the case of vertical arrangement for both of the soft and hard particles.Two vertically placed spheroidal particles (a) and the stress distributions (b).
(a)
(b)


### 4.2. Stress Distributions of Three Particles

The three horizontally arranged prolate ellipsoidal particles (3a=3b=c; see Figure 3(a)) are shown in Figure 7(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It can be seen from Figure 7(b) that the two results are also in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from Figure 7(b) that the stress components σ33 have jumps on the interfaces of either the soft particle or of the hard particle. It can be seen from Figure 7(b) by careful observations that there are slight differences of the stresses in the two prolate particles since the states of the middle particle and the right particle have some differences.Three horizontally placed prolate ellipsoidal particles (a) and the stress distributions (b).
(a)
(b)The three horizontally arranged spheroidal particles (a=b=c; see Figure 3(a)) are shown in Figure 8(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It is seen also from Figure 8(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It is seen from the comparison between the results of Figures 8(b) and 7(b) the effect of the geometrical shape and the hardness of particles on the stress distributions. For hard particles, the stress components σ33 in the particle are higher than those of prolate particles. In contrary, the stress components σ33 in the matrix are lower. However, the stress components σ33 in the matrix are somewhat higher for soft particles.Three horizontally placed spheroidal particles (a) and the stress distributions (b).
(a)
(b)


### 4.3. Stress Distributions of Five Particles

The five spheroidal particles (a=b=c; see Figure 3(a)) arranged in a vertical plane are shown in Figure 8(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It is seen also from Figure 9(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from Figure 9(b) that the stress components σ33 have jumps on the interfaces of either the soft particle or of the hard particle. As there are three particles, each of them placed on the left, above, and below, the state of the middle particle is different from that of the right particle. It can be seen from Figure 9 that the difference of stresses within the two particles is bigger than that of corresponding values in Figure 8.The horizontal computing line for five spheroidal particles (a) and the stress distributions (b).
(a)
(b)As shown in Figure10 for the five particles, the stress distributions along the line connecting the centers of the middle and the upper particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM. It is seen also from Figure 10(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from the comparison between the results of Figures 10 and 6 the effect of the particle arrangement on the stress distributions.The vertical computing line for five spheroidal particles (a) and the stress distributions (b).
(a)
(b)As is well known, solids containing particles belong to inhomogeneous problems, where some of the stress components on interfaces are continuous while others are discontinuous; that is, there are jumps. In the use of the subdomain BEM, the stresses on the interfaces need to be carried out on both sides of the boundaries as in this case the interfaces correspond to the boundaries of subdomains. However, in the use of the proposed eigenstrain BIE, the stresses on the interfaces are computed directly on the interfaces as in this case the domain in computation becomes apparently homogeneous owing to the replacement of equivalent inclusions. It is shown from the computed results that in the cases of jumps, the stress components computed on the interfaces take the mean values of those of the two sides using the proposed algorithm of the eigenstrain BIE, for example,σ33 in Figures 4, 5, 7, 8, and 9 and σ11 and σ22 in Figures 6 and 10, the same as those in the previous work [17].


### 4.4. Computational Efficiency

The computational efficiencies have big differences between the two algorithms, the eigenstrain BIE and the subdomain BEM in full space. The mean CPU times are compared in Table1 for the two algorithms, showing that the differences of CPU times are as high as two orders of magnitude. The computational efficiency of the proposed algorithm of the eigenstrain BIE is very high. It can be seen from the ratio of CPU times of the two algorithms that the differences of computational efficiencies between the two algorithms will grow with the increase of the number of total particles. This is because the CPU time of the subdomain BEM depends mainly on the solution of the total system matrix. The scale of the total system matrix will grow with the increase of the number of total particles in solids so that the CPU time increases geometrically. In sharp contrast, however, in the algorithm of the eigenstrain BIE, the CPU time increases arithmetically once the scale of the near-field particles is given.Table 1
The mean CPU times of two algorithms.
N
I
2
3
5
BEM
190
422
1265
Eigen
1.69
2.38
6.27
CPU ratio
112
177
202


## 4.1. Stress Distributions of Two Particles

The two horizontally arranged oblate ellipsoidal particles (a=b=3c; see Figure 3(a)) are shown in Figure 4(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the two centers of particles are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It can be seen from Figure 4(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen also from Figure 4(b) that the stress components σ33 have jumps on the interfaces, no matter how hard the particles are.Two horizontally placed oblate ellipsoidal particles (a) and the stress distributions (b).
(a)
(b)The two horizontally arranged spheroidal particles (a=b=c; see Figure 3(a)) with different modulus are shown in Figure 5(a) under a unit uniform load in z direction in full space, where the soft particle is placed in the left while the hard particle is in the right. The stress distributions along the line connecting the two centers of particles are compared in Figure 5(b), showing also that the results of the two algorithms are in good agreement.Two horizontally placed spheroidal particles with different modulus (a) and the stress distributions (b).
(a)
(b)The two vertically arranged spheroidal particles (a=b=c; see Figure 3(a)) are shown in Figure 6(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the two centers of particles are compared in Figure 6(b), showing also that the results of the two algorithms are in good agreement. It can be seen from Figure 6(b) that the stress components σ11 have jumps on the interfaces in the case of vertical arrangement for both of the soft and hard particles.Two vertically placed spheroidal particles (a) and the stress distributions (b).
(a)
(b)


## 4.2. Stress Distributions of Three Particles

The three horizontally arranged prolate ellipsoidal particles (3a=3b=c; see Figure 3(a)) are shown in Figure 7(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It can be seen from Figure 7(b) that the two results are also in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from Figure 7(b) that the stress components σ33 have jumps on the interfaces of either the soft particle or of the hard particle. It can be seen from Figure 7(b) by careful observations that there are slight differences of the stresses in the two prolate particles since the states of the middle particle and the right particle have some differences.Three horizontally placed prolate ellipsoidal particles (a) and the stress distributions (b).
(a)
(b)The three horizontally arranged spheroidal particles (a=b=c; see Figure 3(a)) are shown in Figure 8(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It is seen also from Figure 8(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It is seen from the comparison between the results of Figures 8(b) and 7(b) the effect of the geometrical shape and the hardness of particles on the stress distributions. For hard particles, the stress components σ33 in the particle are higher than those of prolate particles. In contrary, the stress components σ33 in the matrix are lower. However, the stress components σ33 in the matrix are somewhat higher for soft particles.Three horizontally placed spheroidal particles (a) and the stress distributions (b).
(a)
(b)


## 4.3. Stress Distributions of Five Particles

The five spheroidal particles (a=b=c; see Figure 3(a)) arranged in a vertical plane are shown in Figure 8(a) under a unit uniform load in z direction in full space. The stress distributions along the line connecting the centers of the middle and the right particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM, where the modulus ratio for soft particles is EI/EM=0.01 and that for hard particles is EI/EM=10. It is seen also from Figure 9(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from Figure 9(b) that the stress components σ33 have jumps on the interfaces of either the soft particle or of the hard particle. As there are three particles, each of them placed on the left, above, and below, the state of the middle particle is different from that of the right particle. It can be seen from Figure 9 that the difference of stresses within the two particles is bigger than that of corresponding values in Figure 8.The horizontal computing line for five spheroidal particles (a) and the stress distributions (b).
(a)
(b)As shown in Figure10 for the five particles, the stress distributions along the line connecting the centers of the middle and the upper particle are computed using, respectively, the proposed algorithm of the eigenstrain BIE with Eshelby matrices as well as the subdomain BEM. It is seen also from Figure 10(b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices and the feasibility of the proposed algorithm. It can be seen from the comparison between the results of Figures 10 and 6 the effect of the particle arrangement on the stress distributions.The vertical computing line for five spheroidal particles (a) and the stress distributions (b).
(a)
(b)As is well known, solids containing particles belong to inhomogeneous problems, where some of the stress components on interfaces are continuous while others are discontinuous; that is, there are jumps. In the use of the subdomain BEM, the stresses on the interfaces need to be carried out on both sides of the boundaries as in this case the interfaces correspond to the boundaries of subdomains. However, in the use of the proposed eigenstrain BIE, the stresses on the interfaces are computed directly on the interfaces as in this case the domain in computation becomes apparently homogeneous owing to the replacement of equivalent inclusions. It is shown from the computed results that in the cases of jumps, the stress components computed on the interfaces take the mean values of those of the two sides using the proposed algorithm of the eigenstrain BIE, for example,σ33 in Figures 4, 5, 7, 8, and 9 and σ11 and σ22 in Figures 6 and 10, the same as those in the previous work [17].


## 4.4. Computational Efficiency

The computational efficiencies have big differences between the two algorithms, the eigenstrain BIE and the subdomain BEM in full space. The mean CPU times are compared in Table1 for the two algorithms, showing that the differences of CPU times are as high as two orders of magnitude. The computational efficiency of the proposed algorithm of the eigenstrain BIE is very high. It can be seen from the ratio of CPU times of the two algorithms that the differences of computational efficiencies between the two algorithms will grow with the increase of the number of total particles. This is because the CPU time of the subdomain BEM depends mainly on the solution of the total system matrix. The scale of the total system matrix will grow with the increase of the number of total particles in solids so that the CPU time increases geometrically. In sharp contrast, however, in the algorithm of the eigenstrain BIE, the CPU time increases arithmetically once the scale of the near-field particles is given.Table 1
The mean CPU times of two algorithms.
N
I
2
3
5
BEM
190
422
1265
Eigen
1.69
2.38
6.27
CPU ratio
112
177
202


## 5. Conclusion and Expectation

The concept of local Eshelby matrix has been introduced into the computational model of the eigenstrain boundary integral equation (BIE) in order to solve the problem of interactions among particles for solving the large scale numerical simulation of particle reinforced materials. The local Eshelby matrix can be considered as an extension of the concepts of Eshelby tensor and the equivalent inclusion in numerical form. With the introduction of the local Eshelby matrix into the eigenstrain BIE, the original algorithm has been improved and the interaction among the near-field particles has been overcome, which will affect the convergence of the iteration procedures. The three-dimensional stress analyses are carried out for some ellipsoidal particles in full space with the proposed computational model, taking the subdomain BEM as the control. Through the numerical examples, it is verified not only the correctness and feasibility but also the high efficiency of the present algorithm, showing the potential of solving the problem of large scale numerical simulation of particle reinforced materials.As a preliminary work, the eigenstrains are assumed to be constants in particles, which is a limitation of the present work. In the future work, the constant assumption of eigenstrains should be avoided by some appropriate interpolations for the eigenstrains so that the proposed algorithm can be extended to the analysis of particles with arbitrary geometrical shapes.



---

*Source: 947205-2014-04-07.xml*