## Abstract

The plastic zone developed during elastoplastic contact may be effectively modeled as an inclusion in an isotropic half space. This paper proposes a simple but efficient computational method to analyze the stresses caused by near surface inclusions of arbitrary shape. The solution starts by solving a corresponding full space inclusion problem and proceeds to annul the stresses acting normal and tangential to the surface, where the numerical computations are processed by taking advantage of the fast Fourier transform techniques with a parallel computing strategy. The extreme case of a cuboidal inclusion with one facet on the surface of the half space is chosen to validate the method. When the surface truncation domain is extended sufficiently and the grids are dense enough, the results based on the new approach are in good agreement with the exact solutions. When solving a typical elastoplastic contact problem, the present analysis is roughly two times faster than the image inclusion approach and six times faster than the direct method. In addition, the present work demonstrates that a significant enhancement in the computational efficiency can be achieved through the introduction of parallel computation.

## Introduction

Inhomogeneous and inelastic strains occur during a loading process and affect the elastic fields of contacting elements. Mura [1] notes that the inelastic strains can be treated as eigenstrains and the regions with eigenstrains termed as inclusions whose elastic moduli are the same as the matrix region [2]. Moreover, the stress disturbance due to inhomogeneities can be equivalent to the eigenstress resulting from corresponding inclusions with proper eigenstrains, known as the equivalent inclusion method (EIM) [2]. Understanding the elastic fields caused by inclusions provides an effective and convenient approach to solving complex contact problems that require a large amount of computation using conventional methods.

The elastic fields caused by the inclusions in a half space were investigated by a number of researchers involving several types of inclusions, such as spherical thermal inclusions [3], ellipsoidal thermal inclusions [4], cuboidal inclusions [5], and periodic distribution of dislocations [6,7]. For cases involving inclusions of more complicated geometry or structure, an analytical solution is generally difficult to obtain, and therefore a numerical solution is sought. A cuboidal inclusion solution is fundamental to conducting numerical computation since an arbitrarily shaped inclusion can be divided into several small cuboids, and the elastic fields calculated by superposing the contributions from each individual cuboidal inclusion. Chiu [5] employed the method of images to obtain closed form solutions to stresses and surface normal displacements. In his method, the relation between the eigenstress and eigenstrain is expressed as a recursive formula, without explicit demonstration of the convolution/correlation properties in the depth direction. Based on Chiu's method, Jacq et al. [8] used a two-dimensional fast Fourier transform (2D-FFT) to calculate the residual stresses caused by the plastic strains, where the plastic zone was uniformly discretized into cuboidal inclusions. The stresses in their method were calculated layer by layer through the depth direction. This method has been widely used to analyze several elasoplastic contact problems [9–14]. Liu and Wang [15] demonstrated the explicit correlation properties in the depth direction with respect to the image inclusion; therefore, the discrete correlation using fast Fourier transform (DCR-FFT) can be employed to solve for the stresses related to the image inclusion, while the stress arising from its own part, i.e., the inclusion in a full space, can be obtained by discrete convolution, using fast Fourier transform (DC-FFT). According to the convolution/correlation properties in the depth direction, Zhou et al. [16] used two three-dimensional fast Fourier transform (3D-FFT) procedures to solve the stresses due to the cuboidal inclusions and its image part in a full space, and one 2D-FFT procedure to solve the stress caused by surface normal traction. By using 3D-FFT, the method proposed by Zhou et al. [16] is significantly faster than the layer-by-layer 2D-FFT of Jacq et al. [8]; however the computational domain required appropriate extension to reduce the numerical error caused by normal traction cancellation.

Recently, Liu et al. [17] employing Galerkin vectors derived analytical solutions to the elastic fields caused by eigenstrains in a half space. Explicit closed-form results were obtained for the influence coefficients relating eigenstrains to displacements or stresses. By using the influence coefficients, a numerical implementation based on 3D-FFT could be employed to speed up the computation, and the elastic field caused by arbitrarily shaped inclusions calculated. This method has been successfully used to analyze the contact involving joined quarter spaces [18]. The analytical solutions of influence coefficients can bypass the additional complexity of surface domain truncation but requires more computational time than the Zhou et al. [16] since four 3D-FFT schemes must be used.

Considering the rapid development of multicore processors in recent years, parallel computing has been increasingly used in solving eigenstrains and related problems, such as dislocation dynamics [19], crack propagation [20], and composite materials [21]. To employ parallel computing, a large problem is divided into smaller independent components, where each processor executes its own part of the algorithm simultaneously, leading to a significant increase in computational speed; and the theoretical maximum acceleration can be predicted by Amdahl's law [22]. In the present problem, calculating influence coefficients and performing the FFT consume most of the numerical implementation time. Since influence coefficients are matrix components independent of each other, parallel computing can be employed to accelerate the computation speed. Moreover, the parallel version of FFT is now well developed and is adopted in this work.

The purpose of this paper is to develop a simpler and more rapid method to calculate the stresses caused by the presence of arbitrarily shaped inclusions in contact problems. Since the new method is only based on one full space solution, it is significantly faster and more efficient than previous methods. In addition, parallel computing is employed in the new solver to accelerate computation speed.

## Theory of Inclusions in Materials Under Contact

### Theoretical Basis of the New Method.

Contact problems can be analyzed in a half space, as shown in Fig. 1. Here the stresses are the sums of those of a full space solution and those caused by tractions required to create the free surface. The full space solutions can be obtained by using Chiu's method [23] or Liu et al. [17].

Fig. 1
Fig. 1
Close modal

#### Chiu's Full Space Solution.

The gradient of displacement components at a target point $x$ in a full space caused by a cuboidal inclusion of uniform eigenstrains $εij*$ can be expressed as [23]
$2μui,q(x)= (8π3)-1∑n=18(-1)n[1-2ν1-νλεkk*·D,iqmm(cn) + 4μεij*·D,jqmm(cn)-2μ1-νεmj*·D,iqmj(cn)]$
(1)
where λ and μ are the Lamé constants, ν is Poisson's ratio, and cn are the vectors linking the eight corners of the cuboidal inclusion to the target point, the comma indicates partial differentiation, and D is as follows [23]:
$D(cn)=∫ξ3=-∞∞∫ξ2=-∞∞∫ξ1=-∞∞iexp(iξ·cn)(ξ12+ξ22+ξ32)2ξ1ξ2ξ3dξ1dξ2dξ3$
(2)
Once the displacement gradients are calculated from Eq. (1), the eigenstress can be obtained from Hooke's law,
$σij*(x)=Cijkl[uk.l(x)+ul,k(x)]/2, x∉Ω$
(3a)
$σij*(x)=Cijkl{[uk.l(x)+ul,k(x)]/2-εkl*}, x∈Ω$
(3b)

where the cuboidal inclusion is in domain Ω.

#### Full Space Solution of Liu et al.

