A refined numerical solution to the inclusion problem

Eigenstrains is a generic name given by Mura [1] to nonelastic strains such as plastic strains, misfits strains, thermal expansion or phase transformation, which generate a linear elastic stress field in an elastic isotropic media, referred to as the surrounding matrix. It should be noted that the region with eigenstrains, namely the inclusion, is assumed to have the same elastic parameters as the surrounding matrix. Solution for the stresses due to inclusions in an infinite space can be achieved using Eshelby’s method [2]. However, in many practical situations, sites of eigenstrains generated by contact loading, heating or metallurgical transformations are confined to a near surface region. Although inclusion problem has received a great deal of attention in the last four decades [3], closed-form solutions exist only in a few cases of regular shaped inclusions. Solutions of elastic state due to a spherical or ellipsoidal inclusion containing pure dilatational eigenstrains have been obtained by Mindlin and Cheng [4] and by Seo and Mura [5], respectively, through integration of Green’s functions of a point force acting in the interior of a semiinfinite solid [6]. The closed form solution for stresses due to uniform cuboidal eigenstrains in an elastic isotropic infinite space, as well as in the half-space, was advanced by Chiu [7] and by Chiu [8] respectively, using the mirror image method. This solution was later used by Jacq et al. [9] to solve the problem of residual stresses in elasticplastic contact due to an arbitrarily shaped plastic region, approximated by the reunion of a finite number of elementary cuboids, each containing uniform plastic strains. Wang and Keer [10] used a similar approach, implementing the Discrete Convolution Fast Fourier Transform (DCFFT) technique [11] in layers of constant depth. Liu and Wang [12] approached the inclusion problem by superposition, differentiation and integration of Mindlin and Cheng’s [13] solutions for the point force in the interior of the semiinfinite solid, but their formulation led to complex calculations involving derivatives of four key integrals. More recent refinements of the numerical method for the inclusion problem involve a different approach [14] in derivation of influence coefficients, which allows coupling with three-dimensional spectral methods, leading to a further decrease of order of computation. When solving more complex eigenstrains configurations, or when the inclusion problem has to be solved repeatedly, as in the case of elastic-plastic contact solvers, the main challenge consists in correlating the available computational resources with the imposed precision goals. In many practical situations, problems are considered unsolvable due to prohibitive requirement of memory or processor speed. The goal of this paper is to advance a refined numerical method for the inclusion problem, incorporating the latest enhancements in a robust, efficient algorithm, which can be used to solve a wider variety of inclusionrelated problems. The strong point of this method consists in the hybrid convolution-correlation multi-dimensional algorithm, which allows the use of three-dimensional DCFFT. The simplified method for enforcing free surface relief in the mirror image method allows extension to the case of general (not only deviatoric) eigenstrains.


Introduction
Eigenstrains is a generic name given by Mura [1] to nonelastic strains such as plastic strains, misfits strains, thermal expansion or phase transformation, which generate a linear elastic stress field in an elastic isotropic media, referred to as the surrounding matrix.It should be noted that the region with eigenstrains, namely the inclusion, is assumed to have the same elastic parameters as the surrounding matrix.Solution for the stresses due to inclusions in an infinite space can be achieved using Eshelby's method [2].However, in many practical situations, sites of eigenstrains generated by contact loading, heating or metallurgical transformations are confined to a near surface region.Although inclusion problem has received a great deal of attention in the last four decades [3], closed-form solutions exist only in a few cases of regular shaped inclusions.
Solutions of elastic state due to a spherical or ellipsoidal inclusion containing pure dilatational eigenstrains have been obtained by Mindlin and Cheng [4] and by Seo and Mura [5], respectively, through integration of Green's functions of a point force acting in the interior of a semiinfinite solid [6].The closed form solution for stresses due to uniform cuboidal eigenstrains in an elastic isotropic infinite space, as well as in the half-space, was advanced by Chiu [7] and by Chiu [8] respectively, using the mirror image method.This solution was later used by Jacq et al. [9] to solve the problem of residual stresses in elasticplastic contact due to an arbitrarily shaped plastic region, approximated by the reunion of a finite number of elementary cuboids, each containing uniform plastic strains.Wang and Keer [10] used a similar approach, implementing the Discrete Convolution Fast Fourier Transform (DCFFT) technique [11] in layers of constant depth.Liu and Wang [12] approached the inclusion problem by superposition, differentiation and integration of Mindlin and Cheng's [13] solutions for the point force in the interior of the semiinfinite solid, but their formulation led to complex calculations involving derivatives of four key integrals.
More recent refinements of the numerical method for the inclusion problem involve a different approach [14] in derivation of influence coefficients, which allows coupling with three-dimensional spectral methods, leading to a further decrease of order of computation.When solving more complex eigenstrains configurations, or when the inclusion problem has to be solved repeatedly, as in the case of elastic-plastic contact solvers, the main challenge consists in correlating the available computational resources with the imposed precision goals.In many practical situations, problems are considered unsolvable due to prohibitive requirement of memory or processor speed.
The goal of this paper is to advance a refined numerical method for the inclusion problem, incorporating the latest enhancements in a robust, efficient algorithm, which can be used to solve a wider variety of inclusionrelated problems.The strong point of this method consists in the hybrid convolution-correlation multi-dimensional algorithm, which allows the use of three-dimensional DCFFT.The simplified method for enforcing free surface relief in the mirror image method allows extension to the case of general (not only deviatoric) eigenstrains.

