An Efficient Numerical Method With a Parallel Computational Strategy for Solving Arbitrarily Shaped Inclusions in Elastoplastic Contact Problems

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.

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].

where the cuboidal inclusion is in domain Ω.

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.

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:

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)$⁠.

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.

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).

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.

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.

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.

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.

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.

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)).

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.

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.

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].

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].

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.

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.

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).

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.

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

a =

Hertzian contact radius, mm

c_n =

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

p_h =

maximum Hertzian pressure, MPa

R =

radius of the sphere, mm

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

influence coefficients relating eigenstrain to stress

u_i =

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

* =

convolution