The full space solution can be expressed as follows [17]:
$σij*(x)=-μ4π(1-ν)∫Ω-R¯¯,ij I[ε*]dx', x∉Ω$
(4a)
$σij*(x)=-μ4π(1-ν)∫Ω-R¯¯,ij I[ε*]dx'-2μεij*-2μ1-ννεkk*δij, x∈Ω$
(4b)
where $[ε*]=[ε11*; ε22*; ε33*; 2ε23*; 2ε13*; 2ε12*]T$, the semicolons are used to separate members in the vector; here the vector $R¯¯ I$ is defined as
$[R¯¯ I]=[R,11I; R,22I; R,33I; R,23I; R,13I; R,12I]$
(5)
where the potential $RI$ is introduced as follows:
$RI=|x-x'|$
(6)

Equation (4) is the general equation characterizing inclusions of arbitrary shape. When the domain Ω is a cuboid containing uniform eigenstrains, integration can be performed and explicit relations between the eigenstresses and eigenstrains obtained.

Here Eq. (4) is used to obtain the full space solutions. The closed-form cuboidal solutions can be found in Liu et al. [17]. The stresses due to an arbitrarily shaped inclusion can be written as a summation of its cuboidal contributions, i.e.,
$σij*(α,β,γ)=-μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(0) ×(α-ζ,β-η,γ-ϑ)εkl*(ζ,η,ϑ) (0≤α≤M-1,0≤β≤N-1,0≤γ≤L-1)$
(7)
where $(α,β,γ)$ or $(ζ,η,ϑ)$ is the central point of the cuboid. Equation (7), written in a matrix form, is as follows:
$[σ11*σ22*σ33*σ23*σ13*σ12*]= -μ4π(1-ν)[T1111(0)T1122(0)T1133(0)T1123(0)T1113(0)T1112(0)T2211(0)T2222(0)T2233(0)T2223(0)T2213(0)T2212(0)T3311(0)T3322(0)T3333(0)T3323(0)T3313(0)T3312(0)T2311(0)T2322(0)T2333(0)T2323(0)T2313(0)T2312(0)T1311(0)T1322(0)T1333(0)T1323(0)T1313(0)T1312(0)T1211(0)T1222(0)T1233(0)T1223(0)T1213(0)T1212(0)]*[ε11*ε22*ε33*2ε23*2ε13*2ε12*]$
(8)

where the asterisk symbolizes convolution; T is the influence coefficient matrix, whose explicit expression is listed in Appendix  A. It is worth mentioning that for a cuboidal inclusion with uniform eigenstrains inside, the internal stresses in Ω and external stresses outside Ω can be expressed by the same equation, analogous to the results of the two-dimensional problem [24,25].

After the full space solution is obtained, the half space problem is considered. Figure 2 shows the computational region, which is meshed into cuboidal elements of the same size, where the stresses at each element center and the surface tractions at the cuboidal surface center are determined. In addition, the surface tractions for each rectangular patch are treated as constant. The stress in a half space can be expressed as the contributions of each cuboidal inclusion as follows:

Fig. 2
Fig. 2
Close modal
$σij*(α,β,γ)=-μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(0)(α-ζ,β-η,γ-ϑ)εkl*(ζ,η,ϑ) +∑η=0N-1∑ζ=0M-1Cijτzz(α-ζ,β-η,γ)τzz(ζ,η) +∑η=0N-1∑ζ=0M-1Cijτzx(α-ζ,β-η,γ)τzx(ζ,η) +∑η=0N-1∑ζ=0M-1Cijτzy(α-ζ,β-η,γ)τzy(ζ,η) (0≤α≤M-1,0≤β≤N-1,0≤γ≤L-1)$
(9)
where $Cijτzz$, $Cijτzx$, $Cijτzy$ are influence coefficients relating the surface tractions $τzz$, $τzx$, and $τzy$ to stresses, where explicit forms can be found in Appendix  B. Similarly, Eq. (9), written in matrix form, is as follows:
$[σ11*σ22*σ33*σ23*σ13*σ12*]= -μ4π(1-ν)[T1111(0)T1122(0)T1133(0)T1123(0)T1113(0)T1112(0)T2211(0)T2222(0)T2233(0)T2223(0)T2213(0)T2212(0)T3311(0)T3322(0)T3333(0)T3323(0)T3313(0)T3312(0)T2311(0)T2322(0)T2333(0)T2323(0)T2313(0)T2312(0)T1311(0)T1322(0)T1333(0)T1323(0)T1313(0)T1312(0)T1211(0)T1222(0)T1233(0)T1223(0)T1213(0)T1212(0)]*[ε11*ε22*ε33*2ε23*2ε13*2ε12*]+[C11τzzC11τzxC11τzyC22τzzC22τzxC22τzyC33τzzC33τzxC33τzyC23τzzC23τzxC23τzyC13τzzC13τzxC13τzyC12τzzC12τzxC12τzy]*[τzzτzxτzy]$
(10)

The first term of Eq. (9) shows three-dimensional discrete convolutions of influence coefficients and eigenstrains, and the others are the two-dimensional discrete convolutions of influence coefficients and surface tractions. The discrete convolutions can be converted to circular discrete convolutions and the FFT can be employed to implement the computation. For the three-dimensional discrete convolution, when the size of the target domain is M × N × L, the extended computational domain must be 2M × 2N × 2L to avoid any FFT error. The calculation of influence coefficients is performed on the extended domain and implemented by the technique of wrap-around order while the eigenstrains in the extended domain need only zero padding [26]. The stresses can be obtained by taking an FFT of each sequence, i.e., influence coefficients and eigenstrains, multiplying point wise, and then performing an inverse FFT. Once completed, the stresses in the extended domain are discarded. Thus, the time to perform both FFT and inverse FFT in the first term of Eq. (9) is $O(864MNLln8MNL)$ and that for remaining terms is $O(216MNLln4MN)$.

Tractions at a surface point $(α,β)$, shown in Fig. 2, can be obtained by
$τzz(α,β)=μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1T33kl(0)(α-ζ,β-η,-ϑ)εkl*(ζ,η,ϑ)$
(11a)
$τzx(α,β)=μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1T13kl(0)(α-ζ,β-η,-ϑ)εkl*(ζ,η,ϑ)$
(11b)
$τzy(α,β)=μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1T23kl(0)(α-ζ,β-η,-ϑ)εkl*(ζ,η,ϑ)$
(11c)

The DC–FFT in the x and y directions can be employed to calculate the shear tractions; computing Eq. (11) requires 3 × L × 6 × 3 times 2D-FFT and inverse 2D-FFT. Thus the total time for performing FFT and inverse FFT in the computation for shear tractions is $O(216MNLln4MN)$. After obtaining the shear tractions, the stresses caused by shear tractions can be calculated by using the second, third, and fourth term of Eq. (9).