Numerical model
Considering existing closed-form solutions, it is convenient to divide the inclusion domain into a collection of non-intersecting cuboids of uniform eigenstrains, using a rectangular three-dimensional mesh, and to approximate stresses due to arbitrarily shaped inclusions by superposing the individual contributions of each elementary cuboid.The starting point in the numerical procedure is the closedform solution advanced by Chiu [7] for the elastic stress field due to a cuboid of uniform eigenstrains located in an infinite, elastic and isotropic medium, referred to as the influence coefficient (IC).The main features of this solution are briefly outlined.A Cartesian coordinate system 1 2 3 ( , , ) x x x is attached to the center of the cuboid contain- ( , , ) x x x (i.e. the observation point), due to the eigenstrains tensor ij   (i.e. the source cuboid), can be ex- pressed as the contribution of all components of the ij   tensor, yielding 36 different ICs:  .
Consequently, A ijkl , i, j, k,  = 1, 2, 3, denotes the stress tensor component ij induced in the infinite elastic matrix by the k component of the eigenstrains tensor, assumed uniform and equal to unity in a cuboid centered in origin.It should be noted that, in case of shear strains, i.e. k  , as AA  , contribution of both corresponding strains, i.e. of both 1 kk    , is accounted for in Eq. (1)   by means of the multiplier 2.

