Eigenstrain boundary integral equations with local Eshelby matrix for stress analysis of ellipsoidal particles

: 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 sub-domain 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. ____________________________________________________________


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][4][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][11][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][16][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, 2 Mathematical Problems in Engineering 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.

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 Ω  is the domain of particle with the boundary or interface Γ  (Γ  = Ω  ∩ Ω).The displacement and the traction eigenstrain boundary integral equations are, respectively, as follows: where represents the free term resulted from the domain integral of (2).Ω  with the boundary Γ  is the infinitesimal region  in Ω  .  =   () −   () stands for the projection of the distance between the field and source points,  and .In (1) and ( 2),  *  ,  *  , and  *  are the Kelvin fundamental solutions of the displacement, the traction, and the stress, respectively. *  ,  *  , and  *  are the corresponding derivations of the fundamental solutions.  is the total particle number in Ω.  0  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.[1], the tress-free constrained strain    and the eigenstrain  0  of a single particle in full space are correlated with Eshelby tensor,   ; that is,

Eshelby Tensor and Equivalent Inclusions. According to the work of Eshelby
The value of Eshelby tensor depends on the shape of Ω  .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: 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 Ω  .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  =   /  , where the subscripts  and  denote the particle and matrix, respectively.According to Hooke's law, if a particle under the applied strain   is replaced by an equivalent inclusion without altering its stress state, the following relation should be satisfied: where 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

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][16][17].As shown in Figure 1, the group of the near-field particles having a number   is defined as those within the circle of dashed line for the current particle  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.
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 particle  can be expressed as follows, using the stress BIE (2) by combining the constitutive relation as well as the integral type transformation [17]: where   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 If the superscripts  ̸ = , the Eshelby tensor    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: Combining ( 10) and ( 11) and after discretization, the matrix form of ( 12) can be written as 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 Figure 1, 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 where the matrix S  in ( 14) is obtained by the inversion then reduction of Eshelby matrix, which corresponds to the eigenstrain vector,  0 () , of the current particle. () = { 1 ,  2 , . . .,    }  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  in (14) goes through all of the particles in the solid; the total number of them being denoted as   .
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 (  = 1), the proposed algorithm would be reduced to that in the previous work [15][16][17], while ( 14) becomes (8).

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 twodimensional 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 Ω  in full space is considered with its boundary or interface Γ  under far-field uniform load as shown in Figure 2 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  0  and the traction − 0  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.
For a full space containing   particles decomposed in such a way, any of the particles (Figure 2(b)) need to be described independently using   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 socalled subdomain BEM in full space is the discrete form of the equations combined from these above mentioned   + 1 boundary integral equations through the interface conditions.The boundary integral equation is written as 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 where U  and T  stand for the coefficient matrices of the displacement and the traction fundamental solutions, respectively.The traction vector on the Γ  can be expressed as follows if the matrix U  is invertible: The matrix form of the boundary integral equation can be written as follows for the   holes in full space (Figure 2(e)) after discretization: where Inserting (17) into (18) with arrangement, the following is obtained: where U  and T  are the submatrices in U and T in (18), while the subscripts  and  mean that the source point  and the field point  are placed on Γ  and Γ  , respectively.When  =  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: The premise of (21) is the assumption that the particle and matrix have the same Poisson's ratio, where   =   /  is the Young's modulus ratio between the particle and the matrix.
After the displacements u on the interfaces are obtained using (19), the tractions   on each interface of the particles can be computed step by step using (17).

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 nearfield group is chosen; the number of them being the same as that of the total number of particles (  =   ) 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][16][17], that is,   = 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 two horizontally arranged spheroidal particles ( =  = ; see Figure 3(a)) with different modulus are shown in Figure 5(a) under a unit uniform load in  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.
The two vertically arranged spheroidal particles ( =  = ; see Figure 3(a)) are shown in Figure 6(a) under a unit uniform load in  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.

Stress Distributions of Three Particles.
The three horizontally arranged prolate ellipsoidal particles (3 = 3 = ; see Figure 3(a)) are shown in Figure 7(a) under a unit uniform load in  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   /  = 0.01 and that for hard particles is   /  = 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.
The three horizontally arranged spheroidal particles ( =  = ; see Figure 3(a)) are shown in Figure 8(a) under a unit uniform load in  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   /  = 0.01 and that for hard particles is   /  = 10.It is seen also from   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.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.As shown in Figure 10 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.
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  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 Table 1 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.

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 can be

IFigure 1 :
Figure 1: The group definitions for multiple particles.

Figure 2 :
Figure 2: Schematics of the subdomain BEM in full space.
(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   and the traction   on the boundary of the particle in Figure 2(b) without far-field load.In addition to the displacement   and the traction −  on the boundary of the hole in Figure

Figure 3 :
Figure 3: The definition of the axes of ellipsoid (a) and the boundary elements in one octant (b).

4. 1 .
Stress Distributions of Two Particles.The two horizontally arranged oblate ellipsoidal particles ( =  = 3; see Figure3(a)) are shown in Figure4(a) under a unit uniform load in  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   /  = 0.01 and that for hard particles is   /  = 10.It can be seen from Figure4(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 Figure4(b) that the stress components  33 have jumps on the interfaces, no matter how hard the particles are.

Figure 4 :
Figure 4: Two horizontally placed oblate ellipsoidal particles (a) and the stress distributions (b).

Figure 5 :
Figure 5: Two horizontally placed spheroidal particles with different modulus (a) and the stress distributions (b).

Figure 8 (
b) that the two results are in good agreement with each other, showing the correctness of the Eshelby matrices

Figure 6 :
Figure 6: Two vertically placed spheroidal particles (a) and the stress distributions (b).

Figure 7 :
Figure 7: Three horizontally placed prolate ellipsoidal particles (a) and the stress distributions (b).

Figure 9 :
Figure 9: The horizontal computing line for five spheroidal particles (a) and the stress distributions (b).

Figure 10 :
Figure 10: The vertical computing line for five spheroidal particles (a) and the stress distributions (b).

Table 1 :
The mean CPU times of two algorithms. 11 and  22 in Figures