The time burden for performing the total FFT and inverse FFT in calculating stresses in a half space are listed in Table 1. The calculation of influence coefficients also requires considerable computing time. Because $Tijkl(0)$ is a symmetric matrix; only 21 terms of $Tijkl(0)$ need to be calculated in the new method.

Table 1

Time to perform FFT and inverse FFT in computing stresses in a half space with different methods

MethodsBased on Jacq et al. [8]Based on Zhou et al. [16]Based on Liu et al. [17]New method
Time$432MNL2ln4MN$$1728MNLln8MNL +144MNLln4MN$$3456MNLln8MNL$$864MNLln8MNL +432MNLln4MN$
MethodsBased on Jacq et al. [8]Based on Zhou et al. [16]Based on Liu et al. [17]New method
Time$432MNL2ln4MN$$1728MNLln8MNL +144MNLln4MN$$3456MNLln8MNL$$864MNLln8MNL +432MNLln4MN$

### Direct Method.

Recently, Liu et al. [17] derived the complete set of elastic field solutions caused by an eigenstrain and further obtained the influence coefficients to relate eigenstrain to displacement or stresses. Eigenstress $σij*$ due to eigenstrain in a half space can be expressed as [17]
$σij*(x)=-μ4π(1-ν)∫Ω(Θij-R¯¯,ij I-R¯¯,ij-2zR¯¯,3ij+2z2Φ¯¯,ij)[ε*]dx', x∉Ω$
(12a)
$σij*(x)= -μ4π(1-ν)∫Ω(Θij-R¯¯,ij I-R¯¯,ij-2zR¯¯,3ij+2z2Φ¯¯,ij)[ε*]dx'-2μεij*-2μ1-ννεkk*δij, x∈Ω$
(12b)
where vectors $Θ$, $R¯¯$, and $Φ¯¯$ are the combination of potentials. Similarly, an arbitrarily shaped inclusion can be divided into several small cuboidal inclusions with uniform eigenstrains inside of each one, where the total stresses are the superposition of the solutions caused by each cuboidal inclusion.
$σij*(α,β,γ)=-μ4π(1-ν)(∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(0) ×(α-ζ,β-η,γ-ϑ)εkl*(ζ,η,ϑ) +∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(1)(α-ζ,β-η,γ+ϑ)εkl*(ζ,η,ϑ) +z∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(2)(α-ζ,β-η,γ+ϑ)εkl*(ζ,η,ϑ) +z2∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(3)(α-ζ,β-η,γ+ϑ)εkl*(ζ,η,ϑ)) (0≤α≤M-1,0≤β≤N-1,0≤γ≤L-1)$
(13)

Here $Tijkl(0)$, $Tijkl(1)$, $Tijkl(2)$, and $Tijkl(3)$ are the influence coefficients relating eigenstrains to eigenstresses. More details on the influence coefficients $Tijkl(0)$, $Tijkl(1)$, $Tijkl(2)$, and $Tijkl(3)$ are available in Liu et al. [17], Appendix (A2). The first term of Eq. (13) includes 3D convolutions and the others include 3D convolution-correlation combinations. The convolution-correlation combination terms can be calculated by combining discrete convolution-discrete correlation and FFT (DC-DCR-FFT) (DC–FFT in the x and y directions and DCR–FFT in z direction). Likewise, 3D-FFT is employed to speed up the computation. Since the stress calculation needs 36 × 4 × 3 times 3D-FFT and inverse 3D-FFT, the total time for performing 3D-FFT is $O(3456MNLln8MNL)$ (see Table 1).

### Image Inclusion Approach.

The half space problem, shown in Fig. 3, can be solved by using the image method [5]. The cuboidal inclusion has uniform eigenstrains and its center is located at (0, 0, $z'$). Suppose the cuboidal inclusion has an image domain with its center located at (0, 0, −$z'$), and the eigenstrain components $εxx*$, $εyy*$, $εzz*$, and $εxy*$ in the two domains are equal while components $εxz*$, $εyz*$ are equal but opposite. Then, the plane of symmetry is free from shear tractions, and the normal stress in this plane due to the two domains is double the stress caused by a single domain. The stresses in a half space can be obtained by adding the stresses caused by two cuboidal inclusions in a full space and the stresses caused by surface normal traction, which can make the surface free of normal traction.

Fig. 3
Fig. 3
Close modal
An arbitrarily shaped inclusion is thus divided into several small cuboidal inclusions of the same size. The total stress caused by arbitrary inclusions can be calculated from the equation below [16],
$σij*(α,β,γ)=-μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(0) ×(α-ζ,β-η,γ-ϑ)εkl*(ζ,η,ϑ) +-μ4π(1-ν)∑ϑ=0L-1∑η=0N-1∑ζ=0M-1Tijkl(0)(α-ζ,β-η,γ+ϑ)εkl'*(ζ,η,ϑ) +2∑η=0N-1∑ζ=0M-1Cijτzz(α-ζ,β-η,γ)τzz(ζ,η) (0≤α≤M-1,0≤β≤N-1,0≤γ≤L-1)$
(14)

where $Tijkl(0)$ are the influence coefficients relating eigenstrain to eigenstress discussed above; $Cijτzz$ are the influence coefficients relating surface normal traction $τzz$ to stresses, and the explicit form of $Cijτzz$ can be found in Appendix  B. The image eigenstrain is $εij'*$, where $ε11'*=ε11*$, $ε22'*=ε22*$, $ε33'*=ε33*$, $ε12'*=ε12*$, $ε13'*=-ε13*$, and $ε23'*=-ε23*$. The surface normal traction $τzz$ can be obtained by using Eq. (11a). In Eq. (14), the first term is a 3D convolution group; the second term is a 3D convolution-correlation combination; and the third term is a two-dimensional convolution set.

The 3D-FFT can be employed to accelerate the calculation of first and second terms of Eq. (14) [15,16], and the 2D-FFT for the third term [16,27]. For the first and second terms in Eq. (14), the time burden to perform one set 3D-FFT or inverse 3D-FFT is $O(8MNLln8MNL)$. The total time consumption is listed in Table 1. In order to capture the surface normal traction $τzz$ accurately, the computational domain must be extended appropriately and covered by fine grids, which are discussed in Sec. 3 below.

## Numerical Results and Discussion

### Full Space Solutions.