ijk
A is therefore a fourth- order tensor, whose components depend on the size of the source cuboid and on the distance between the center of the cuboid and the observation point.Eq. ( 1) can be written in a simplified form, using Einstein summation convention (which will be used from this point forward): where i, j, k,  = 1, 2, 3. Detailed derivation of ijk A in closed-form expression is discussed in [9] or [14].
As pointed out by Mura [15], in the presence of initial strains, a finite body with a traction-free surface can be treated as an infinitely extended body, if equal but opposite normal and shear stresses are applied on the boundary, compensating for the ones corresponding to the full space solution.Based on this assertion, Chiu [8] derives the elastic field due to a cuboidal inclusion in an elastic half-space, by superposition of three solutions, corresponding to three elastic states, denoted by superscripts I, II and III in Fig. 1.The coordinate system is moved with its origin on the half-space boundary, laying on the normal axis passing through the centre of the cuboid.State I employs the eigenstrains of the original problem, i.e.
but the surrounding elastic matrix is considered infinite in all directions.Let (0, 0, ) h denote the position of the centre of the cuboidal inclusion.The same type of boundary condition is imposed in state II, except the eigenstrains region is chosen as the mirror image of that in state I (i.e. a cuboid having the centre at (0,0, ) h  ), in such a manner that summation of elastic fields induced by both cuboidal regions leaves the half-space boundary free of shear tractions.Using Eq. ( 2) with the new coordinate system, superposition principle applied to states I and II yields the stress induced in any point of the infinite space: ( , , ) ( , , ) (0, 0, ) ( , , ) (0, 0, ), where i, j, k,  = 1, 2, 3.The choice of eigenstrains depicted in Fig. 1  ( , ,0) . According to the method indicated by Mura [15], stresses equal but opposite to ( , ,0) should be applied on the boundary x 3 = 0, and their contribution subtracted from summation of solutions I and II, i.e. from the ( , , ) computed in Eq. ( 3), in order to fully satisfy the traction free surface condition in the original inclusion problem.Let ( , ,0) yielding from the states I and II.The solution of the inclusion problem with cuboidal uniform eigenstrains in a finite body can then be obtained by superposition of the three elastic states, as depicted in Fig. 1: ( , , ) ( , , ) ( , , ) The latter solution, also referred to as the influence coefficient for the inclusion problem in a finite medium, can be further employed to solve numerically any type of inclusion.In the numerical formulation, which involves digitization of problem parameters, the inclusion domain, which can be arbitrary shaped, containing known but otherwise arbitrarily distributed eigenstrains, is divided into elementary cuboids of uniform eigenstrains, as shown in Fig. 2. Superposition principle is subsequently applied to evaluate the resulting stress state.As all distributions are assumed piece-wise constant and equal, in each cuboid, to the value in its centre, it is convenient to index all cuboids in the rectangular computational domain by a sequence of three integers ranging from 1 to N 1 , N 2 and N 3 respectively, ( , ,0) and to express all distributions as functions of these integers instead of continuous coordinates.For example, the value of any continuous distribution f in the centre C ijk of the cuboid (i, j, k) will be denoted by f (i, j, k), and will be computed as f(x 1 (C ijk ), x 2 (C ijk ), x 3 (C ijk )).In this manner, the integral equation yielding stresses due to an arbitrary distribution of eigenstrains, occupying an arbitrary domain, can be estimated numerically by summating the individual contributions of all elementary cuboids with non-vanishing eigenstrains.

Fig. 2 Inclusion problem digitization
However, in order to implement spectral methods, which compute multiple elementary contributions simulta-neously, the multi-summation process must cover all the cells in the computational domain.Applying this technique to Eq. (3) yields: ) ( , , ), where 1 2 3 , , , , , In this relation, the cuboid ( , , ) mn    , of centre mn C , is the mirror image with respect to the plane x 3 = 0 of the cuboid ( , , ) mn centred at mn C .Consequently, Eq. ( 4) expresses the  stress tensor component induced in any observation cuboid ( , , ) i j k by a digitised distribution of eigenstrains located in the half-space 3 0 x  (the first term on the right-hand side of Eq. ( 4)), as well as by its mirror image with respect to the plane 3 0 x  (the second term on the right-hand side of Eq. ( 4)), both laying in an infinite elastic matrix.Due to symmetry, 11 ( ) ( ) . Consequently, the term accounting for the contribution of the mirror images can be expressed using the set of the indexes employed in state I, provided the correct coordinate system transformation is performed, yielding: It is also convenient to match the set of observation points to the set of centres of all elementary cuboids, i.e.  1 kN  .In this manner, the number of different ICs to be computed is restricted to the number of different distances between observation and source nodes (possibly with a changed sign), yielding at most 1 2 3 8N N N different terms for each of the two mem- bers on the right-hand side of Eq. ( 5).
The fictitious normal traction 33 1 2 ( , ,0) acting on the half-space boundary, needed for the solution of state III, yields from Eq. ( 5) by setting k = 1 and 3   .The stress induced by this traction in the half- space 3 0 x  , denoted ( , , ) such as the plastic strains.Moreover, the resulting influence coefficients depend not only on the distance between the observation and the source points, but also on the depth of the source point, leading to a total of  The approach newly proposed in this paper differs from that of Chiu in the way III   is related to    , allow- ing extension of the numerical solution to the case of general (not only deviatoric) eigenstrains, and providing an important reduction of computational complexity and of memory requirements.According to [16], the stress induced in an elastic half-space by a digitised distribution of normal tractions can be approximated by the summation: The influence coefficient S  expresses the  stress ten- sor component induced in the cuboid (i, j, m) by a unit uniform pressure acting on the side of the cuboid (k, , 1) included in the half-space boundary.Closed-form expressions for these ICs were derived in [16].Finally, the solution for the stress due to digitised arbitrary shaped eigenstrains in an elastic isotropic half-space results from superposition of solutions ( 5) and ( 6), as shown in Eq. ( 3) for a single cuboidal inclusion.

