Extending Esh3D Code to Solve Interacting Eshelby's Inhomogeneity Problems

I present a new development of an open source code Esh3D in solving interacting Eshelby's inhomogeneity problems in whole, half, and finite spaces. The code resolves the interaction between multiple inhomogeneity bodies by assembling and solving a global linear system. For half or finite space problems, the code imposes Neumann and Dirichlet boundary conditions and resolves the surface‐body interaction by hybridizing Eshelby's analytical solutions with a numerical part. I demonstrate the impact of multibody interaction and surface‐body interaction by comparing the results of different problems in whole and finite spaces. This analytical‐numerical hybrid approach, compared to conventional finite element method, is inexpensive and flexible in modeling multiple and evolving inhomogeneity bodies.


Introduction
The Eshelby's (1957Eshelby's ( , 1959Eshelby's ( , 1961 inclusion problem describes the perturbed elastic field in an infinite homogeneous elastic body as a result of the presence of an ellipsoidal inclusion, subject to stress free (eigen)strain. The (constrained) interior strain is uniform and linearly related with the eigenstrain by where * kl is the eigenstrain and S ijkl is named Eshelby tensor. Mura (1987, pp. 22.8-22.13) provides the formulation for the exterior elastic strain Unlike the constant Eshelby tensor, D ijkl (x) is a function of the evaluation location x relative to the ellipsoid center.
The stress perturbation is given by the product of the elastic strain and stiffness tensors For the interior stress, the eigen (inelastic) strain needs to be subtracted from the total strain. Meng et al. (2012) modify (Mura, 1987, pp. 11.30) to obtain the displacement field where is Poisson's ratio and and are given by integrals where Ω is the inclusion region. Operator (·) ,i denotes spatial gradient in i-th dimension.

10.1029/2019EA000594
Eshelby solution has wide application in earth science. Rudnicki (1977Rudnicki ( , 1999 uses the solution to investigate the inception of faulting and to calculate the alteration of regional stress by inhomogeneities. Healy et al. (2006) use the solution to model the evolution of polymodal fracture patterns. Eidelman and Reches (1992) use the solution to interpret the tectonic stress state from fracture patterns in pebbles embedded in a conglomerate. Sharma and Ganti (2004) modify the original solution to simulate nano-inclusions. Katsman (2010) uses the solution to model compaction bands. Jiang (2007Jiang ( , 2014Jiang ( , 2016) extends Eshelby's formalism to allow viscous and anisotropic materials and reviews the efforts that make the solution more explicit, among which Ju and Sun (1999) provide an alternative expression for the exterior Eshelby tensor.

Esh3D Code and Its Recent Development
To evaluate Eshelby solution, in general, one has to apply numerical method to evaluate all the integral components. Reches (1998) and Healy (2009) present closed form solution with computer codes, where two of the ellipsoid's semiaxes are equal.
A computer code that evaluates the triaxial Eshelby solution has not appeared in the public domain until Meng et al. (2012) present their code Esh3D (originally in Matlab™) and make it open source. Esh3D has been used in a number of researches. Meng and Pollard (2014) use it to study deformation band formation in sand bodies. Meng et al. (2013) use it to study elastic field around mixed mode fracture tip. Townsend et al. (2015) use it to study elastic perturbation due to volcanic dike intrusion. Bedayat and Dahi Taleghani (2015) use it to study poroelastic deformations. Guido et al. (2015) and Niu et al. (2017) use it to study ground deformation above depleting/inflating reservoirs. Carter et al. (2018) use it to study fatigue crack deflections.
Later, Meng (2019) translates the original code to Fortran 90 and enables half-space capability for inclusion problems by superposing the whole space solution with numerical solution. In this study, I extend Esh3D to obtain finite space solution, where the surface can interact with the inclusions for inhomogeneity problems.
The numerical code is derived from an open source finite element code Defmod (https://bitbucket.org/stali/ defmod) introduced by Ali (2014) and benchmarked by Meng (2017). To distinguish the whole space (semianalytical) code from this numerical code, I simply call them the "analytical part" and "numerical part," respectively.
The new Esh3D evaluates incomplete elliptical integrals with the Fukushima's (2011) routine. The code includes a custom (Okada, 1985) routine to allow multiple Eshelby's inclusions and Okada's faults in a single model. The code is openly available online (https://github.com/Chunfang/Esh3D).

Equivalent Inclusion Method for Inhomogeneity Problems
Perturbations due to an ellipsoidal inhomogeneity subject to uniform remote strain ∞ can be treated as an equivalent inclusion problem. The effective (fictitious) eigenstrain is given by (Eshelby, 1957;Mura, 1987) (C − ΔCS) * * = ΔC ∞ + C in * , where C and C in are the stiffness matrices of the host and inclusion, respectively; ΔC = C − C in ; and * is an arbitrary prescribed (intrinsic) eigenstrain.
The total exterior strain is the remote strain plus the perturbation strain: When the inclusion has the same stiffness matrix as the host, that is, ΔC = 0, ** is equal to * . Therefore, the original inclusion problem is a special case of inhomogeneity problem with prescribed eigenstrain.
In Esh3D, I form a unified approach that allows coexistence of three types of inclusions: (1) pure inclusion with prescribed eigenstrain, (2) pure inhomogeneity with elastic moduli different than the host, and (3) inhomogeneity inclusion with both different moduli and prescribed eigenstrain. In the rest of this paper, I simply use the word "inclusions," without specifying the types.
Unlike the original inclusion problem, the inhomogeneity problem requires the eigenstrain to be resolved from the remote strain and differential elastic moduli. In presence of multiple inclusions, the effective host strain around one inclusion is perturbed by the neighbor inclusions and is in general different than the remote strain. In this study I formulate a global linear system to resolve such multibody interaction.

Linear System Resolving Multibody Interaction
When there are n inclusions, the effective host strain around ith inclusion can be substantially different than the remote strain, due to the presence of neighboring n − 1 inclusions. I modify the equivalent eigenstrain formulation in equation (6) as where (·) i denotes the same quantity as in equation (6) for i-th inclusion; D ij denotes the exterior Eshelby's tensor of j-th inclusion evaluated at the centroid of i-th inclusion. The summation term adds the interacting perturbation to the remote strain.
Equation (8) treats the eigenstrains as piecewise uniform within each inclusion while varying across different inclusions. This formulation approximates the effective host stress that considers the influence of neighboring inclusions. Novák et al. (2012) suggest that this approximation is accurate, even if the inclusions are closely spaced, and propose an iterative "self-compatibility" algorithm to resolve these coupled eigenstrains. This iterative approach reevaluates the host stress and perturbs the eigenstrain of individual inclusions to seek a set of converging eigenstrains. In this study, I rearrange equation (8) to have all the **(·) terms on one side, forming a global linear system: where all the 6 × 6 matrices, Ss and Ds, need to be evaluated in the same global coordinate, before being assembled.
For a small number of inclusions, for example, n < 10 5 , one can solve equation (9) efficiently. For a large number, one can sparsify the coefficient matrix by limiting the interaction range between two inclusions, which will accelerate any solving techniques.

Numerically Imposing Neumann and Dirichlet Boundary Conditions
I generalize the method of obtaining half-space Eshelby solution to obtain finite space solution, where traction (Neumann) boundary conditions are imposed by traction cancellation and displacement (Dirichlet) boundary conditions are imposed by displacement cancellation. For inhomogeneity problems, I consider the boundary-inclusion interaction by adjusting ∞ in equation (9) to where surfi is the strain perturbation at the centroid of ith inclusion due to the presence of half or finite space surface.
The numerical routine for finite space solution is as follows: 1. Solve the linear system of equation (9) for interacting eigenstrains and evaluate stress at the centroids of Neumann surface elements and displacement at the vertices of Dirichlet surface elements. 2. Estimate the nodal forces that cancel unwanted traction at the Neumann surface and nodal displacements that cancel unwanted displacement at the Dirichlet surface. 3. Place the canceling nodal forces and displacements as Neumann and Dirichlet boundary conditions to the numerical domain, respectively, and solve this boundary condition problem. 4. Superpose the tractions and displacements at the Neumann and Dirichlet surfaces, respectively, and the strains at the centroids of the inclusions, that is, surfi in equation (10). If | surfi | max < tol, where tol is an arbitrarily small value, check out with the latest superposed solution. 5. If not, adjust equation (9) by equation (10), add the opposite residual traction and displacement as Neumann and Dirichlet boundary conditions to the numerical domain, respectively, and repeat this routine from the beginning. This approach is functionally similar to the hybrid-Treffz approach by Novák et al. (2012, section 4), except the numerical part here is optional. If a problem can be idealized as a whole space problem, one can solve equation (9) and evaluate the solution directly, without engaging this iterative numerical correction. Therefore, a whole space problem does not require any numerical grid. Zheng et al. (2016) apply similar formalism to waveform modeling in layered 2-D medium that contains multiple elliptical voids.

Example Problems
I compare the results for five inhomogeneity problems: 1. Whole space containing closely spaced inhomogeneity ellipsoids under uniaxial loading. The multibody interaction is captured by directly solving equation (9). 2. Same as (1), except the interaction is neglected. This is obtained by setting the off-diagonal submatrices in equation (9)   The comparison between the first two problems demonstrates the effect of inclusion-inclusion interaction. The comparison between the last three problems demonstrates the effect of surface-inclusion interaction.
I omit the half-space problems, but the technique to obtain finite space solution is equally applicable to a half-space problem, where only the top (Neumann) surface will interact with the inclusions. Table 1 lists all the ellipsoidal inclusions involved in these problems, where the first five inclusions appear in the Problems 1 and 2, and the last three, same as in Novák et al. (2012, table 2), appear in the Problems 3-5.
In all the problems, the host matrix has Young's modulus 80 GPa and Poisson's ratio 0.25, while all the inclusions have twice the Young's modulus and the same Poisson's ratio. Figure 1 illustrates the problem layout. I place uniaxial remote loading yy = 10 7 Pa to excite inhomogeneity perturbation which is evaluated along a line 15 km above the inclusion cluster. I set this evaluation line such that the perturbation is sensitive to multibody interaction and is not dominated by a single inclusion. inhomogeneities. The displacement, u (·) , and stress, (··) , are evaluated along a line 15 km above the inclusions (see Figure 1). Figure 2 compares the displacement and stress perturbations of interacting and noninteracting inhomogeneities, where each row, that is, displacement u i , principal stress ii , and shear stress ij , is set to the same scale. Both the solutions are obtained by solving equation (9), except the noninteracting case has the off-diagonal submatrices set to zeros.

Inhomogeneities in Whole Space
The outstanding differences in u z and three principal stress components suggest that the multibody interaction has significant impact to the displacement and stress perturbations. Due to the dipping angles of the second and fourth inclusions, see Table 1, displacement u x and shear stress xz are slightly above zero. Figures 3a and 3b illustrate the hexahedron and tetrahedron grids used by the hybrid method and pure finite element method, respectively. Ideally, I should also use hexahedron grid in the purely numerical case, for better precision and equivalence. But the presence of the inclusions makes it difficult to discretize the domain with conforming hexahedra. This "hex versus tet" difference may contribute to the disagreement between the hybrid and purely numerical results, which is discussed later. Figure 3c illustrates The Dirichlet (roller) boundary conditions that prevent movements in normal directions but allow movements in tangential directions. On the opposite sides of the roller boundaries, I place isotropic compression traction, xx = yy = zz = −1 MPa, to excite inhomogeneity perturbation.

Inhomogeneities in Finite Space
I set the evaluation line, y = z = 0 in Figure 3d, such that it intersects with two of the three inclusions. Figure 4 compares the results produced by the hybrid method, purely analytical method (for full space), and purely numerical, finite element, method, in terms of displacement and stress perturbations.
As expected, the finite space correction makes the hybrid result agree better with the numerical result. The hybrid result has the middle portion of the vertical displacement, u z , slightly overshooting the numerical (b) tetrahedron grid of 2,328,524 cells used by the "pure" finite element method, where both the ellipsoids and host are discretized; (c) traction (Neumann) boundary condition at x + , y + , and z + sides and roller (Dirichlet) boundary condition, that is, zero normal displacement, at x − , y − , and z − sides; (d) cross section of (b) exposing the evaluation line, y = z = 0, along which displacement and stress by different methods are compared.
result. The stress plateaus, resulting from the intersection between the evaluation line and two inclusions (see Figure 3d), are also slightly different.
The possible causes for these discrepancies include the following: • Equation (8) approximates the eigenstrains as piecewise uniform within each inclusion. In this problem, the close spacing of the inclusions may slightly affect the accuracy of this approximation. • The numerical (finite element) method describes the ellipsoidal inclusions with small tetrahedra. The smooth shapes of the ellipsoids cannot be preserved exactly. This may explain the gaps between the stress plateaus.
Given their fundamental differences, I consider the hybrid and numerical methods reasonably agree.
Interestingly, the full space result agrees with the other two in most of the stress components, except it shows different trends in zz , when approaching the x + (Neumann) surface. Note that the full space case still considers the inclusion-inclusion interaction by keeping all the entries of the matrix in equation (9) while neglecting the surface-inclusion interaction.

Discussion and Conclusion
The advantages of this hybrid method, over purely numerical method, are flexibility and efficiency. One can use the same numerical grid for different and evolving inclusion settings. The numerical part of the hybrid method requires considerably less number of unknowns than its purely numerical counterpart, where the inclusions need to be discretized by much finer grid than the host (see Figure 3d). Figure 5 compares the computational cost in terms of degrees of freedom (DoFs) required by the hybrid and numerical methods for the finite space problem in section 4.2, when gradually increasing the number of inclusions. When there is no inclusion, the hexahedron (hybrid) grid, shown in Figure 3a, has even more DoFs than an equivalent tetrahedron (numerical) grid. As soon as I add one of the inclusions, the DoFs of the tetrahedron grid increase dramatically. When I add all the three inclusions, the total number of numerical DoFs becomes more than one order larger than in the homogeneous case. The additional number of DoFs required by the hybrid method is only 6n (n being the total number of inclusions) in form of the unknowns in equation (9).
When dealing with evolving inclusions, the numerical method requires meshing the domain repetitively.
The hybrid method, in contrast, only requires to update corresponding ellipsoidal parameters, that is, the columns in Table 1.
The distinctive features of this hybrid method, compared to the hybrid-Treffz method (Novák et al., 2012, section 4), include the following: • One can choose part of the domain surface to interact with the inclusion bodies. To solve a half-space problem, for example, only the top (topographical) surface needs to interact with the inclusions. To solve a whole space problem, there is no need to engage the numerical part at all, which makes the model preparation and execution trivially inexpensive. • The loading source can be partially or totally implemented by the analytical part, instead of solely via traction (Neumann) boundary conditions of the numerical part. The Example Problems (1), (2), and (4) have the loading sources prescribed via remote stresses. • Also, an inclusion can behave as a loading source by having "intrinsic eigenstrain," that is, *(·) in equation (8). At the same time, the source inclusion can also behave as an inhomogeneity body that interacts with others. • Alternatively, one can use Okada's (1985) fault slip as loading mechanism. In this case, the fault slip will excite perturbation in the presence of inhomogeneities, but the perturbation will not in return affect the fault slip. This is because the Okada's fault, as apposed to Eshelby's inclusion, does not consider source-host compatibility (Meng, 2019).
In conclusion, I extend the Esh3D code to solve interacting Eshelby's inhomogeneity problems in whole, half, and finite spaces. The new functionalities improve the code's applicability to geophysical problems.