First, the stresses caused by a cuboidal inclusion in a full space are considered. The stresses in this surface are influenced by the distance between the target surface and the inclusion. Figure 4 shows the stresses along the x axis caused by the cuboidal inclusion at different positions. The three sides of the cuboid are each of length 2a in the x, y, and z directions. The cuboid center is located at the z axis, and the xy plane is selected as the target plane, which is parallel to one of the cuboidal faces. Assuming uniform eigenstrains, $[ε*]=[εxx*;εyy*;εzz*;2εyz*;2εxz*;2εxy*]=10-3×[1,1,1,1,1,1]T$ inside and zero outside of the cuboidal inclusion. If the inclusion touches the target surface (Fig. 4(b)), the values of surface stresses $σzz$, $σzx$, and $σzy$ are large and dramatic variations at the left and right interfaces can be observed. As the inclusion is further from the target surface, the values of surface stresses decrease (Figs. 4(b)–4(d)). If three tractions $τzz=-σzz$, $τzx=-σzx$, and $τzy=-σzy$ are imposed on the target plane, the target plane becomes traction-free, and the half space solution can be obtained. In the new method, the stresses caused by tractions $τzz$, $τzx$, $τzy$ are calculated by using the numerical approach. Theoretically, tractions have nonzero values on the entire target plane; however, as they decrease drastically when the distance is further from the inclusions, the computational domain must be extended to reduce truncation errors caused by shear tractions. On the other hand, the computational domain is meshed by many cuboidal elements with uniform eigenstrains inside; thus, the mesh should be fine enough to represent arbitrarily shaped inclusions.

Fig. 4
Fig. 4
Close modal

### Half Space Solutions

#### A Single Inclusion.

In this section half space problems are considered, where the extreme cases, i.e., the inclusion touching the surface, are chosen. The inclusion size and the eigenstrains values inside are the same as those mentioned in Sec. 3.1. Six cases with different computational domains and grid sizes are listed in Table 2 and are used as benchmarks, where the analytical solution by Liu et al. [17] is employed.

Table 2

Cases of different computational domains and grid sizes

Case 1Case 2Case 3Case 4Case 5Case 6
Computational domain[−a,a] × [−a,a] × [0,3a][−3a,3a] × [−3a,3a] × [0,3a][−5a,5a] × [−5a,5a] × [0,3a][−3a,3a] × [−3a,3a] × [0,3a]
Cuboidal elements9 × 9 × 1527 × 27 × 1545 × 45 × 153 × 3 × 159 × 9 × 15243 × 243 × 15
Grid intervals(0.222, 0.222, 0.2)(2, 2, 0.2)(0.667, 0.667, 0.2)(0.025, 0.025, 0.2)
Case 1Case 2Case 3Case 4Case 5Case 6
Computational domain[−a,a] × [−a,a] × [0,3a][−3a,3a] × [−3a,3a] × [0,3a][−5a,5a] × [−5a,5a] × [0,3a][−3a,3a] × [−3a,3a] × [0,3a]
Cuboidal elements9 × 9 × 1527 × 27 × 1545 × 45 × 153 × 3 × 159 × 9 × 15243 × 243 × 15
Grid intervals(0.222, 0.222, 0.2)(2, 2, 0.2)(0.667, 0.667, 0.2)(0.025, 0.025, 0.2)
##### Effect of computational domain.

Cases 1 to 3 are analyzed, where the mesh size is held constant in the three directions, i.e., (Δx, Δy, Δz) = (0.222, 0.222, 0.2). Figure 5(a) shows the von Mises stress along the line passing through point (0, 0.1a, 0) and parallel to the x axis, and Fig. 5(b) shows the von Mises stress along the z axis. As the computational domain is increased, the results become more accurate and the computational domain of [−3a, 3a] × [−3a, 3a] × [0, 3a] in three directions is seen to be large enough to obtain a stable and accurate solution. Note that the computational domain is affected by the position of the inclusion. The extreme case occurs as the inclusion touches the surface, and the extended computational domain is larger in the extreme case than that for the inclusion inside the body.

Fig. 5
Fig. 5
Close modal
##### Effects of grid density.

A set of simulations on a fixed computational domain [−3a, 3a] × [−3a, 3a] × [0, 3a] with varying grid sizes is performed to study the effect of grid size. Figure 6 shows the von Mises stress along the two target lines, i.e., the line passing through point (0, 0.1a, 0) and parallel to the x axis and a line along the z axis. As the grid size increases, the results from the new method agree well with the analytical solutions. In case 6, the peaks of the stress are captured by using fine grids (Fig. 6(a)).

Fig. 6
Fig. 6
Close modal
##### Comparison of different methods.

The results obtained by the new method are compared to those from the method of Zhou et al. [16], where case 2 is taken as an example. Figure 7 shows the stresses $σxx$, $σyy$, and $σzz$ along the two target lines, similar to Figs. 5 and 6. The main source of error in the method of Zhou et al. [16] arises from the elimination of the surface normal traction; in the new approach, error is caused by elimination of both the surface normal and the shear tractions. Note that the normal traction in method of Zhou et al. [16] is twice that of the normal traction in the new method. In most cases, the method based on Zhou et al. [16] can achieve a better solution when the inclusions are near to the surface (Fig. 7(a)). When the computational domain is extended sufficiently large with grids dense enough along the x and y directions (as discussed above), the numerical results by using the two methods are almost the same as those obtained from the analytical solution.

Fig. 7
Fig. 7
Close modal

#### Multiple Inclusions.

A multi-inclusion problem is further considered, shown in Fig. 8(a), a half space having two layered inclusions with the top one touching the surface. Similarly, each inclusion size and the eigenstrains values inside the inclusions are the same as those mentioned in Sec. 3.1. The distance between two adjacent inclusions are all assumed to be a. A fixed computational domain [−6a, 6a] × [−6a, 6a] × [0, 6a] with 128 × 128 × 64 grids is used for the numerical solutions. The contours of the dimensionless von Mises stress in the xz plane are shown in Figs. 8(b)–8(d), where those in Figs. 8(b)–8(d) are the results from the new method, the methods of Zhou et al. [16] and Liu et al. [17], respectively. Multi-inclusion problem can be solved by using the three methods and the results are in good agreement with each other.

Fig. 8
Fig. 8
Close modal

### Elastoplastic Contact.

The contact of an elastoplastic sphere with elastic perfectly plastic behavior with a rigid plane is shown in Fig. 9 and is analyzed to evaluate the accuracy and efficiency of different methods. The loading conditions and the material properties of the sphere are listed in Table 3. In the numerical simulation, the loading path is divided into 10 steps. Here the plastic strains are treated as eigenstrains and the residual stresses as eigenstresses. The residual stresses due to plastic strains are calculated with an elastoplastic contact solver. In this section, the new method is compared with the image inclusion approach of Zhou et al. [16] and the direct method Liu et al. [17]. More details for the numerical approaches to solve the elastoplastic contact can be found in the related works [8–14].

Fig. 9
Fig. 9
Close modal
Table 3

Material parameters and contact conditions

ParameterCase 1
Normal load, W (N)800
Ball's radius, R (mm)20
Young's modulus of sphere, E (GPa)100
Poisson's ratios of sphere, ν0.3
Yield stress, σy (MPa)600
Maximum Hertzian contact pressure, ph (MPa)1671.92
Hertzian contact radius, a (mm)0.4780
Computational domain[−2a, 2a] × [−2a, 2a] × [0, 2a]
ParameterCase 1
Normal load, W (N)800
Ball's radius, R (mm)20
Young's modulus of sphere, E (GPa)100
Poisson's ratios of sphere, ν0.3
Yield stress, σy (MPa)600
Maximum Hertzian contact pressure, ph (MPa)1671.92
Hertzian contact radius, a (mm)0.4780
Computational domain[−2a, 2a] × [−2a, 2a] × [0, 2a]