Computation acceleration
The multi-summation process in Eq. ( 5) is very computationally intensive, having an order of computation of O(N 2 ) for a grid with N elements when classical direct summation is used.In order to circumvent this constraint, which limits dramatically the resolution that can be imposed when solving numerically the inclusion problem, many authors [9,10,12,14] employed spectral methods, based on the convolution theorem.The main idea is to transfer the convolution calculation from the space domain (SD) to the frequency domain (FD), where the computational complexity is reduced to an improved order of O(N log N).The source of this reduction is the convolution theorem [17], which states that convolution in the SD can be computed as an element-wise product in the FD, in O(N) operations.Additional computational effort is required to apply the direct (FFT) and inverse (IFFT) Fast Fourier Transform, to accomplish the transfer back and forth between the SD and the FD.
In [9], formulas for the global ICs of the inclusion problem were derived by superposition of the ICs for a single cuboid corresponding to states I, II and III, and superposition of all cuboids contributions was subsequently performed.This course of action limits the use of spectral methods to the layer-by-layer two-dimensional case, as the global ICs depend explicitly not only on the distance between the observation and the source, but also on the source depth.It should be noted that, in order to apply DCFFT in the direction of 3 x , the global ICs should de- pend only on the distance between the source and the observation points.Using a two-dimensional algorithm to solve an intrinsically three-dimensional problem is an imperfect solution, which limit dramatically the resolution that can be achieved in the numerical approach.In the newly proposed algorithm, superposition principle with respect to individual contributions of elementary cuboids of the same state is applied first, and summation of solutions I, II and III is subsequently performed.In this manner, the global ICs derived in [9] are no longer computed, as solution of state III is achieved in a manner different from that of Chiu [8], and the new algorithm can benefit from implementation of three-dimensional DCFFT.
The two terms in the right-hand side of Eq. ( 5) involve multi-summation in the three-dimensional space, with both source and observation domains threedimensional.While the first term in Eq. ( 5) is a threedimensional convolution, the second term is a twodimensional convolution with respect to directions of 1 x and of 2 x , and a one-dimensional correlation with respect to direction of 3 x .The DCFFT technique presented in [11]  for the one-dimensional case extends naturally to three dimensions, while for the second term, a special hybrid convolution-correlation algorithm is advanced herein.The main issue is to adapt the DCFFT algorithm, which is very efficient in terms of computational complexity, to allow calculation of correlation products as well.
The base for this new approach is the correlation theorem, which states [17] that a correlation product can be assessed as the spectral (i.e. in the FD) convolution between one member of the correlation and the complex conjugate of the other.Therefore, classical DCFFT technique can be applied with respect to direction of 3 x too, provided the spectral eigenstrains array is replaced by its complex conjugate.This substitution can be achieved by reordering the terms in the eigenstrains array in upturned order with respect to direction of 3 x , prior to transfer to the FD.In- deed, when FFT is applied to a series s of real terms in the SD to acquire its spectral counterpart ŝ , one can obtain the complex conjugate of ŝ , denoted by ŝ , by rearranging the original series s in reversed order before applying the FFT.This feature allows combining convolutions and correlations products with respect to different directions in a hybrid convolution-correlation multi-dimensional algorithm, whose flow-chart is presented in Fig. 3.The algorithm steps are described in the following section.
Firstly, compute the three-dimensional arrays of ICs entering Eq. ( 5).It should be noted that the same operations apply to all the ICs in Eq. ( 5), regardless of the referred stress or strain tensor component, i.e. regardless of the instance of indices , , ,      . Therefore, the subscripts denoting the tensor components will be omitted in this section for brevity, and classical matrix notation will be employed instead, i.e.  ( , , ) i j m element of any of the 36 three-dimensional arrays in Eq. ( 1), relating any stress to any eigenstrains tensor component.
As the algorithm steps match along any direction corresponding to the same type of product (i.e.convolution or correlation), algorithm description will be limited to the one-dimensional case for brevity.The notation will be simplified accordingly, thus let us denote the threedimensional arrays B i, j, m as vectors B  , 12 k N  , having as components two-dimensional arrays (i.e.matrices), where k can be any of the three directions, i.e. k = 1, 2 or 3, and N k is the number of grids in the direction of k x .The algorithm advanced herein can also be applied to any multidimensional configurations, i.e. when k > 3.As stated before, computation of B  components depends weather k corresponds to a convolution or a correlation.In case of convolution, B 1 corresponds to the negative greatest distance between two grid nodes taken along the x k -axis, i.e the observation mark is indexed with 1 and the source mark with N k .The situation when the observation and the source points are interchanged corresponds to 2 k N B .Alt- hough there are apparently 2N k different terms to be computed, their number is actually lower due to evenness of A  , i.e. to the double of the greatest distance between two grid points, taken along direction of k x .In this case, there are exactly 2 k N different compo- nents of B .For example, the ICs arrays for the first term on the right-hand side of Eq. ( 5) result from , and those for the sec- ) . This treatment, requested by the classical DCFFT algorithm [11], is employed to avoid the so-called periodicity error related to transfer to and from the FD, and its base is detailed in [17].In case of correlation, no additional rearrangement is necessary.
Input the eigenstrains in the computational domain and extend it by zeropadding in all directions Compute the ICs in the SD using existing closed-form expressions [9,14] Select the product type (  The same notation convention can be applied to the digitized eigenstrains, i.e. ,, i j m   denotes any of the eigenstrain tensor components in the cell (i, j, m), 1 mN  , and   is a vector of two-dimensional arrays of eigenstrains, where 1 The treatment for the eigenstrains also depends on the product type.In case of convolution, eigenstrains are simply extended by zero-padding to match the size of the . In case of correlation, after zeropadding, the eigenstrains are rearranged in reversed order with respect to the corresponding direction, i.e. 0 In the following step, the ICs and the eigenstrains arrays are transferred to FD by means of FFT:    . The spectral array of elastic stresses can be computed as element-by-element product in the FD: ˆB . Its SD counterpart is then retrieved by means of an inverse Fourier transform: IFFT( )   , and only the terms in the original computational domain are retained, i.e.   , 1 . By applying the multi-dimensional spectral algorithm described above, the order of computation for solving the inclusion problem in infinite, elastic and isotropic space (states I and II in Fig. 1) is reduced from the existing The efficiency of the newly pro- posed algorithm also stems from the solution of state III, which is achieved directly from Eq. ( 6) as a twodimensional convolution, which can be efficiently computed using a layer-by-layer two-dimensional DCFFT algorithm, as detailed in [16].In this case, although the observation domain is three-dimensional, the source is twodimensional only, comprised in the half-space boundary.The computational complexity for achieving solution of state III is therefore only of order

Program validation
Numerical predictions of the newly advanced algorithm are benchmarked against existing solutions for the inclusion problem involving regularly shaped domains containing eigenstrains.The solution of an ellipsoidal inclusion with uniform dilatational eigenstrains (i.e. , where ij  is the Kronecker delta and e a pre- scribed strain) in an elastic isotropic half-space is presented by Mura [1].Ellipsoid half-axis along direction of , , is denoted by x , depicted in Fig. 5 for the case 3 2 ha  , and in Fig. 6 for the case 3 ha  , match well the results pre- sented in [1], pp.117-119.
Subsequent simulations are performed employing a hemispherical (a 1 = a 2 = a 3 ) inclusion with the base parallel to the free surface.The depth of the centre of the sphere is denoted by h, and the eigenstrains are located in the half-space 3 xh  , as shown in Fig. 4, b.The closed- form solution to this problem is presented in [18], and the imposed stress normaliser is

Extension of results
The newly advanced computer program is used next to predict spherical inclusion interaction, and a critical interaction distance between two neighboring inclusions is advanced for the case when dilatational eigenstrains are uniform or vary linearly with the radius of the sphere.
These particular examples are meant to prove the potential In a first scenario, two spherical inclusions of radius 4 Ra  (with a fixed but otherwise arbitrarily cho- sen), containing uniform dilatational eigenstrains, are considered in a computational domain of side lengths 4aaa , divided into 800 200 200  elementary cuboids.The coordinate system origin matches the center of the top side of the rectangular computational domain, and the 3 x -axis points inwards.Dimensionless coordinates are defined as ratio to a ,  .This proves that, in the investigated configuration, the free surface has a negligibly weak effect on the stress distribution.Consequently, the critical interaction distance found in this case also holds when the sphere of center 2 C is translated along an axis parallel to 3 x .It has been veri- fied that at smaller depths, the 22  and 33  stress tensor components no longer overlap, and the effect of the free surface must be accounted for, although the same critical interaction distance holds.
Curves 1 and 2 in Fig. 9 prove a strong interaction between the tangent spherical inclusions, leading to an increase of almost 100 % in 11  at the contact point.
When the second inclusion is moved so that the distance between 1 C and 2 C is eight times the radius of the sphere, the interaction is small enough to be neglected, as confirmed by the symmetry of left and right branches of curves 3 and 4 with respect to the centers of the spheres, located at 1 1 x  .In this case, stresses are close to uni- form inside the inclusions, in agreement with existing closed-form results [1], and vanish asymptotically at the edges of the computational domain,  ( 1,0,1 2)  .Curves 1 and 2, corresponding to tangent spherical inclusions, show the intensity of the interaction, leading to loss of stress symmetry in both inclusions and to higher stresses at the interface.When the second inclusion is located eight radii away (curves 5 and 6 in Fig. 10), stresses regain symmetry with respect to the center of each sphere.The vanishing stresses on the sides of the computational domain, x  , prove that at this distance, interaction is minimal, and stresses vary almost linearly inside the inclusion.The latter behavior was also observed by Liu and Wang [12] when studying a single spherical inclusion with linearly distributed dilatational eigenstrains.An intermediate configuration is depicted by curves 3 and 4, when maximum stresses at the edge of the inclusions are weakly affected, but stresses in the middle region are distorted due to proximity of the neighboring inclusion.

Conclusions
A refined numerical solution to the inclusion problem, based on the mirror image method, advanced in this paper.Its strong point consists in implementation of spectral methods implementing the convolution and correlation theorems in superposition of effects in three dimensions.An existing technique for rapid computation of convolution products is enhanced by joining convolution and correlation in a hybrid algorithm, and used to efficiently compute stresses due to inclusions in infinite, elastic and isotropic space.
Additional computational efficiency yields from the robust manner of imposing free surface relief, allowing for solution of inclusion problem involving general, not only deviatoric, eigenstrains.A good agreement of the numerical predictions with existing solutions for regularly shaped domains containing uniform dilatational eigenstrains was found.
The computationally efficient new method allows fine resolutions capable of simulating interaction of multiple inclusions.The effect of a second inclusion proximity is investigated through several practical examples.A critical interaction distance, above which the interference between two identical spherical inclusions can be neglected, is advanced.

S. Spinu A REFINED NUMERICAL SOLUTION TO THE INCLUSION PROBLEM
S u m m a r y A refined numerical method for the inclusion problem involving known, but otherwise arbitrary shaped multiple regions with eigenstrains is advanced in this paper.The newly advanced algorithm, employing the mirror image method, derives its effectiveness from implementation of a hybrid multi-dimensional algorithm, based on convolution and correlation theorems.A good validation for various regularly shaped inclusions containing uniform eigenstrains is found.The method is subsequently applied to predict a critical interaction distance between two neighboring spherical inclusions, containing uniform or non-uniform dilatational eigenstrains.


denote the elastic state induced in the half-space x 3  0 by the fictitious normal traction 33 1 2

Fig. 1
Fig. 1 Solution of the inclusion problem employing the mirror image method,

1 2 3 8N
N N different terms for each of the 3 N layers of constant depth.The amount of memory required to store these arrays limits dramatically the achievable resolution.

Fig. 3
Fig. 3 Flow-chart of the algorithm for the hybrid convolution-correlation multi-dimensional product of the ellipsoid center by h , as shown in Fig.4, a. Dimensionless coordinates are defined as ratio to half-axis in the corresponding direction, and dimensionless normal stresses  and  are the elastic constants of the elastic matrix.Normal stresses along direction of 3


e is a prescribed strain.Position of the center of one sphere is kept fixed at 1 stress tensor component is omitted, because its distribution matches closely that of 33

1 2 1 0
x  , as well as in the middle, where x  .