# A refined numerical solution to the inclusion problem/Istobulintas skaitinis intarpu uzdavinio sprendimo metodas.

1. 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 semi-infinite 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 semi-infinite 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 inclusion-related 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.

2. 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 closed-form 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 ([x.sub.1], [x.sub.2], [x.sub.3]) is attached to the center of the cuboid containing uniform eigenstrains [[epsilon].sup.*.sub.ij], i, j = 1,2,3. The six components of the stress tensor [[sigma].sub.ij], i, j = 1,2,3 in a point of coordinates ([x.sub.1], [x.sub.2], [x.sub.3]) (i.e. the observation point), due to the eigenstrains tensor [[epsilon].sup.*.sub.ij] (i.e. the source cuboid), can be expressed as the contribution of all components of the [[epsilon].sup.*.sub.ij] tensor, yielding 36 different ICs:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

Consequently, [A.sub.ijkl], i, j, k, l = 1, 2, 3, denotes the stress tensor component ij induced in the infinite elastic matrix by the k l 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 [not equal to] l, as [A.sub.ijkl] = [A.sub.ijlk], contribution of both corresponding strains, i.e. of both [[epsilon].sup.*.sub.kl] = [[epsilon].sup.*.sub.lk] = 1, is accounted for in Eq. (1) by means of the multiplier 2. [A.sub.ijkl] is therefore a fourt-horder 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):

[[sigma].sub.ij]([x.sub.1], [x.sub.2], [x.sub.3]) = [A.sub.ijkl]([x.sub.1], [x.sub.2], [x.sub.3]) [[epsilon].sup.*.sub.kl] (0,0,0), (2)

where i, j, k, l = 1, 2, 3. Detailed derivation of [A.sub.ijkl] 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. [[epsilon].sup.I.sub.ij] = [[epsilon].sup.*.sub.ij], 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:

[[sigma].sup.I+II.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) = [A.sub.ijkl] ([x.sub.1], [x.sub.2], [x.sub.3]-h) [[epsilon].sup.I.sub.kl] (0,0, h) + ... + + [A.sub.ijkl] ([x.sub.1], [x.sub.2], [x.sub.3] + h) [[epsilon].sup.II.sub.kl] (0,0, -h), (3)

where i, j, k, l = 1, 2, 3. The choice of eigenstrains depicted in Fig. 1 leads to [[sigma].sup.I+II.sub.i3] ([x.sub.1],[x.sub.2],0) = 0 for i = 1, 2. Consequently, the half-space boundary is left free of shear tractions, exhibiting only a fictitious normal traction (pressure) [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],0). According to the method indicated by Mura [15], stresses equal but opposite to [[sigma].sup.I+II.sub.33] ([x.sub.1], [x.sub.2],0) should be applied on the boundary [x.sub.3] = 0, and their contribution subtracted from summation of solutions I and II, i.e. from the [[sigma].sup.I+II.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) computed in Eq. (3), in order to fully satisfy the traction free surface condition in the original inclusion problem.

Let [[sigma].sup.III.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) denote the elastic state induced in the half-space [x.sub.3] [greater than or equal to] 0 by the fictitious normal traction [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],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:

[[sigma].sub.ij]([x.sub.1],[x.sub.2], [x.sub.2]) = [[sigma].sup.I+II.sub.ij] ([x.sub.1],[x.sub.2], [x.sub.3])- [[sigma].sup.III.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]). (4)

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.sub.1], [N.sub.2] and [N.sub.3] respectively, with N = [N.sub.1][N.sub.2][N.sub.3], 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.sub.ijk] of the cuboid (i, j, k) will be denoted by f (i, j, k), and will be computed as f([x.sub.1]([C.sub.ijk]), [x.sub.2]([C.sub.ijk]), [x.sub.3]([C.sub.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.

However, in order to implement spectral methods, which compute multiple elementary contributions simultaneously, the multi-summation process must cover all the cells in the computational domain. Applying this technique to Eq. (3) yields:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

where [xi], [zeta], [??], [gamma] = 1,2,3, 1[less than or equal to]i[less than or equal to][N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1[less than or equal to]k[less than or equal to][N.sub.3]. In this relation, the cuboid (l',m',n'), of centre [C'.sub.lmn], is the mirror image with respect to the plane [x.sub.3] = 0 of the cuboid (l, m, n) centred at [C.sub.lmn]. Consequently, Eq. (4) expresses the [xi][zeta] stress tensor component induced in any observation cuboid (i, j, k) by a digitised distribution of eigenstrains located in the half-space [x.sub.3] [greater than or equal to] 0 (the first term on the right-hand side of Eq. (4)), as well as by its mirror image with respect to the plane [x.sub.3] = 0 (the second term on the right-hand side of Eq. (4)), both laying in an infinite elastic matrix. Due to symmetry, [x.sub.1]([C'.sub.lmn]) = [x.sub.1]([C.sub.lmn]), [x.sub.2]([C'.sub.lmn]) = [x.sub.2]([C.sub.lmn]), but [x.sub.3]([C'.sub.lmn]) = -[x.sub.3]([C.sub.lmn]). 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

It is also convenient to match the set of observation points to the set of centres of all elementary cuboids, i.e. 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] k [less than or equal to] [N.sub.3]. 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 8[N.sub.1][N.sub.2][N.sub.3] different terms for each of the two members on the right-hand side of Eq. (5).

The fictitious normal traction [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],0) acting on the half-space boundary, needed for the solution of state III, yields from Eq. (5) by setting k = 1 and [xi] = [zeta] = 3. The stress induced by this traction in the halfspace [x.sub.3] [greater than or equal to] 0, denoted of [[sigma].sup.I+II.sub.[xi][zeta]](i, j, k), [xi], [zeta] = 1,2,3, 1 [less than or equal to] I [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] k [less than or equal to] [N.sub.3], can subsequently be approximated numerically. Chiu [8] expressed of [[sigma].sup.I+II.sub.[xi][zeta]] in terms of [[epsilon].sup.*.sub.[xi][zeta]] using an existing solution for the stresses in an isotropic half-space under periodically distributed normal surface loading. It should be noted that his method holds only for deviatoric eigenstrains (i.e. tr([[epsilon].sup.*.sub.ii) = 0), 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 8[N.sub.1][N.sub.2][N.sup.2.sub.3] different terms, i.e. 8[N.sub.1][N.sub.2][N.sub.3] different terms for each of the [N.sub.3] layers of constant depth. The amount of memory required to store these arrays limits dramatically the achievable resolution.

The approach newly proposed in this paper differs from that of Chiu in the way [[sigma].sup.III.sub.[xi][zeta]] is related to [[epsilon].sup.*.sub.[xi][zeta]], allowing 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

where [xi],[zeta] = 1,2,3 , 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] m [less than or equal to] [N.sub.3]. The influence coefficient [S.sub.[xi][zeta]] expresses the f stress tensor component induced in the cuboid (i, j, m) by a unit uniform pressure acting on the side of the cuboid (k, l, 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.

3. Computation acceleration

The multi-summation process in Eq. (5) is very computationally intensive, having an order of computation of O([N.sup.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 [[??].sub.3], the global ICs should depend 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 three-dimensional. While the first term in Eq. (5) is a three-dimensional convolution, the second term is a two-dimensional convolution with respect to directions of [[??].sub.1] and of [[??].sub.2], and a one-dimensional correlation with respect to direction of [[??].sub.3]. 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 [[??].sub.3] 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 [[??].sub.3], prior to transfer to the FD. Indeed, 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 s , 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 [??], [xi], [??], [gamma] = 1,2,3. 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. [B.sub.i,j,m], 1 [less than or equal to] i [less than or equal to] 2[N.sub.1], 1 [less than or equal to] j [less than or equal to] 2[N.sub.2], 1 [less than or equal to] m [less than or equal to] 2[N.sub.3], is the (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.sub.i,j,m] as vectors [B.sub.l], 1 [less than or equal to] 1 [less than or equal to] 2[N.sub.k], 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.sub.k] is the number of grids in the direction of [[??].sub.k]. The algorithm advanced herein can also be applied to any multidimensional configurations, i.e. when k > 3. As stated before, computation of [B.sub.l] components depends weather k corresponds to a convolution or a correlation. In case of convolution, [B.sub.1] corresponds to the negative greatest distance between two grid nodes taken along the [x.sub.k]-axis, i.e the observation mark is indexed with 1 and the source mark with [N.sub.k]. The situation when the observation and the source points are interchanged corresponds to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Although there are apparently [2N.sub.k] different terms to be computed, their number is actually lower due to evenness of [A.sub.[xi][zeta][??][gamma]], i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for even [A.sub.[xi][zeta][??][gamma]], and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for odd [A.sub.[xi][zeta][??][gamma]]. If, on the other hand, direction of [[??].sub.k] corresponds to a correlation, [B.sub.1] corresponds to the situation when the source and the observation marks coincide, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to the double of the greatest distance between two grid points, taken along direction of [[??].sub.k]. In this case, there are exactly [2N.sub.k] different components of [B.sub.l]. For example, the ICs arrays for the first term on the right-hand side of Eq. (5) result from [B.sub.i,j,m] = [A.sub.[xi][zeta][??][gamma]] (i -[N.sub.1], j-[N.sub.2],m -[N.sub.3]), and those for the second term from [B.sub.i,j,m] = [A.sub.[xi][zeta][??][gamma]] (i - [N.sub.1], j - [N.sub.2], m - 1), where [??], [xi], [??], [gamma] = 1,2,3 , 1 [less than or equal to] i [less than or equal to] 2[N.sub.1], 1 [less than or equal to] j [less than or equal to] 2[N.sub.2], 1 [less than or equal to] m [less than or equal to] 2[N.sub.3].

The next step consists in reordering of [B.sub.l] terms. In case of convolution, a new vector [[B.bar].sub.l] is generated by zero-padding and rearrangement in wrap-around order, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. 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.

The same notation convention can be applied to the digitized eigenstrains, i.e. [[epsilon].sup.*.sub.i,j,m] denotes any of the eigenstrain tensor components in the cell (i, j, m), 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] m [less than or equal to] [N.sub.3], and [[epsilon].sup.*.sub.l] is a vector of two-dimensional arrays of eigenstrains, where 1 [less than or equal to] l [less than or equal to] [N.sub.k], k = 1,2 or 3. 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 [B.sub.l] vector: [[epsilon.bar]*.sub.l], = [[epsilon.bar]*.sub.l], 1 [less than or equal to] l [less than or equal to] [N.sub.k], [[epsilon.bar]*.sub.l] = 0, [N.sub.k]+1[less than or equal to]l[less than or equal to][2N.sub.k]. In case of correlation, after zeropadding, the eigenstrains are rearranged in reversed order with respect to the corresponding direction, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the following step, the ICs and the eigenstrains arrays are transferred to FD by means of FFT: [[??].bar] FFT([[B].bar]), [[??].bar]* = FFT([[epsilon].bar]*). The spectral array of elastic stresses can be computed as element-by-element product in the FD: [[sigma].bar].sub.l] =[[??].bar].sub.l] [[??]*.bar].sub.l], 1 [less than or equal to] l [less than or equal to][2N.sub.k]. Its SD counterpart is then retrieved by means of an inverse Fourier transform: [[sigma].bar] = IFFT([[??].bar]), and only the terms in the original computational domain are retained, i.e. [[sigma].sub.l] , = [[sigma].sub.l], 1 [less than or equal to] l [less than or equal to] [N.sub.k].

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 O([N.sup.2.sub.3][N.sup.1][N.sup.2] log [N.sup.1] [N.sup.2]) in [9, 10], to a further improved O([N.sup.1][N.sup.2][N.sup.3] log [N.sup.1][N.sup.2][N.sup.3]). The efficiency of the newly proposed algorithm also stems from the solution of state III, which is achieved directly from Eq. (6) as a two-dimensional 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 two-dimensional only, comprised in the half-space boundary. The computational complexity for achieving solution of state III is therefore only of order O([N.sup.3][N.sup.1][N.sup.2] log [N.sup.1][N.sup.2]).

4. 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. [[epsilon].sup.*.sub.ij] = [[delta].sub.ij]e, where [[delta].sub.ij] is the Kronecker delta and e a prescribed strain) in an elastic isotropic half-space is presented by Mura [1]. Ellipsoid half-axis along direction of [[??].sub.i], i = 1,2,3, is denoted by [a.sub.i] ([a.sub.1] = [a.sub.2] in all simulations), and the depth 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 [[??].sub.ii] as ratio to the normalizer [[sigma].sub.0] = 2[mu]e(1 + v)/(1 - v), where [mu] and v are the elastic constants of the elastic matrix. Normal stresses along direction of [x.sub.3], depicted in Fig. 5 for the case h = 2[a.sub.3], and in Fig. 6 for the case h = [a.sub.3], match well the results presented in [1], pp. 117-119.

Subsequent simulations are performed employing a hemispherical ([a.sub.1] = [a.sub.2] = [a.sub.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 [x.sub.3] [greater than or equal to] h, as shown in Fig. 4, b. The closed-form solution to this problem is presented in [18], and the imposed stress normaliser is [[sigma].sub.0] = 2[mu]e. In Fig. 7, uniform dilatational eigenstrains, i.e. [[epsilon].sup.*.sub.ij] = [[delta].sub.ij]e, are imposed, and in Fig. 8, [[epsilon].sup.*.sub.11] = [[epsilon].sup.*.sub.33] = e, [[epsilon].sup.*.sub.22] = 0.5e, [[epsilon].sup.*.sub.ij] = 0, i [not equal to] j.

5. 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 of the method to simulate a large variety of scenarios involving inclusions.

In a first scenario, two spherical inclusions of radius R = a/4 (with a fixed but otherwise arbitrarily chosen), containing uniform dilatational eigenstrains, are considered in a computational domain of side lengths 4a x a x a, divided into 800 x 200 x 200 elementary cuboids. The coordinate system origin matches the center of the top side of the rectangular computational domain, and the [x.sub.3]-axis points inwards. Dimensionless coordinates are defined as ratio to a, [[bar.x].sub.i] = [x.sub.i]/a, i = 1,2 ,3 , and stresses are normalized by the quantity [[sigma].sub.0] = 2[mu]e(1 + v)/(1-v), where e is a prescribed strain. Position of the center of one sphere is kept fixed at [C.sub.1](-1,0,1/2), while the other one is moved along an axis [x'.sub.1], parallel to [x.sub.1], passing through [C.sub.1]. Stresses along the [x'.sub.1]-axis are presented in Fig. 9, for various positions of [C.sub.2]. The [[sigma].sub.22] stress tensor component is omitted, because its distribution matches closely that of [[sigma].sub.33]. 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 [C.sub.2] is translated along an axis parallel to [x.sub.3]. It has been verified that at smaller depths, the [[sigma].sub.22] and [[sigma].sub.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 [[bar.[sigma]].sub.11] at the contact point. When the second inclusion is moved so that the distance between [C.sub.1] and [C.sub.2] 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 [[bar.x]'.sub.1] = [+ or -] 1. In this case, stresses are close to uniform inside the inclusions, in agreement with existing closed-form results [1], and vanish asymptotically at the edges of the computational domain, [[bar.x]'.sub.1] = [+ or -]2, as well as in the middle, where [[bar.x]'.sub.1] = 0.

In a second scenario, the eigenstrains vary linearly with the radius of the sphere, i.e. [[epsilon].sup.*.sub.ij] (r) = [[delta].sub.ij]e (R - r)/R in a spherical coordinate system linked to the center of each sphere. Again, [[delta].sub.ij] is the Kronecker delta and e a prescribed strain. The results for this configuration are presented in Fig. 10, considering three positions of [C.sub.2], while [C.sub.1] is kept fixed at coordinates (-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, [[bar.x]'.sub.1] = [+ or -] 2, and in-between the spherical inclusions, [[bar.x]'.sub.1] = 0, 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.

6. Conclusions

A refined numerical solution to the inclusion problem, based on the mirror image method, is 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.

10.5755/j01.mech.19.3.4657

Received December 06, 2011 Accepted May 15, 2013

References

[1.] Mura, T. 1987. Micromechanics of Defects in Solids, Martinus Nijhoff: Kluwer Academic Publishers, 587 p.

[2.] Eshelby, J.D. 1959. The elastic field outside an ellipsoidal inclusion, Proceedings of the Royal Society of London A252:561-569.

[3.] Mura, T. 1988. Inclusion problem, ASME Applied Mechanics Review 41(1): 15-20. http://dx.doi.org/10.1115/1.3151875.

[4.] Mindlin, R.D.; Cheng, D.H. 1950. Thermoelastic stress in the semi-infinite solid, Journal of Applied Physics 21(9): 931-933. http://dx.doi.org/10.1063/1.1699786.

[5.] Seo, K.; Mura, T. 1979. The elastic field in half space due to ellipsoidal inclusions with uniform dilatational eigenstrains, ASME Journal of Applied Mechanics 46(3): 568-572. http://dx.doi.org/10.1115/1.3424607.

[6.] Mindlin, R.D. 1936. Force at a point in the interior of a semi-infinite solid, Physics 7(5): 195-202. http://dx.doi.org/10.1063/1.1745385.

[7.] Chiu, Y.P. 1977. On the stress field due to initial strains in a cuboid surrounded by an infinite elastic space, ASME Journal of Applied Mechanics 44(4): 587-590. http://dx.doi.org/10.1115/1.3424140.

[8.] Chiu, Y.P. 1978. On the stress field and surface deformation in a half space with cuboidal zone in which initial strains are uniform, ASME Journal of Applied Mechanics 45(2): 302-306. http://dx.doi.org/10.1115/1.3424292.

[9.] Jacq, C.; Nelias, D.; Lormand, G.; Girodin, D. 2002. Development of a three-dimensional semi-analytical elastic-plastic contact code, ASME Journal of Tribology 124(4): 653-667. http://dx.doi.org/10.1115/1.1467920.

[10.] Wang, F.; Keer, L.M. 2005. Numerical simulation for three dimensional elastic-plastic contact with hardening behavior, ASME Journal of Tribology 127(3): 494-502. http://dx.doi.org/10.1115/1.1924573.

[11.] Liu, S.B.; Wang, Q.; Liu, G. 2000. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses, Wear 243(1-2): 101-111. http://dx.doi.org/10.1016/S0043-1648(00)00427-0.

[12.] Liu, S.; Wang, Q. 2005. Elastic fields due to eigenstrains in a half-space, ASME Journal of Applied Mechanics 72(6): 871-878. http://dx.doi.org/10.1115/1.2047598.

[13.] Mindlin, R.D.; Cheng, D.H. 1950. Nuclei of strain in the semi-infinite solid, Journal of Applied Physics 21(9): 926-930. http://dx.doi.org/10.1063/L1699785.

[14.] Zhou, K.; Chen, W.W.; Keer, L.M.; Wang, Q.J. 2009. A fast method for solving three-dimensional arbitrarily shaped inclusions in a half-space, Computer Methods in Applied Mechanics and Engineering 198(9 12): 885-892. http://dx.doi.org/10.1016/j.cma.2008.10.021.

[15.] Mura, T. 1968. The continuum theory of dislocation, in Advances in Material Research, Vol. 3.-Herbert Herman: Interscience Publisher: 1-108.

[16.] Liu, S.B.; Wang, Q. 2002. Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm, ASME Journal of Tribology 124(1): 36-45. http://dx.doi.org/10.1115/1.1401017.

[17.] Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. 1992. Numerical Recipes in C: The Art of Scientific Computing, Second Edition, Cambridge University Press. 994p.

[18.] Wu, L. 2003. The elastic field induced by a hemispherical inclusion in the half-space, Acta Mechanica Sinica 19(3): 253-262. http://dx.doi.org/10.1007/BF02484488.

S. Spinu

University of Suceava, 13th University Street, 720229, Suceava, Romania, E-mail: sergiu.spinu@fim.usv.ro

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 semi-infinite 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 semi-infinite 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 inclusion-related 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.

2. 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 closed-form 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 ([x.sub.1], [x.sub.2], [x.sub.3]) is attached to the center of the cuboid containing uniform eigenstrains [[epsilon].sup.*.sub.ij], i, j = 1,2,3. The six components of the stress tensor [[sigma].sub.ij], i, j = 1,2,3 in a point of coordinates ([x.sub.1], [x.sub.2], [x.sub.3]) (i.e. the observation point), due to the eigenstrains tensor [[epsilon].sup.*.sub.ij] (i.e. the source cuboid), can be expressed as the contribution of all components of the [[epsilon].sup.*.sub.ij] tensor, yielding 36 different ICs:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

Consequently, [A.sub.ijkl], i, j, k, l = 1, 2, 3, denotes the stress tensor component ij induced in the infinite elastic matrix by the k l 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 [not equal to] l, as [A.sub.ijkl] = [A.sub.ijlk], contribution of both corresponding strains, i.e. of both [[epsilon].sup.*.sub.kl] = [[epsilon].sup.*.sub.lk] = 1, is accounted for in Eq. (1) by means of the multiplier 2. [A.sub.ijkl] is therefore a fourt-horder 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):

[[sigma].sub.ij]([x.sub.1], [x.sub.2], [x.sub.3]) = [A.sub.ijkl]([x.sub.1], [x.sub.2], [x.sub.3]) [[epsilon].sup.*.sub.kl] (0,0,0), (2)

where i, j, k, l = 1, 2, 3. Detailed derivation of [A.sub.ijkl] 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. [[epsilon].sup.I.sub.ij] = [[epsilon].sup.*.sub.ij], 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:

[[sigma].sup.I+II.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) = [A.sub.ijkl] ([x.sub.1], [x.sub.2], [x.sub.3]-h) [[epsilon].sup.I.sub.kl] (0,0, h) + ... + + [A.sub.ijkl] ([x.sub.1], [x.sub.2], [x.sub.3] + h) [[epsilon].sup.II.sub.kl] (0,0, -h), (3)

where i, j, k, l = 1, 2, 3. The choice of eigenstrains depicted in Fig. 1 leads to [[sigma].sup.I+II.sub.i3] ([x.sub.1],[x.sub.2],0) = 0 for i = 1, 2. Consequently, the half-space boundary is left free of shear tractions, exhibiting only a fictitious normal traction (pressure) [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],0). According to the method indicated by Mura [15], stresses equal but opposite to [[sigma].sup.I+II.sub.33] ([x.sub.1], [x.sub.2],0) should be applied on the boundary [x.sub.3] = 0, and their contribution subtracted from summation of solutions I and II, i.e. from the [[sigma].sup.I+II.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) computed in Eq. (3), in order to fully satisfy the traction free surface condition in the original inclusion problem.

Let [[sigma].sup.III.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]) denote the elastic state induced in the half-space [x.sub.3] [greater than or equal to] 0 by the fictitious normal traction [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],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:

[[sigma].sub.ij]([x.sub.1],[x.sub.2], [x.sub.2]) = [[sigma].sup.I+II.sub.ij] ([x.sub.1],[x.sub.2], [x.sub.3])- [[sigma].sup.III.sub.ij] ([x.sub.1], [x.sub.2], [x.sub.3]). (4)

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.sub.1], [N.sub.2] and [N.sub.3] respectively, with N = [N.sub.1][N.sub.2][N.sub.3], 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.sub.ijk] of the cuboid (i, j, k) will be denoted by f (i, j, k), and will be computed as f([x.sub.1]([C.sub.ijk]), [x.sub.2]([C.sub.ijk]), [x.sub.3]([C.sub.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.

However, in order to implement spectral methods, which compute multiple elementary contributions simultaneously, the multi-summation process must cover all the cells in the computational domain. Applying this technique to Eq. (3) yields:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

where [xi], [zeta], [??], [gamma] = 1,2,3, 1[less than or equal to]i[less than or equal to][N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1[less than or equal to]k[less than or equal to][N.sub.3]. In this relation, the cuboid (l',m',n'), of centre [C'.sub.lmn], is the mirror image with respect to the plane [x.sub.3] = 0 of the cuboid (l, m, n) centred at [C.sub.lmn]. Consequently, Eq. (4) expresses the [xi][zeta] stress tensor component induced in any observation cuboid (i, j, k) by a digitised distribution of eigenstrains located in the half-space [x.sub.3] [greater than or equal to] 0 (the first term on the right-hand side of Eq. (4)), as well as by its mirror image with respect to the plane [x.sub.3] = 0 (the second term on the right-hand side of Eq. (4)), both laying in an infinite elastic matrix. Due to symmetry, [x.sub.1]([C'.sub.lmn]) = [x.sub.1]([C.sub.lmn]), [x.sub.2]([C'.sub.lmn]) = [x.sub.2]([C.sub.lmn]), but [x.sub.3]([C'.sub.lmn]) = -[x.sub.3]([C.sub.lmn]). 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

It is also convenient to match the set of observation points to the set of centres of all elementary cuboids, i.e. 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] k [less than or equal to] [N.sub.3]. 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 8[N.sub.1][N.sub.2][N.sub.3] different terms for each of the two members on the right-hand side of Eq. (5).

The fictitious normal traction [[sigma].sup.I+II.sub.33] ([x.sub.1],[x.sub.2],0) acting on the half-space boundary, needed for the solution of state III, yields from Eq. (5) by setting k = 1 and [xi] = [zeta] = 3. The stress induced by this traction in the halfspace [x.sub.3] [greater than or equal to] 0, denoted of [[sigma].sup.I+II.sub.[xi][zeta]](i, j, k), [xi], [zeta] = 1,2,3, 1 [less than or equal to] I [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] k [less than or equal to] [N.sub.3], can subsequently be approximated numerically. Chiu [8] expressed of [[sigma].sup.I+II.sub.[xi][zeta]] in terms of [[epsilon].sup.*.sub.[xi][zeta]] using an existing solution for the stresses in an isotropic half-space under periodically distributed normal surface loading. It should be noted that his method holds only for deviatoric eigenstrains (i.e. tr([[epsilon].sup.*.sub.ii) = 0), 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 8[N.sub.1][N.sub.2][N.sup.2.sub.3] different terms, i.e. 8[N.sub.1][N.sub.2][N.sub.3] different terms for each of the [N.sub.3] layers of constant depth. The amount of memory required to store these arrays limits dramatically the achievable resolution.

The approach newly proposed in this paper differs from that of Chiu in the way [[sigma].sup.III.sub.[xi][zeta]] is related to [[epsilon].sup.*.sub.[xi][zeta]], allowing 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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

where [xi],[zeta] = 1,2,3 , 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] m [less than or equal to] [N.sub.3]. The influence coefficient [S.sub.[xi][zeta]] expresses the f stress tensor component induced in the cuboid (i, j, m) by a unit uniform pressure acting on the side of the cuboid (k, l, 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.

3. Computation acceleration

The multi-summation process in Eq. (5) is very computationally intensive, having an order of computation of O([N.sup.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 [[??].sub.3], the global ICs should depend 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 three-dimensional. While the first term in Eq. (5) is a three-dimensional convolution, the second term is a two-dimensional convolution with respect to directions of [[??].sub.1] and of [[??].sub.2], and a one-dimensional correlation with respect to direction of [[??].sub.3]. 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 [[??].sub.3] 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 [[??].sub.3], prior to transfer to the FD. Indeed, 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 s , 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 [??], [xi], [??], [gamma] = 1,2,3. 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. [B.sub.i,j,m], 1 [less than or equal to] i [less than or equal to] 2[N.sub.1], 1 [less than or equal to] j [less than or equal to] 2[N.sub.2], 1 [less than or equal to] m [less than or equal to] 2[N.sub.3], is the (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.sub.i,j,m] as vectors [B.sub.l], 1 [less than or equal to] 1 [less than or equal to] 2[N.sub.k], 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.sub.k] is the number of grids in the direction of [[??].sub.k]. The algorithm advanced herein can also be applied to any multidimensional configurations, i.e. when k > 3. As stated before, computation of [B.sub.l] components depends weather k corresponds to a convolution or a correlation. In case of convolution, [B.sub.1] corresponds to the negative greatest distance between two grid nodes taken along the [x.sub.k]-axis, i.e the observation mark is indexed with 1 and the source mark with [N.sub.k]. The situation when the observation and the source points are interchanged corresponds to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Although there are apparently [2N.sub.k] different terms to be computed, their number is actually lower due to evenness of [A.sub.[xi][zeta][??][gamma]], i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for even [A.sub.[xi][zeta][??][gamma]], and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for odd [A.sub.[xi][zeta][??][gamma]]. If, on the other hand, direction of [[??].sub.k] corresponds to a correlation, [B.sub.1] corresponds to the situation when the source and the observation marks coincide, and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] to the double of the greatest distance between two grid points, taken along direction of [[??].sub.k]. In this case, there are exactly [2N.sub.k] different components of [B.sub.l]. For example, the ICs arrays for the first term on the right-hand side of Eq. (5) result from [B.sub.i,j,m] = [A.sub.[xi][zeta][??][gamma]] (i -[N.sub.1], j-[N.sub.2],m -[N.sub.3]), and those for the second term from [B.sub.i,j,m] = [A.sub.[xi][zeta][??][gamma]] (i - [N.sub.1], j - [N.sub.2], m - 1), where [??], [xi], [??], [gamma] = 1,2,3 , 1 [less than or equal to] i [less than or equal to] 2[N.sub.1], 1 [less than or equal to] j [less than or equal to] 2[N.sub.2], 1 [less than or equal to] m [less than or equal to] 2[N.sub.3].

The next step consists in reordering of [B.sub.l] terms. In case of convolution, a new vector [[B.bar].sub.l] is generated by zero-padding and rearrangement in wrap-around order, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. 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.

The same notation convention can be applied to the digitized eigenstrains, i.e. [[epsilon].sup.*.sub.i,j,m] denotes any of the eigenstrain tensor components in the cell (i, j, m), 1 [less than or equal to] i [less than or equal to] [N.sub.1], 1 [less than or equal to] j [less than or equal to] [N.sub.2], 1 [less than or equal to] m [less than or equal to] [N.sub.3], and [[epsilon].sup.*.sub.l] is a vector of two-dimensional arrays of eigenstrains, where 1 [less than or equal to] l [less than or equal to] [N.sub.k], k = 1,2 or 3. 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 [B.sub.l] vector: [[epsilon.bar]*.sub.l], = [[epsilon.bar]*.sub.l], 1 [less than or equal to] l [less than or equal to] [N.sub.k], [[epsilon.bar]*.sub.l] = 0, [N.sub.k]+1[less than or equal to]l[less than or equal to][2N.sub.k]. In case of correlation, after zeropadding, the eigenstrains are rearranged in reversed order with respect to the corresponding direction, i.e. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

In the following step, the ICs and the eigenstrains arrays are transferred to FD by means of FFT: [[??].bar] FFT([[B].bar]), [[??].bar]* = FFT([[epsilon].bar]*). The spectral array of elastic stresses can be computed as element-by-element product in the FD: [[sigma].bar].sub.l] =[[??].bar].sub.l] [[??]*.bar].sub.l], 1 [less than or equal to] l [less than or equal to][2N.sub.k]. Its SD counterpart is then retrieved by means of an inverse Fourier transform: [[sigma].bar] = IFFT([[??].bar]), and only the terms in the original computational domain are retained, i.e. [[sigma].sub.l] , = [[sigma].sub.l], 1 [less than or equal to] l [less than or equal to] [N.sub.k].

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 O([N.sup.2.sub.3][N.sup.1][N.sup.2] log [N.sup.1] [N.sup.2]) in [9, 10], to a further improved O([N.sup.1][N.sup.2][N.sup.3] log [N.sup.1][N.sup.2][N.sup.3]). The efficiency of the newly proposed algorithm also stems from the solution of state III, which is achieved directly from Eq. (6) as a two-dimensional 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 two-dimensional only, comprised in the half-space boundary. The computational complexity for achieving solution of state III is therefore only of order O([N.sup.3][N.sup.1][N.sup.2] log [N.sup.1][N.sup.2]).

4. 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. [[epsilon].sup.*.sub.ij] = [[delta].sub.ij]e, where [[delta].sub.ij] is the Kronecker delta and e a prescribed strain) in an elastic isotropic half-space is presented by Mura [1]. Ellipsoid half-axis along direction of [[??].sub.i], i = 1,2,3, is denoted by [a.sub.i] ([a.sub.1] = [a.sub.2] in all simulations), and the depth 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 [[??].sub.ii] as ratio to the normalizer [[sigma].sub.0] = 2[mu]e(1 + v)/(1 - v), where [mu] and v are the elastic constants of the elastic matrix. Normal stresses along direction of [x.sub.3], depicted in Fig. 5 for the case h = 2[a.sub.3], and in Fig. 6 for the case h = [a.sub.3], match well the results presented in [1], pp. 117-119.

Subsequent simulations are performed employing a hemispherical ([a.sub.1] = [a.sub.2] = [a.sub.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 [x.sub.3] [greater than or equal to] h, as shown in Fig. 4, b. The closed-form solution to this problem is presented in [18], and the imposed stress normaliser is [[sigma].sub.0] = 2[mu]e. In Fig. 7, uniform dilatational eigenstrains, i.e. [[epsilon].sup.*.sub.ij] = [[delta].sub.ij]e, are imposed, and in Fig. 8, [[epsilon].sup.*.sub.11] = [[epsilon].sup.*.sub.33] = e, [[epsilon].sup.*.sub.22] = 0.5e, [[epsilon].sup.*.sub.ij] = 0, i [not equal to] j.

5. 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 of the method to simulate a large variety of scenarios involving inclusions.

In a first scenario, two spherical inclusions of radius R = a/4 (with a fixed but otherwise arbitrarily chosen), containing uniform dilatational eigenstrains, are considered in a computational domain of side lengths 4a x a x a, divided into 800 x 200 x 200 elementary cuboids. The coordinate system origin matches the center of the top side of the rectangular computational domain, and the [x.sub.3]-axis points inwards. Dimensionless coordinates are defined as ratio to a, [[bar.x].sub.i] = [x.sub.i]/a, i = 1,2 ,3 , and stresses are normalized by the quantity [[sigma].sub.0] = 2[mu]e(1 + v)/(1-v), where e is a prescribed strain. Position of the center of one sphere is kept fixed at [C.sub.1](-1,0,1/2), while the other one is moved along an axis [x'.sub.1], parallel to [x.sub.1], passing through [C.sub.1]. Stresses along the [x'.sub.1]-axis are presented in Fig. 9, for various positions of [C.sub.2]. The [[sigma].sub.22] stress tensor component is omitted, because its distribution matches closely that of [[sigma].sub.33]. 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 [C.sub.2] is translated along an axis parallel to [x.sub.3]. It has been verified that at smaller depths, the [[sigma].sub.22] and [[sigma].sub.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 [[bar.[sigma]].sub.11] at the contact point. When the second inclusion is moved so that the distance between [C.sub.1] and [C.sub.2] 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 [[bar.x]'.sub.1] = [+ or -] 1. In this case, stresses are close to uniform inside the inclusions, in agreement with existing closed-form results [1], and vanish asymptotically at the edges of the computational domain, [[bar.x]'.sub.1] = [+ or -]2, as well as in the middle, where [[bar.x]'.sub.1] = 0.

In a second scenario, the eigenstrains vary linearly with the radius of the sphere, i.e. [[epsilon].sup.*.sub.ij] (r) = [[delta].sub.ij]e (R - r)/R in a spherical coordinate system linked to the center of each sphere. Again, [[delta].sub.ij] is the Kronecker delta and e a prescribed strain. The results for this configuration are presented in Fig. 10, considering three positions of [C.sub.2], while [C.sub.1] is kept fixed at coordinates (-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, [[bar.x]'.sub.1] = [+ or -] 2, and in-between the spherical inclusions, [[bar.x]'.sub.1] = 0, 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.

6. Conclusions

A refined numerical solution to the inclusion problem, based on the mirror image method, is 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.

10.5755/j01.mech.19.3.4657

Received December 06, 2011 Accepted May 15, 2013

References

[1.] Mura, T. 1987. Micromechanics of Defects in Solids, Martinus Nijhoff: Kluwer Academic Publishers, 587 p.

[2.] Eshelby, J.D. 1959. The elastic field outside an ellipsoidal inclusion, Proceedings of the Royal Society of London A252:561-569.

[3.] Mura, T. 1988. Inclusion problem, ASME Applied Mechanics Review 41(1): 15-20. http://dx.doi.org/10.1115/1.3151875.

[4.] Mindlin, R.D.; Cheng, D.H. 1950. Thermoelastic stress in the semi-infinite solid, Journal of Applied Physics 21(9): 931-933. http://dx.doi.org/10.1063/1.1699786.

[5.] Seo, K.; Mura, T. 1979. The elastic field in half space due to ellipsoidal inclusions with uniform dilatational eigenstrains, ASME Journal of Applied Mechanics 46(3): 568-572. http://dx.doi.org/10.1115/1.3424607.

[6.] Mindlin, R.D. 1936. Force at a point in the interior of a semi-infinite solid, Physics 7(5): 195-202. http://dx.doi.org/10.1063/1.1745385.

[7.] Chiu, Y.P. 1977. On the stress field due to initial strains in a cuboid surrounded by an infinite elastic space, ASME Journal of Applied Mechanics 44(4): 587-590. http://dx.doi.org/10.1115/1.3424140.

[8.] Chiu, Y.P. 1978. On the stress field and surface deformation in a half space with cuboidal zone in which initial strains are uniform, ASME Journal of Applied Mechanics 45(2): 302-306. http://dx.doi.org/10.1115/1.3424292.

[9.] Jacq, C.; Nelias, D.; Lormand, G.; Girodin, D. 2002. Development of a three-dimensional semi-analytical elastic-plastic contact code, ASME Journal of Tribology 124(4): 653-667. http://dx.doi.org/10.1115/1.1467920.

[10.] Wang, F.; Keer, L.M. 2005. Numerical simulation for three dimensional elastic-plastic contact with hardening behavior, ASME Journal of Tribology 127(3): 494-502. http://dx.doi.org/10.1115/1.1924573.

[11.] Liu, S.B.; Wang, Q.; Liu, G. 2000. A versatile method of discrete convolution and FFT (DC-FFT) for contact analyses, Wear 243(1-2): 101-111. http://dx.doi.org/10.1016/S0043-1648(00)00427-0.

[12.] Liu, S.; Wang, Q. 2005. Elastic fields due to eigenstrains in a half-space, ASME Journal of Applied Mechanics 72(6): 871-878. http://dx.doi.org/10.1115/1.2047598.

[13.] Mindlin, R.D.; Cheng, D.H. 1950. Nuclei of strain in the semi-infinite solid, Journal of Applied Physics 21(9): 926-930. http://dx.doi.org/10.1063/L1699785.

[14.] Zhou, K.; Chen, W.W.; Keer, L.M.; Wang, Q.J. 2009. A fast method for solving three-dimensional arbitrarily shaped inclusions in a half-space, Computer Methods in Applied Mechanics and Engineering 198(9 12): 885-892. http://dx.doi.org/10.1016/j.cma.2008.10.021.

[15.] Mura, T. 1968. The continuum theory of dislocation, in Advances in Material Research, Vol. 3.-Herbert Herman: Interscience Publisher: 1-108.

[16.] Liu, S.B.; Wang, Q. 2002. Studying contact stress fields caused by surface tractions with a discrete convolution and fast fourier transform algorithm, ASME Journal of Tribology 124(1): 36-45. http://dx.doi.org/10.1115/1.1401017.

[17.] Press, W.H.; Teukolsky, S.A.; Vetterling, W.T.; Flannery, B.P. 1992. Numerical Recipes in C: The Art of Scientific Computing, Second Edition, Cambridge University Press. 994p.

[18.] Wu, L. 2003. The elastic field induced by a hemispherical inclusion in the half-space, Acta Mechanica Sinica 19(3): 253-262. http://dx.doi.org/10.1007/BF02484488.

S. Spinu

University of Suceava, 13th University Street, 720229, Suceava, Romania, E-mail: sergiu.spinu@fim.usv.ro

Printer friendly Cite/link Email Feedback | |

Author: | Spinu, S. |
---|---|

Publication: | Mechanika |

Article Type: | Report |

Geographic Code: | 4EXRO |

Date: | May 1, 2013 |

Words: | 6299 |

Previous Article: | Investigation of hardening and tempering deformations of carburized low alloy steel/Ianglinto legiruoto plieno grudinimo ir atleidimo deformaciju... |

Next Article: | Influence of cantilever length on stress distribution in fixation screws of All-on-4 full-arch bridge/ Gembes ilgio itaka viso dantu lanko protezo... |

Topics: |