Figure 10 plots the dimensionless pressure along the x axis and the equivalent plastic strain along the z axis, and Fig. 11 shows the contours of the dimensionless residual von Mises stress $σvonr/ph$ in the xz plane. More computational results comparing the three different methods are listed in Table 4. Excellent agreement from the three methods is observed. Figure 12 shows the computational time consumed by different methods using a computer of a 64-bit windows operating system with Intel® Core™ i7-920 Processor @ 2.67 GHz. As shown in Fig. 12, the computational speed of the new method is about twice faster than that of the image inclusion approach [16] and six times faster than the direct method [17].

Fig. 10
Fig. 10
Close modal
Fig. 11
Fig. 11
Close modal
Fig. 12
Fig. 12
Close modal
Table 4

Results for elastoplastic contact with different methods

New methodBased on Zhou et al. [16]Based on Liu et al. [17]
Max pressure, p/ph0.73290.73260.7322
Max equivalent plastic strain, $εeqp$0.00670.00670.0067
Max residual von Mises stress, $σvonr/ph$0.11570.11540.1151
Max von Mises stress, σvon/ph0.35890.35890.3589
New methodBased on Zhou et al. [16]Based on Liu et al. [17]
Max pressure, p/ph0.73290.73260.7322
Max equivalent plastic strain, $εeqp$0.00670.00670.0067
Max residual von Mises stress, $σvonr/ph$0.11570.11540.1151
Max von Mises stress, σvon/ph0.35890.35890.3589

## Parallel Computing Scheme

The computational efficiency can be significantly increased by using a parallel computing technique. There are two types of parallel programming models, i.e., open multiprocessing (OpenMP) and message passing interface (MPI). The former refers to computing on shared-memory computers with multiprocessing, while the latter is for distributed memory computers with multiple independent processors, and the data on each processor are inaccessible to other processors except by explicit message passing.

