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.
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 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 . 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 and the eigenstrain of a single particle in full space are correlated with Eshelby tensor, ; that is,
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
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 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 in (14) is obtained by the inversion then reduction of Eshelby matrix, which corresponds to the eigenstrain vector, , of the current particle. 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–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 in full space is considered with its boundary or interface 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 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 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 and the traction 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.
(a)
(b)
(c)
(d)
(e)
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 so-called subdomain BEM in full space is the discrete form of the equations combined from these above mentioned 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 and 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 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 and are the submatrices in and 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).
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 () 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, 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.
(a)
(b)
4.1. Stress Distributions of Two Particles
The two horizontally arranged oblate ellipsoidal particles (; see Figure 3(a)) are shown in Figure 4(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 and that for hard particles is . 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 have jumps on the interfaces, no matter how hard the particles are.
(a)
(b)
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.
(a)
(b)
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 have jumps on the interfaces in the case of vertical arrangement for both of the soft and hard particles.
(a)
(b)
4.2. Stress Distributions of Three Particles
The three horizontally arranged prolate ellipsoidal particles (; 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 and that for hard particles is . 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 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.
(a)
(b)
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 and that for hard particles is . 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 in the particle are higher than those of prolate particles. In contrary, the stress components in the matrix are lower. However, the stress components in the matrix are somewhat higher for soft particles.
(a)
(b)
4.3. Stress Distributions of Five Particles
The five spheroidal particles (; see Figure 3(a)) arranged in a vertical plane 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 and that for hard particles is . 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 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.
(a)
(b)
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.
(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, in Figures 4, 5, 7, 8, and 9 and and 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 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.
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work was supported by National Natural Science Foundation of China (nos. 11272195, 11332005, and 10972131).