The parallel strategy for calculating eigenstress is shown in Fig. 13. The influence coefficients matrix T is a symmetric matrix for the full space problem. $Tijkl(0)$ is a block matrix, where every component of $Tijkl(0)$ is independent of each other, and can be calculated simultaneously through parallel computing. After obtained the influence coefficients, the discrete convolution of the eigenstrains and influence coefficients can employ the FFT algorithms via the circular convolution theorem. The eigenstresses can be obtained after performing 36 × 3 times FFT or inverse FFT. In the present study, a free software package, i.e., FFTW [28], is employed for performing the FFT. FFTW is a C subroutine library for computing the FFT in one or more dimensions of an arbitrary input size, and for both real and complex data (http://www.fftw.org). Parallel computing based on OpenMP is used when performing the FFT. The source code in this paper is written in standard Fortran 90 and compiled using Intel Intel® Visual Fortran Studio XE 2011 for Windows. The Intel® Core i7-920 Processor @ 2.67 GHz with four cores is used for speed comparison. In each of the examples below, all data types are in double precision.

Fig. 13
Fig. 13
Close modal

The efficiency of parallel FFTW is examined first. Figure 14 shows the execution time for different 3D transform size with the double-precision complex data type. To get more accurate results, each execution time is obtained by repeating 3D-FFT transforms 10 times. The results show that when the processor number is increased, the execution time decreases dramatically. With the same array size, the computational efficiency by using four processors is almost three times faster than a single processor. The execution time for calculating the eigenstresses are shown in Fig. 15 for different grid sizes. Furthermore, the elastoplastic case described in the Sec. 3 is analyzed, and the execution time is shown in Fig. 16. For a finer grid, i.e., 256 × 256 × 64, the time is 3295 s by using four processors and 10,043 s by one processor. Note that the calculation of eigenstrains, eigenstresses, or the residual stresses, requires roughly 95% of the computation time during the whole contact analysis. Several other portions of elastoplastic code, such as the calculation of influence coefficients relating pressure displacement and pressure stress, 2D-FFT involving stress and displacement are also performed similarly with the parallel implementation strategy. Parallel computing is well suited for the SAMs based on the influence coefficients and FFT.

Fig. 14
Fig. 14
Close modal
Fig. 15
Fig. 15
Close modal
Fig. 16
Fig. 16
Close modal

It is well worth mentioning that many researchers have studied the full space solutions for polyhedral inclusions with uniform eigenstrains [29–31]. The half space solution for tetrahedral or other types of inclusions can be calculated with a numerical method from the full space solution, as shown in Fig. 1. In general tetrahedral meshes are better for the arbitrary inclusions; however, FFT can only be employed for the cuboidal meshes of equal size (Fig. 2).

## Conclusions

A simple but efficient method is proposed to solve the stresses caused by inclusions in a half space, based on full space solutions. The verification calculation shows the solutions based on the new method agree well with the corresponding analytical solutions from the direct method of Liu [17]. The new method is also applicable to the extreme cases, i.e., when the inclusion is touching the surface. Here the computational domain along x and y directions must be sufficiently extended, e.g., three times larger than that of the inclusion domain, and the grid fine enough to discretize the inclusions. The computational speed based on the new method is about twice as fast as the imaging method of Zhou et al. [16] and six times as fast as the direct method by Liu et al. [17].

Furthermore, parallel computing is introduced to perform the FFT to obtain the influence coefficients. When solving for the eigenstresses in a half space, the computational efficiency is almost three times faster by using four processors than by using a single processor.

## Acknowledgment

Z. Wang would like to express sincere gratitude to the support from the National Science Foundation of China under Grant No. 51105391 and the Ph.D. Programs Foundation of Ministry of Education of China under Grant No. 20110191120009. The authors would like to acknowledge supports from Center for Surface Engineering and Tribology at Northwestern University, USA, and State Key Laboratory of Mechanical Transmission, Chongqing University, China. In addition, the authors like to thank Professor Nenzi Wang at Chang Gung University for inspiratory discussion.

### Nomenclature

Nomenclature

• a =

Hertzian contact radius, mm

•
• cn =

vectors linking the eight corners of the cuboidal inclusion to the target point

•
• $Cijkl$ =

elastic moduli, MPa

•
• $Cijτzz$, $Cijτzx$, $Cijτzy$ =

influence coefficients relating surface tractions $τzz$, $τzx$, and $τzy$ to stresses

•
• E =

Young's modulus, GPa

•
• p =

pressure, MPa

•
• ph =

maximum Hertzian pressure, MPa

•
• R =

radius of the sphere, mm

•
• $Tijkl(0)$, $Tijkl(1)$, $Tijkl(2)$, $Tijkl(3)$ =

influence coefficients relating eigenstrain to stress

•
• ui =

displacement, mm

•
• W =

applied normal load, N

•
• x, y, z =

space coordinates, mm

•
• δij =

Kronecker delta

•
• Δx, Δy, Δz =

grid size in the x, y, and z direction, mm

•
• $εij*$ =

eigenstrain

•
• ν =

Possion's ratio

•
• $σij$ =

stress, MPa

•
• $σij*$ =

eigenstress due to eigenstrain, MPa

•
• μ =

shear modulus, MPa

•
• $τzz$, $τzx$, $τzy$ =

surface tractions, MPa

### Special Symbols

Special Symbols

• * =

convolution

### Influence Coefficients Relating Eigenstrains to Eigenstresses in a Full Space

Consider a cuboidal inclusion with uniform eigenstrains $εij*$ in a full space. The center of the cuboid is located at the origin of the Cartesian system (x, y, z), and the three sides of the cuboid have length Δx, Δy, Δz in the x, y, and z directions, respectively. The eigenstress at the target point $Q(α,β,γ)$ can be calculated by
$σij*(α,β,γ)=-μ4π(1-ν)Tijkl(0)(α,β,γ)εkl*$
(A1)
where
$Tijkl(0)(α,β,γ)= -fijkl(x1,y1,z1)+fijkl(x2,y1,z1)-fijkl(x2,y2,z1)+fijkl(x1,y2,z1)-fijkl(x1,y2,z2)+fijkl(x1,y1,z2)-fijkl(x2,y1,z2)+fijkl(x2,y2,z2)$
(A2)
$x1=α-Δx/2, x2=α+Δx/2,y1=β-Δy/2,y2=β+Δy/2, z1=γ-Δz/2, z2=γ+Δz/2$
(A3)
and
$f1111(ξ1,ξ2,ξ3)=ξ1ξ3r(r+ξ2)+ξ1ξ2r(r+ξ3)+2[tan-1(ξ1ξ3rξ2)+tan-1(ξ1ξ2rξ3)]$
(A4)
$f1122(ξ1,ξ2,ξ3)=-ξ1ξ2r(r+ξ3)+2νtan-1(ξ1ξ2rξ3)$
(A5)
$f1133(ξ1,ξ2,ξ3)=-ξ1ξ3r(r+ξ2)+2νtan-1(ξ1ξ3rξ2)$
(A6)
$f1123(ξ1,ξ2,ξ3)=2νlog(r+ξ1)-ξ1r$
(A7)
$f1113(ξ1,ξ2,ξ3)=log(r+ξ2)-ξ12r(r+ξ2)$
(A8)
$f1112(ξ1,ξ2,ξ3)=log(r+ξ3)-ξ12r(r+ξ3)$
(A9)
$f2211(ξ1,ξ2,ξ3)=f1122(ξ1,ξ2,ξ3)$
(A10)
$f2222(ξ1,ξ2,ξ3)=ξ2ξ3r(r+ξ1)+ξ1ξ2r(r+ξ3)+2[tan-1(ξ2ξ3rξ1)+tan-1(ξ1ξ2rξ3)]$
(A11)
$f2233(ξ1,ξ2,ξ3)=-ξ2ξ3r(r+ξ1)+2νtan-1(ξ2ξ3rξ1)$
(A12)
$f2223(ξ1,ξ2,ξ3)=log(r+ξ1)-ξ22r(r+ξ1)$
(A13)
$f2213(ξ1,ξ2,ξ3)=2νlog(r+ξ2)-ξ2r$
(A14)
$f2212(ξ1,ξ2,ξ3)=log(r+ξ3)-ξ22r(r+ξ3)$
(A15)
$f3311(ξ1,ξ2,ξ3)=f1133(ξ1,ξ2,ξ3)$
(A16)
$f3322(ξ1,ξ2,ξ3)=f2233(ξ1,ξ2,ξ3)$
(A17)
$f3333(ξ1,ξ2,ξ3)=ξ1ξ3r(r+ξ2)+ξ2ξ3r(r+ξ1)+2[tan-1(ξ2ξ3rξ1)+tan-1(ξ1ξ3rξ2)]$
(A18)
$f3323(ξ1,ξ2,ξ3)=log(r+ξ1)-ξ32r(r+ξ1)$
(A19)
$f3313(ξ1,ξ2,ξ3)=log(r+ξ2)-ξ32r(r+ξ2)$
(A20)
$f3312(ξ1,ξ2,ξ3)=2νlog(r+ξ3)-ξ3r$
(A21)
$f2311(ξ1,ξ2,ξ3)=f1123(ξ1,ξ2,ξ3)$
(A22)
$f2322(ξ1,ξ2,ξ3)=f2223(ξ1,ξ2,ξ3)$
(A23)
$f2333(ξ1,ξ2,ξ3)=f3323(ξ1,ξ2,ξ3)$
(A24)
$f2323(ξ1,ξ2,ξ3)=-ξ2ξ3r(r+ξ1)+(1-ν)tan-1(ξ2ξ3rξ1)$
(A25)
$f2313(ξ1,ξ2,ξ3)=(1-ν)log(r+ξ3)-ξ3r$
(A26)
$f2312(ξ1,ξ2,ξ3)=(1-ν)log(r+ξ2)-ξ2r$
(A27)
$f1311(ξ1,ξ2,ξ3)=f1113(ξ1,ξ2,ξ3)$
(A28)
$f1322(ξ1,ξ2,ξ3)=f2213(ξ1,ξ2,ξ3)$
(A29)
$f1333(ξ1,ξ2,ξ3)=f3313(ξ1,ξ2,ξ3)$
(A30)
$f1323(ξ1,ξ2,ξ3)=f2313(ξ1,ξ2,ξ3)$
(A31)
$f1313(ξ1,ξ2,ξ3)=-ξ1ξ3r(r+ξ2)+(1-ν)tan-1(ξ1ξ3rξ2)$
(A32)
$f1312(ξ1,ξ2,ξ3)=(1-ν)log(r+ξ1)-ξ1r$
(A33)
$f1211(ξ1,ξ2,ξ3)=f1112(ξ1,ξ2,ξ3)$
(A34)
$f1222(ξ1,ξ2,ξ3)=f2212(ξ1,ξ2,ξ3)$
(A35)
$f1233(ξ1,ξ2,ξ3)=f3312(ξ1,ξ2,ξ3)$
(A36)
$f1223(ξ1,ξ2,ξ3)=f2312(ξ1,ξ2,ξ3)$
(A37)
$f1213(ξ1,ξ2,ξ3)=f1312(ξ1,ξ2,ξ3)$
(A38)
$f1212(ξ1,ξ2,ξ3)=-ξ1ξ2r(r+ξ3)+(1-ν)tan-1(ξ1ξ2rξ3)$
(A39)
where
$r=ξ12+ξ22+ξ32 and tan-1(ξ1ξ2rξ3)|ξ3=0=π2sign(ξ1ξ2), tan-1(ξ1ξ3rξ2)|ξ2=0=π2sign(ξ1ξ3), tan-1(ξ2ξ3rξ1)|ξ1=0=π2sign(ξ2ξ3)$
(A40)

### Influence Coefficients Relating Surface Tractions to Stresses in a Half Space

Consider uniform tractions $τzz$, $τzx$, and $τzy$ over a rectangular patch in a half space. The center of a rectangular patch is located at the origin of the Cartesian system (x, y, z), and each side of a rectangle has length Δx, Δy in the x and y directions, respectively. The stress at the target point $Q(α,β,γ)$ can be calculated by [27]
$σij(α,β,γ)=Cijτzz(α,β,γ)τzz+Cijτzx(α,β,γ)τzx+Cijτzy(α,β,γ)τzy$
(B1)
where
$Cijτ(α,β,γ)=gijτ(x1,y1,γ)+gijτ(x2,y2,γ)-gijτ(x1,y2,γ)-gijτ(x2,y1,γ)$
(B2)
$x1=α-Δx/2, x2=α+Δx/2, y1=β-Δy/2,y2=β+Δy/2$
(B3)
and
$g11τzz(ξ1,ξ2,ξ3)=12π[-2νtan-1(ξ1ξ2rξ3)+ 2(1-2ν)tan-1(ξ1r+ξ2+ξ3)-ξ1ξ3r(r+ξ2)]$
(B4)
$g22τzz(ξ1,ξ2,ξ3)=g11τzz(ξ2,ξ1,ξ3)$
(B5)
$g33τzz(ξ1,ξ2,ξ3)=12π[-tan-1ξ1ξ2rξ3+ξ1ξ3r(r+ξ2)+ξ2ξ3r(r+ξ1)]$
(B6)
$g23τzz(ξ1,ξ2,ξ3)=12π[-ξ32r(r+ξ1)]$
(B7)
$g13τzz(ξ1,ξ2,ξ3)=g23τzz(ξ2,ξ1,ξ3)$
(B8)
$g12τzz(ξ1,ξ2,ξ3)=12π[(2ν-1)log(r+ξ3)-ξ3r]$
(B9)
$g11τzx(ξ1,ξ2,ξ3)= 12π{2log(r+ξ2)+ξ3(1-2ν) ×[ξ2r(r+ξ3)+ξ3r(r+ξ2)]-2νξ12r(r+ξ2)}$
(B10)
$g22τzx(ξ1,ξ2,ξ3)=12π[2νlog(r+ξ2)-ξ3(1-2ν)ξ2r(r+ξ3)-2νξ2r]$
(B11)
$g33τzx(ξ1,ξ2,ξ3)=12π[-ξ32r(r+ξ2)]$
(B12)
$g23τzx(ξ1,ξ2,ξ3)=12π(-ξ3r)$
(B13)
$g13τzx(ξ1,ξ2,ξ3)=12π[-ξ1ξ3r(r+ξ2)-tan-1ξ1ξ2rξ3]$
(B14)
$g12τzx(ξ1,ξ2,ξ3)=12π[log(r+ξ1)-ξ3(1-2ν)ξ1r(r+ξ3)-2νξ1r]$
(B15)
$g11τzy(ξ1,ξ2,ξ3)=g22τzx(ξ2,ξ1,ξ3)$
(B16)
$g22τzy(ξ1,ξ2,ξ3)=g11τzx(ξ2,ξ1,ξ3)$
(B17)
$g33τzy(ξ1,ξ2,ξ3)=g33τzx(ξ2,ξ1,ξ3)$
(B18)
$g23τzy(ξ1,ξ2,ξ3)=g13τzx(ξ2,ξ1,ξ3)$
(B19)
$g13τzy(ξ1,ξ2,ξ3)=g23τzx(ξ2,ξ1,ξ3)$
(B20)
$g12τzy(ξ1,ξ2,ξ3)=g12τzx(ξ2,ξ1,ξ3)$
(B21)
where
$r=ξ12+ξ22+ξ32 and tan-1(ξ1ξ2rξ3)|ξ3=0=π2sign(ξ1ξ2)$
(B22)

## References

1.
Mura
,
T.
,
1993
,
Micromechanics of Defects in Solids
, 2nd and revised edition,
Kluwer Academic
,
Dordrecht, Netherlands
.
2.
Eshelby
,
J. D.
,
1957
, “
The Determination of the Elastic Field of an Ellipsoidal Inclusion, and Related Problems
,”
Proc. R. Soc. London Ser. A
,
241
, pp.
376
396
.10.1098/rspa.1957.0133
3.
Mindlin
,
R. D.
, and
Cheng
,
D. H.
,
1950
, “
Thermo Elastic Stress in the Semi-Infinite Solid
,”
J. Appl. Phys.
,
21
(
9
), pp.
931
933
.10.1063/1.1699786
4.
Seo
,
K.
, and
Mura
T.
,
1979
, “
The Elastic Field in a Half Space Due to Ellipsoidal Inclusions With Uniform Dilatational Eigenstrains
,”
ASME J. Appl. Mech.
,
46
(
3
), pp.
568
572
.10.1115/1.3424607
5.
Chiu
,
Y. P.
,
1978
, “
On the Stress Field and Surface Deformation in a Half Space With a Cuboidal Zone in Which Initial Strains Are Uniform
,”
ASME J. Appl. Mech.
,
45
(
2
), pp.
302
306
.10.1115/1.3424292
6.
Owen
,
D. R. J.
, and
Mura
,
T.
,
1967
, “
Periodic Dislocation Distributions in a Half-Space
,”
J. Appl. Phys.
,
38
(
5
), pp.
1999
2009
.10.1063/1.1709818
7.
Owen
,
D. R. J.
,
1971
, “
Solutions to Arbitrarily Oriented Periodic Dislocation and Eigenstrain Distributions in a Half-Space
,”
Int. J. Solids Struct.
,
7
(
10
), pp.
1343
1361
.10.1016/0020-7683(71)90050-3
8.
Jacq
,
C.
,
Nélias
,
D.
,
Lormand
,
G.
, and
Girodin
,
D.
,
2002
, “
Development of a Three-Dimensional Semi-Analytical Elastic–Plastic Contact Code
,”
ASME J. Tribol.
,
124
(
4
), pp.
653
667
.10.1115/1.1467920
9.
Boucly
,
V.
,
Nélias
,
D.
,
Liu
,
S.
,
Wang
,
Q.
, and
Keer
L. M.
,
2005
, “
Contact Analyses for Bodies With Frictional Heating and Plastic Behavior
,”
ASME J. Tribol.
,
127
(
2
), pp.
355
364
.10.1115/1.1843851
10.
Nélias
,
D.
,
Boucly
,
V.
, and
Brunet
M.
,
2006
, “
Elastic-Plastic Contact Between Rough Surfaces: Proposal for a Wear or Running-In Model
,”
ASME J. Tribol.
,
128
(
2
), pp.
236
244
.10.1115/1.2163360
11.
Nélias
,
D.
,
Antaluca
,
E.
,
Boucly
,
V.
, and
Cretu
,
S.
,
2007
, “
A Three-Dimensional Semi-Analytical Model for Elastic–Plastic Sliding Contacts
,”
ASME J. Tribol.
,
129
(
4
), pp.
761
771
.10.1115/1.2768076
12.
Chen
,
W. W.
, and
Wang
,
Q.
,
2008
, “
Thermomechanical Analysis of Elasto-Plastic Bodies in a Sliding Spherical Contact and the Effects of Sliding Speed, Heat Partition, and Thermal Softening
,”
ASME J. Tribol.
,
130
(
4
), p.
041402
.10.1115/1.2959110
13.
Chen
,
W. W.
,
Liu
,
S. B.
, and
Wang
,
Q.
,
2008
, “
Fast Fourier Transform Based Numerical Methods for Elasto-Plastic Contacts With Nominally Flat Surface
,”
ASME J. Appl. Mech.
,
75
(
1
), p.
011022
.10.1115/1.2755158
14.
Wang
,
Z. J.
,
Wang
,
W. Z.
,
Hu
,
Y. Z.
, and
Wang
H.
,
2010
, “
A Numerical Elastic-Plastic Contact Model for Rough Surfaces
,”
Tribol. Trans.
,
53
(
2
), pp.
224
238
.10.1080/10402000903177908
15.
Liu
,
S. B.
, and
Wang
,
Q.
,
2005
, “
Elastic Fields Due to Eigenstrains in a Half-Space
,”
ASME J. Appl. Mech.
,
72
(
6
), pp.
871
878
.10.1115/1.2047598
16.
Zhou
,
K.
,
Chen
,
W. W.
,
Keer
,
L. M.
, and
Wang
,
Q. J.
,
2009
, “
A Fast Method for Solving Three-Dimensional Arbitrarily Shaped Inclusions in a Half Space
,”
Comput. Methods Appl. Mech. Eng.
,
198
(
9–12
), pp.
885
892
.10.1016/j.cma.2008.10.021
17.
Liu
,
S. B.
,
Jin
,
X. Q.
,
Wang
,
Z. J.
,
Keer
,
L. M.
, and
Wang
,
Q.
,
2012
, “
Analytical Solution for Elastic Fields Caused by Eigenstrains in a Half-Space and Numerical Implementation Based on FFT
,”
Int. J. Plasticity
,
35
, pp.
135
154
.10.1016/j.ijplas.2012.03.002
18.
Wang
,
Z. J.
,
Jin
,
X. Q.
,
Keer
,
L. M.
, and
Wang
,
Q.
,
2012
, “
Numerical Methods for Contact Between Two Joined Quarter Spaces and a Rigid Sphere
,”
Int. J. Solids Struct.
,
49
(
18
), pp.
2515
2527
.10.1016/j.ijsolstr.2012.05.027
19.
Wang
,
Z. Q.
,
Ghoniem
,
N. M.
,
Swaminarayan
,
S.
, and
LeSar
,
R.
,
2006
, “
A Parallel Algorithm for 3D Dislocation Dynamics
,”
J. Comput. Phys.
,
219
(
2
), pp.
608
621
.10.1016/j.jcp.2006.04.005
20.
Chung
,
S. W.
,
Lee
,
C. S.
, and
Kim
,
S. K.
,
2006
, “
Large-Scale Simulation of Crack Propagation Based on Continuum Damage Mechanics and Two-Step Mesh Partitioning
,”
Mech. Mater.
,
38
(
1–2
), pp.
76
87
.10.1016/j.mechmat.2005.05.011
21.
Lee
,
H.
,
Gillman
,
A. S.
, and
Matouîs
,
K.
,
2011
, “
Computing Overall Elastic Constants of Polydisperse Particulate Composites From Microtomographic Data
,”
J. Mech. Phys. Solids.
,
59
(
9
), pp.
1838
1857
.10.1016/j.jmps.2011.05.010
22.
Amdahl
,
G. M.
,
1967
, “
Validity of the Single Processor Approach to Achieving Large-Scale Computing Capabilities
,”
AFIPS Conf. Proc.
,
30
, pp.
483
485
.10.1145/1465482.1465560
23.
Chiu
,
Y. P.
,
1977
, “
Stress-Field Due to Initial Strains in a Cuboid Surrounded by an Infinite Elastic Space
,”
ASME J. Appl. Mech.
,
44
(
4
), pp.
587
590
.10.1115/1.3424140
24.
Chiu
,
Y. P.
,
1980
, “
On the Internal Stresses in a Half-Plane and a Layer Containing Localized Inelastic Strains or Inclusions
,”
ASME J. Appl. Mech.
,
47
(
2
), pp.
313
318
.10.1115/1.3153661
25.
Jin
,
X.
,
Keer
,
L. M.
, and
Wang
,
Q.
,
2009
, “
New Green's Function for Stress Field and a Note of Its Application in Quantum-Wire Structures
,”
Int. J. Solids Struct.
,
46
(
21
), pp.
3788
3798
.10.1016/j.ijsolstr.2009.07.005
26.
Liu
,
S. B.
,
Wang
,
Q.
, and
Liu
,
G.
,
2000
, “
A Versatile Method of Discrete Convolution and FFT (DC-FFT) for Contact Analyses
,”
Wear
,
243
(
1–2
), pp.
101
111
.10.1016/S0043-1648(00)00427-0
27.
Liu
,
S. B.
, and
Wang
,
Q.
,
2002
, “
Study Contact Stress Fields Caused by Surface Tractions With a Discrete Convolution and Fast Fourier Transform Algorithm
,”
ASME J. Tribol.
,
124
(
1
), pp.
36
45
.10.1115/1.1401017
28.
FFTW Package, available at http://www.fftw.org/
29.
Rodin
,
G. J.
,
1996
, “
Eshelby's Inclusion Problem for Polygons and Polyhedral
,”
J. Mech. Phys. Solids
,
44
(
12
), pp.
1977
1995
.10.1016/S0022-5096(96)00066-X
30.
Nozaki
,
H.
, and
Taya
,
M.
,
2001
, “
Elastic Fields in a Polyhedral Inclusion With Uniform Eigenstrains and Related Problems
,”
ASME J. Appl. Mech.
,
68
(
3
), pp.
441
452
.10.1115/1.1362670
31.
Gao
,
X. L.
, and
Liu
,
M. Q.
,
2012
, “
Strain Gradient Solution for the Eshelby-Type Polyhedral Inclusion Problem
,”
J. Mech. Phys. Solids
,
60
(
2
), pp.
261
276
.10.1016/j.jmps.2011.10.010